/* Constants*/
#define constant WAD = 0xde0b6b3a7640000
#define constant HALF_WAD = 0x6f05b59d3b20000
#define constant TYPE_UINT_MAX =  0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff


#define macro WAD_MUL() = takes (2) returns (1) {
    // Input Stack:             // [b, a]
    [HALF_WAD]                  // [HALF_WAD, b, a]
    [TYPE_UINT_MAX]             // [UINT_MAX, HALF_WAD, b, a]
    sub                         // [UINT_MAX-HALF_WAD, b, a]
    dup2                        // [b,UINT_MAX-HALF_WAD,b, a]    
    dup2                        // [UINT_MAX-HALF_WAD,b,UINT_MAX-HALF_WAD,b, a]   
    div                         // [(UINT_MAX-HALF_WAD / b),UINT_MAX-HALF_WAD,b, a]
    dup4                        // [a, (UINT_MAX-HALF_WAD / b),UINT_MAX-HALF_WAD,b, a]
    swap1                       // [(UINT_MAX-HALF_WAD / b),a,UINT_MAX-HALF_WAD,b, a]
    gt return_result jumpi      // [(UINT_MAX-HALF_WAD / b) > a,UINT_MAX-HALF_WAD,b, a]
    0x00 0x00 revert
   
    
    return_result:
        pop                     // [b, a]  
        mul                     // [b * a]  
        [HALF_WAD]              // [HALF_WAD, b*a]
        add                     // [HALF_WAD + (b*a)]
        [WAD]                   // [1e18, HALF_WAD + (b*a)]
        swap1                   // [HALF_WAD + (b*a),1e18]
        div                     // [HALF_WAD + (b*a) / 1e18]
        // Return Stack         // [HALF_WAD + (b*a) / 1e18]
}

//https://github.com/pentagonxyz/huffmate/blob/main/src/math/FixedPointMath.huff
#define macro SQRT() = takes (1) returns (1) {
    // Input stack:                [x]

    // We start y at x, which will help us make our initial estimate.
    dup1                        // [y, x]

    // The "correct" value is 1, but this saves a multiplication later.
    0xb5                        // [0xb5, y, x]

    0x10000000000000000000000000000000000
    dup3                        // [y, 0x1000..., z, y, x]
    lt                          // [y < 0x1000..., z, y, x]
    continue_1 jumpi            // [z, y, x]
    
    0x40 shl                    // [z << 0x40, y, x]
    swap1                       // [y, z << 0x40, x]
    0x80 shr                    // [y >> 0x80, z << 0x40, x]
    swap1                       // [z << 0x40, y >> 0x80, x]

    continue_1:

    0x1000000000000000000
    dup3                        // [y, 0x1000..., z, y, x]
    lt                          // [y < 0x1000..., z, y, x]
    continue_2 jumpi            // [z, y, x]
    
    0x20 shl                    // [z << 0x20, y, x]
    swap1                       // [y, z << 0x20, x]
    0x40 shr                    // [y >> 0x40, z << 0x20, x]
    swap1                       // [z << 0x20, y >> 0x40, x]

    continue_2:

    0x10000000000
    dup3                        // [y, 0x1000..., z, y, x]
    lt                          // [y < 0x1000..., z, y, x]
    continue_3 jumpi            // [z, y, x]
    
    0x10 shl                    // [z << 0x10, y, x]
    swap1                       // [y, z << 0x10, x]
    0x20 shr                    // [y >> 0x20, z << 0x10, x]
    swap1                       // [z << 0x10, y >> 0x20, x]

    continue_3:

    0x1000000
    dup3                        // [y, 0x1000..., z, y, x]
    lt                          // [y < 0x1000..., z, y, x]
    continue_4 jumpi            // [z, y, x]
    
    0x08 shl                    // [z << 0x08, y, x]
    swap1                       // [y, z << 0x08, x]
    0x10 shr                    // [y >> 0x10, z << 0x08, x]
    swap1                       // [z << 0x08, y >> 0x10, x]

    continue_4:

    // Goal was to get z*z*y within a small factor of x. More iterations could
    // get y in a tighter range. Currently, we will have y in [256, 256*2^16).
    // We ensured y >= 256 so that the relative difference between y and y+1 is small.
    // That's not possible if x < 256 but we can just verify those cases exhaustively.

    // Now, z*z*y <= x < z*z*(y+1), and y <= 2^(16+8), and either y >= 256, or x < 256.
    // Correctness can be checked exhaustively for x < 256, so we assume y >= 256.
    // Then z*sqrt(y) is within sqrt(257)/sqrt(256) of sqrt(x), or about 20bps.

    // For s in the range [1/256, 256], the estimate f(s) = (181/1024) * (s+1) is in the range
    // (1/2.84 * sqrt(s), 2.84 * sqrt(s)), with largest error when s = 1 and when s = 256 or 1/256.

    // Since y is in [256, 256*2^16), let a = y/65536, so that a is in [1/256, 256). Then we can estimate
    // sqrt(y) using sqrt(65536) * 181/1024 * (a + 1) = 181/4 * (y + 65536)/65536 = 181 * (y + 65536)/2^18.

    // There is no overflow risk here since y < 2^136 after the first branch above.

    // A mul is saved from starting z at 181.
    swap1                       // [y, z, x]
    0x010000 add                // [y + 0x010000, z, x]
    mul                         // [(y + 0x010000) * z, x]
    0x12 shr                    // [((y + 0x010000) * z) >> 0x12, x]

    // Given the worst case multiplicative error of 2.84 above, 7 iterations should be enough.
    dup1 dup3                   // [x, z, z, x]
    div                         // [x / z, z, x]
    add                         // [x / z + z, x]
    0x01 shr                    // [(x / z + z) >> 0x01, x]

    dup1 dup3                   // [x, z, z, x]
    div                         // [x / z, z, x]
    add                         // [x / z + z, x]
    0x01 shr                    // [(x / z + z) >> 0x01, x]

    dup1 dup3                   // [x, z, z, x]
    div                         // [x / z, z, x]
    add                         // [x / z + z, x]
    0x01 shr                    // [(x / z + z) >> 0x01, x]

    dup1 dup3                   // [x, z, z, x]
    div                         // [x / z, z, x]
    add                         // [x / z + z, x]
    0x01 shr                    // [(x / z + z) >> 0x01, x]

    dup1 dup3                   // [x, z, z, x]
    div                         // [x / z, z, x]
    add                         // [x / z + z, x]
    0x01 shr                    // [(x / z + z) >> 0x01, x]

    dup1 dup3                   // [x, z, z, x]
    div                         // [x / z, z, x]
    add                         // [x / z + z, x]
    0x01 shr                    // [(x / z + z) >> 0x01, x]

    dup1 dup3                   // [x, z, z, x]
    div                         // [x / z, z, x]
    add                         // [x / z + z, x]
    0x01 shr                    // [(x / z + z) >> 0x01, x]

    // If x+1 is a perfect square, the Babylonian method cycles between
    // floor(sqrt(x)) and ceil(sqrt(x)). This statement ensures we return floor.
    // See: https://en.wikipedia.org/wiki/Integer_square_root#Using_only_integer_division
    // Since the ceil is rare, we save gas on the assignment and repeat division in the rare case.
    // If you don't care whether the floor or ceil square root is returned, you can remove this statement.
    dup1 dup1 swap3             // [x, z, z, z]
    div                         // [x / z, z, z]
    lt                          // [(x / z) < z, z]
    swap1 sub                   // [z - (x / z) < z]
}

