Biguint division (#159)

* Division draft

* Non working draft

* Fix compile errors

* Use proper pointers for quotient and remainder

* Add fix note, some more changes

* Add comment

* Implement `big_uint_bit_size`

* Increase pointer to prevent it from steping over console_log

* WIP divrem

* Fix loop and zero initializer

* Push test cases

* Add other test case

* Add docs and tests for `big_uint_inplace_or_1`

* Fix bug related to bit shifting

* Fix borrow return in big uint sub function

* Delete playground file used for debugging

* Fix sub with borrow function

* Add playground again to check more big integer division tests

* Remove playground used for testing

* Write documentation for new shift functions

* Improve naming and documentation for new helper functions

* Rename bigUIntOrWith1 to bigUintInPlaceOrWith1

* Add tmp buffer parameters to bigUIntDivRem. Improve docs.

* Simplify subLimbsWithBorrow

Co-authored-by: Ivan Litteri <[email protected]>

* Remove `mul` call from `bigUIntInPlaceOrWith1`

Co-authored-by: Ivan Litteri <[email protected]>

* Remove multiplications from copyBigUint

Co-authored-by: Ivan Litteri <[email protected]>

* Optimize bigUIntBitSize loop

Co-authored-by: Ivan Litteri <[email protected]>

* Simplify zeroWithLimbSizeAt

---------

Co-authored-by: Francisco Krause Arnim <[email protected]>
Co-authored-by: IAvecilla <[email protected]>
Co-authored-by: Ivan Litteri <[email protected]>

precompiles/Modexp.yul

@@ -11,9 +11,44 @@ object "ModExp" {
             function LIMB_SIZE_IN_BITS() -> limbSize {
                 limbSize := 0x100
             }
-
             // HELPER FUNCTIONS
 
+            /// @notice Stores a one in big unsigned integer form in memory.
+            /// @param nLimbs The number of limbs needed to represent the operand.
+            /// @param toAddress The pointer to the MSB of the destination.
+            function oneWithLimbSizeAt(nLimbs, toAddress) {
+                let pointerToOne :=  toAddress
+                mstore(pointerToOne, 0x1)
+                for { let i := sub(nLimbs, 1) } gt(i, 0) { i := sub(i, 1) } {
+                   let offset := add(mul(i, 32), pointerToOne)
+                   mstore(offset, 0x0)
+                }
+             }
+
+            /// @notice Stores a zero in big unsigned integer form in memory.
+            /// @param nLimbs The number of limbs needed to represent the operand.
+            /// @param toAddress The pointer to the MSB of the destination.
+            function zeroWithLimbSizeAt(nLimbs, toAddress) {
+                let overflow := add(toAddress, shl(5, nLimbs))
+                for { } lt(toAddress, overflow) { toAddress := add(toAddress, LIMB_SIZE_IN_BYTES()) } {
+                    mstore(toAddress, 0)
+                }
+            }
+
+            /// @notice Copy a big unsigned integer from one memory location to another.
+            /// @param nLimbs The number of limbs needed to represent the operand.
+            /// @param fromAddress The pointer to the MSB of the number to copy.
+            /// @param toAddress The pointer to the MSB of the destination.
+            function copyBigUint(nLimbs, fromAddress, toAddress) {
+                let offset := 0
+                for { let i := 0 } lt(i, nLimbs) { i := add(i, 1) } {
+                    let fromOffset := add(offset, fromAddress)
+                    let toOffset := add(offset, toAddress)
+                    mstore(toOffset, mload(fromOffset))
+                    offset := add(offset, LIMB_SIZE_IN_BYTES())
+                }
+            }
+
             /// @notice Computes an addition and checks for overflow.
             /// @param augend The value to add to.
             /// @param addend The value to add.
@@ -346,6 +381,19 @@ object "ModExp" {
                 mstore(add(limbPointer, anOffset), aValue)
             }
 
+            /// @notice Computes the difference between two 256 bit number and keeps
+            /// account of the borrow bit
+            /// in lshPointer and rhsPointer.
+            /// @dev Reference: https://github.com/lambdaclass/lambdaworks/blob/main/math/src/unsigned_integer/element.rs#L785
+            /// @param leftLimb The left side of the difference (i.e. the a in a - b).
+            /// @param rightLimb The right side of the difference (i.e. the b in a - b).
+            /// @return subtractionResult i.e. the c in c = a - b.
+            /// @return returnBorrow If there was any borrow on the subtraction, is returned as 1.
+            function subLimbsWithBorrow(leftLimb, rightLimb, limbBorrow) -> subtractionResult, returnBorrow {
+                let rightPlusBorrow := add(rightLimb, limbBorrow)
+                subtractionResult := sub(leftLimb, rightPlusBorrow)
+                returnBorrow := gt(subtractionResult, leftLimb)
+            }
             /// @notice Computes the BigUint subtraction between the number stored
             /// in minuendPtr and subtrahendPtr.
             /// @dev Reference: https://github.com/lambdaclass/lambdaworks/blob/main/math/src/unsigned_integer/element.rs#L795
@@ -357,12 +405,13 @@ object "ModExp" {
                 let minuendCurrentLimb
                 let subtrahendCurrentLimb
                 let differenceCurrentLimb
+                borrow := 0
                 let limbOffset := 0
                 for { let i := nLimbs } gt(i, 0) { i := sub(i, 1) } {
                     limbOffset := mul(sub(i,1), 32)
                     let minuendCurrentLimb := getLimbValueAtOffset(minuendPtr, limbOffset)
                     let subtrahendCurrentLimb := getLimbValueAtOffset(subtrahendPtr, limbOffset)
-                    let differenceCurrentLimb, borrow := overflowingSubWithBorrow(minuendCurrentLimb, subtrahendCurrentLimb, borrow)
+                    differenceCurrentLimb, borrow := overflowingSubWithBorrow(minuendCurrentLimb, subtrahendCurrentLimb, borrow)
                     storeLimbValueAtOffset(differencePtr, limbOffset, differenceCurrentLimb)
                 }
             }
@@ -402,6 +451,156 @@ object "ModExp" {
                 }
             }
 
+            // @notice Computes the bit size of an unsigned integer.
+            // @dev Return value boundary: `0 <= bitSize <= 256`
+            // @param number An unsigned integer value.
+            // @return bitSize Number of bits required to represent `number`.
+            function UIntBitSize(number) -> bitSize {
+                // Increment bitSize until there are no significant bits left.
+                bitSize := 0
+                for { let shift_me := number } lt(0, shift_me) { shift_me := shr(1, shift_me) } {
+                    bitSize := add(bitSize, 1)
+                }
+            }
+
+            /// @notice Computes the bit size of a big unsigned integer.
+            /// @param basePtr Base pointer for a big unsigned integer.
+            /// @param nLimbs The number of limbs needed to represent the operand.
+            /// @return bitSize Number of bits of the big unsigned integer.
+            function bigUIntBitSize(basePtr, nLimbs) -> bitSize {
+                bitSize := shl(8, nLimbs)
+
+                // Iterate until finding the most significant limb or reach the end of the limbs.
+                let limb := 0
+                let offset := 0
+                for { let i := 0 } and(lt(i, nLimbs), iszero(limb)) { i := add(i, 1) } {
+                    bitSize := sub(bitSize, 256) // Decrement one limb worth of bits.
+                    let ptr_i := add(basePtr, offset) // = basePtr + i * 32 bytes
+                    limb := mload(ptr_i)
+                    offset := add(offset, LIMB_SIZE_IN_BYTES())
+                }
+
+                // At this point, `limb == limbs[i - 1]`. Where `i` equals the
+                // last value it took.
+
+                // At this point, `bitSize` equals the amount of bits in the
+                // limbs following the most significant limb.
+
+                bitSize := add(bitSize, UIntBitSize(limb))
+            }
+
+            /// @notice Performs in-place `x | 1` operation.
+            /// @dev This function will mutate the memory space `mem[basePtr...(basePtr + nLimbs * 32)]`
+            /// @dev It consumes constant time, aka `O(1)`.
+            /// @param basePtr Base pointer for a big unsigned integer.
+            /// @param nLimbs Number of 32 Byte limbs composing the big unsigned integer.
+            function bigUIntInPlaceOrWith1(basePtr, nLimbs) {
+                let offset := shl(5, sub(nLimbs, 1))
+                let limbPtr := add(basePtr, offset)
+                let limb := mload(limbPtr)
+                mstore(limbPtr, or(limb, 0x1))
+            }
+
+            /// @notice Performs one shift to the left for a big unsigned integer (<<).
+            /// @dev The shift is performed in-place, mutating the memory space of the number.
+            /// @param numberPtr The pointer to the MSB of the number to shift.
+            /// @param nLimbs The number of limbs needed to represent the operand.
+            function bigUIntOneShiftLeft(numberPtr, nLimbs) {
+                let p := add(numberPtr, shl(5, nLimbs)) // numberPtr + 32 * nLimbs
+                let carryBit := 0
+                for {  } lt(numberPtr, p) {  } {
+                    p := sub(p, 32)
+                    let limb := mload(p)
+                    let msb := shr(255, limb)
+                    limb := or(shl(1, limb), carryBit)
+                    mstore(p, limb)
+                    carryBit := msb
+                }
+            }
+
+            /// @notice Performs one shift to the right for a big unsigned integer (>>).
+            /// @dev The shift is performed in-place, mutating the memory space of the number.
+            /// @param numberPtr The pointer to the MSB of the number to shift.
+            /// @param nLimbs The number of limbs needed to represent the operand.
+            function bigUIntOneShiftRight(numberPtr, nLimbs) {
+                let overflowPtr := add(numberPtr, shl(5, nLimbs))
+                let carryBit := 0
+                for { let p := numberPtr } lt(p, overflowPtr) { p := add(p, 32) } {
+                    let limb := mload(p)
+                    let lsb := and(limb, 1)
+                    limb := or(shr(1, limb), carryBit)
+                    carryBit := shl(255, lsb)
+                    mstore(p, limb)
+                }
+            }
+
+            /// @notice Computes the quotiend and reminder of dividing two big unsigned integers.
+            /// @dev
+            /// @dev Temporary buffers:
+            /// @dev ------------------
+            /// @dev
+            /// @dev This function requires two temporary buffers for internal storage:
+            /// @dev - Both buffers must provide `n_limbs * 32` bytes of writable memory space.
+            /// @dev - Neither buffer should overlap with each other.
+            /// @dev - Neither needs to be initialized to any particular value.
+            /// @dev - Consider the written values as undefined after the function returns.
+            /// @dev
+            /// @dev Return values:
+            /// @dev --------------
+            /// @dev
+            /// @dev - resulting `quotient` will be written `mem[base_ptr, base_ptr + 32 * n_limbs)`
+            /// @dev - resulting `reminder` will be written `mem[base_ptr, base_ptr + 32 * n_limbs)`
+            /// @dev
+            /// @param dividend_ptr Base pointer for a big unsigned integer representing the dividend.
+            /// @param divisor_ptr  Base pointer for a big unsigned integer representing the divisor.
+            /// @param tmp_ptr_1    Base pointer for a contiguous memory space of `n_limbs` for internal usage. Will be overwritten.
+            /// @param tmp_ptr_2    Base pointer for a contiguous memory space of `n_limbs` for internal usage. Will be overwritten.
+            /// @param n_limbs      Amount of limbs for each big unsigned integer.
+            /// @param quotient_ptr Base pointer for a big unsigned integer to write the division quotient.
+            /// @param rem_ptr Base pointer for a big unsigned integer to write the division remainder.
+            function bigUIntDivRem(dividend_ptr, divisor_ptr, tmp_ptr_1, tmp_ptr_2, n_limbs, quotient_ptr, rem_ptr) {
+                // Assign meaningful internal names to the temporary buffers passed as parameters. We use abstract names for
+                // parameters to prevent the leakage of implementation details.
+                let c_ptr := tmp_ptr_1
+                let r_ptr := tmp_ptr_2
+
+                copyBigUint(n_limbs, dividend_ptr, rem_ptr) // rem = dividend 
+
+                // Init quotient to 0.
+                zeroWithLimbSizeAt(n_limbs, quotient_ptr) // quotient = 0
+
+                let mb := bigUIntBitSize(divisor_ptr, n_limbs)
+                let bd := sub(mul(n_limbs, 256), mb)
+                bigUIntShl(bd, divisor_ptr, n_limbs, c_ptr) // c == divisor << bd
+
+                for { } iszero(0) { } {
+                    let borrow := bigUIntSubWithBorrow(rem_ptr, c_ptr, n_limbs, r_ptr)
+
+                    if iszero(borrow) {
+                        copyBigUint(n_limbs, r_ptr, rem_ptr)
+                    }
+
+                    copyBigUint(n_limbs, quotient_ptr, r_ptr) // r = quotient
+                    bigUIntInPlaceOrWith1(r_ptr, n_limbs) // r = quotient | 1
+
+                    if iszero(borrow) {
+                        copyBigUint(n_limbs, r_ptr, quotient_ptr)
+                    }
+
+                    if iszero(bd) {
+                        break
+                    }
+
+                    bd := sub(bd, 1)
+                    bigUIntOneShiftRight(c_ptr, n_limbs) // c = c >> 1
+                    bigUIntOneShiftLeft(quotient_ptr, n_limbs) // q[] = q[] << 1
+                }
+
+                if iszero(mb) {
+                    zeroWithLimbSizeAt(n_limbs, quotient_ptr)
+                }
+            }
+
             ////////////////////////////////////////////////////////////////
             //                      FALLBACK
             ////////////////////////////////////////////////////////////////