Float point arithmetic helpers (#21)

Co-authored-by: Amar Singh <asinghchrony@protonmail.com>

foundry.toml

@@ -18,3 +18,11 @@ gas_limit = 30000000
 gas_price = 1
 block_base_fee_per_gas = 0
 block_gas_limit = 30000000
+
+# Fuzz tests options
+[fuzz]
+# Reduce the numbers of runs if fuzz tests takes too long in your machine.
+runs = 2500
+
+# When debuging fuzz tests, uncomment this seed to make tests reproducible.
+# seed = "0xdeadbeefdeadbeefdeadbeefdeadbeef"

src/utils/BranchlessMath.sol

@@ -1,7 +1,20 @@
 // SPDX-License-Identifier: MIT
 // Analog's Contracts (last updated v0.1.0) (src/utils/BranchlessMath.sol)
 
-pragma solidity ^0.8.20;
+pragma solidity >=0.8.20;
+
+/**
+ * Rounding mode used when divide an integer.
+ */
+enum Rounding {
+    // Rounds towards zero
+    Floor,
+    // Rounds to the nearest value; if the number falls midway,
+    // it is rounded to the value above.
+    Nearest,
+    // Rounds towards positive infinite
+    Ceil
+}
 
 /**
  * @dev Utilities for branchless operations, useful when a constant gas cost is required.
@@ -38,6 +51,16 @@ library BranchlessMath {
         }
     }
 
+    /**
+     * @dev If `condition` is true returns `a`, otherwise returns `b`.
+     * see `BranchlessMath.ternary`
+     */
+    function ternary(bool condition, int256 a, int256 b) internal pure returns (int256 r) {
+        assembly {
+            r := xor(b, mul(xor(a, b), condition))
+        }
+    }
+
     /**
      * @dev If `condition` is true returns `a`, otherwise returns `b`.
      * see `BranchlessMath.ternary`
@@ -179,7 +202,7 @@ library BranchlessMath {
     }
 
     /**
-     * @dev Unsigned saturating left shift, bounds to UINT256 MAX instead of overflowing.
+     * @dev Unsigned saturating left shift, bounds to `2 ** 256 - 1` instead of overflowing.
      */
     function saturatingShl(uint256 x, uint8 shift) internal pure returns (uint256 r) {
         assembly {
@@ -191,6 +214,32 @@ library BranchlessMath {
         }
     }
 
+    /**
+     * @dev Returns the absolute unsigned value of a signed value.
+     *
+     * Ref: https://graphics.stanford.edu/~seander/bithacks.html#IntegerAbs
+     */
+    function abs(int256 a) internal pure returns (uint256 r) {
+        assembly {
+            // Formula from the "Bit Twiddling Hacks" by Sean Eron Anderson.
+            // Since `n` is a signed integer, the generated bytecode will use the SAR opcode to perform the right shift,
+            // taking advantage of the most significant (or "sign" bit) in two's complement representation.
+            // This opcode adds new most significant bits set to the value of the previous most significant bit. As a result,
+            // the mask will either be `bytes(0)` (if n is positive) or `~bytes32(0)` (if n is negative).
+            let mask := sar(255, a)
+
+            // A `bytes(0)` mask leaves the input unchanged, while a `~bytes32(0)` mask complements it.
+            r := xor(add(a, mask), mask)
+        }
+    }
+
+    /**
+     * @dev Computes the absolute difference between x and y.
+     */
+    function absDiff(uint256 x, uint256 y) internal pure returns (uint256) {
+        return abs(int256(x) - int256(y));
+    }
+
     /**
      * @dev Cast a boolean (false or true) to a uint256 (0 or 1) with no jump.
      */
@@ -207,4 +256,138 @@ library BranchlessMath {
     function toUint(address addr) internal pure returns (uint256) {
         return uint256(uint160(addr));
     }
+
+    /**
+     * @dev Count the consecutive zero bits (trailing) on the right.
+     */
+    function trailingZeros(uint256 x) internal pure returns (uint256 r) {
+        assembly {
+            // Compute largest power of two divisor of `x`.
+            x := and(x, sub(0, x))
+
+            // Use De Bruijn lookups to convert the power of 2 to log2(x).
+            // Reference: https://graphics.stanford.edu/~seander/bithacks.html#IntegerLogDeBruijn
+            r := byte(and(div(80, mod(x, 255)), 31), 0x0706050000040000010003000000000000000000020000000000000000000000)
+            r := add(byte(31, div(0xf8f0e8e0d8d0c8c0b8b0a8a09890888078706860585048403830282018100800, shr(r, x))), r)
+        }
+    }
+
+    /**
+     * @dev Return the log in base 2 of a positive value rounded towards zero.
+     * Returns 0 if given 0.
+     */
+    function log2(uint256 x) internal pure returns (uint256 r) {
+        unchecked {
+            // Round down to the closest power of 2
+            // Reference: https://graphics.stanford.edu/~seander/bithacks.html#RoundUpPowerOf2
+            x |= x >> 1;
+            x |= x >> 2;
+            x |= x >> 4;
+            x |= x >> 8;
+            x |= x >> 16;
+            x |= x >> 32;
+            x |= x >> 64;
+            x |= x >> 128;
+            x = (x >> 1) + 1;
+
+            // Use De Bruijn lookups to convert the power of 2 to floor(log2(x)).
+            // Reference: https://graphics.stanford.edu/~seander/bithacks.html#IntegerLogDeBruijn
+            assembly {
+                r :=
+                    byte(and(div(80, mod(x, 255)), 31), 0x0706050000040000010003000000000000000000020000000000000000000000)
+                r :=
+                    add(byte(31, div(0xf8f0e8e0d8d0c8c0b8b0a8a09890888078706860585048403830282018100800, shr(r, x))), r)
+            }
+        }
+    }
+
+    /**
+     * @notice Calculates x * y / denominator with full precision, following the selected rounding direction.
+     * Throws if result overflows a uint256 or denominator == 0.
+     *
+     * @dev This this an modified version of the original implementation by OpenZeppelin SDK, which is released under MIT.
+     * original: https://github.com/OpenZeppelin/openzeppelin-contracts/blob/v5.0.2/contracts/utils/math/Math.sol#L117-L202
+     */
+    function mulDiv(uint256 x, uint256 y, uint256 denominator, Rounding rounding) internal pure returns (uint256) {
+        unchecked {
+            // Compute remainder.
+            // - Rounding.Floor   then remainder is 0
+            // - Rounding.Nearest then remainder is denominator / 2
+            // - Rounding.Ceil    then remainder is denominator - 1
+            uint256 remainder = denominator;
+            remainder *= toUint(rounding != Rounding.Floor);
+            remainder >>= toUint(rounding == Rounding.Nearest);
+            remainder -= toUint(rounding == Rounding.Ceil);
+
+            // 512-bit multiply [prod1 prod0] = x * y + remainder.
+            // Compute the product mod 2²⁵⁶ and mod 2²⁵⁶ - 1, then use the Chinese Remainder Theorem to reconstruct
+            // the 512 bit result. The result is stored in two 256 variables such that product = prod1 * 2²⁵⁶ + prod0.
+            uint256 prod0 = x * y; // Least significant 256 bits of the product
+            uint256 prod1; // Most significant 256 bits of the product
+            assembly {
+                let mm := mulmod(x, y, not(0))
+                prod1 := sub(sub(mm, prod0), lt(mm, prod0))
+
+                // Add 256 bit remainder to 512 bit number.
+                mm := add(prod0, remainder)
+                prod1 := add(prod1, lt(mm, prod0))
+                prod0 := mm
+            }
+
+            // Make sure the result is less than 2**256. Also prevents denominator == 0.
+            require(prod1 < denominator, "muldiv overflow");
+
+            ///////////////////////////////////////////////
+            // 512 by 256 division.
+            ///////////////////////////////////////////////
+
+            // Make division exact by subtracting the remainder from [prod1 prod0].
+            assembly {
+                // Compute remainder using addmod and mulmod.
+                remainder := addmod(remainder, mulmod(x, y, denominator), denominator)
+
+                // Subtract 256 bit number from 512 bit number.
+                prod1 := sub(prod1, gt(remainder, prod0))
+                prod0 := sub(prod0, remainder)
+            }
+
+            // Factor powers of two out of denominator and compute largest power of two divisor of denominator.
+            // Always >= 1. See https://cs.stackexchange.com/q/138556/92363.
+
+            uint256 twos = denominator & (0 - denominator);
+            assembly {
+                // Divide denominator by twos.
+                denominator := div(denominator, twos)
+
+                // Divide [prod1 prod0] by twos.
+                prod0 := div(prod0, twos)
+
+                // Flip twos such that it is 2²⁵⁶ / twos. If twos is zero, then it becomes one.
+                twos := add(div(sub(0, twos), twos), 1)
+            }
+
+            // Shift in bits from prod1 into prod0.
+            prod0 |= prod1 * twos;
+
+            // Invert denominator mod 2²⁵⁶. Now that denominator is an odd number, it has an inverse modulo 2²⁵⁶ such
+            // that denominator * inv ≡ 1 mod 2²⁵⁶. Compute the inverse by starting with a seed that is correct for
+            // four bits. That is, denominator * inv ≡ 1 mod 2⁴.
+            uint256 inverse = (3 * denominator) ^ 2;
+
+            // Use the Newton-Raphson iteration to improve the precision. Thanks to Hensel's lifting lemma, this also
+            // works in modular arithmetic, doubling the correct bits in each step.
+            inverse *= 2 - denominator * inverse; // inverse mod 2⁸
+            inverse *= 2 - denominator * inverse; // inverse mod 2¹⁶
+            inverse *= 2 - denominator * inverse; // inverse mod 2³²
+            inverse *= 2 - denominator * inverse; // inverse mod 2⁶⁴
+            inverse *= 2 - denominator * inverse; // inverse mod 2¹²⁸
+            inverse *= 2 - denominator * inverse; // inverse mod 2²⁵⁶
+
+            // Because the division is now exact we can divide by multiplying with the modular inverse of denominator.
+            // This will give us the correct result modulo 2²⁵⁶. Since the preconditions guarantee that the outcome is
+            // less than 2²⁵⁶, this is the final result. We don't need to compute the high bits of the result and prod1
+            // is no longer required.
+            return prod0 * inverse;
+        }
+    }
 }