diff --git a/docs/utils/fixedpointmathlib.md b/docs/utils/fixedpointmathlib.md index 67d9d7ef3c..f06629fc3c 100644 --- a/docs/utils/fixedpointmathlib.md +++ b/docs/utils/fixedpointmathlib.md @@ -147,10 +147,7 @@ Equivalent to `(x * y) / WAD` rounded down. ### rawMulWad(uint256,uint256) ```solidity -function rawMulWad(uint256 x, uint256 y) - internal - pure - returns (uint256 z) +function rawMulWad(uint256 x, uint256 y) internal pure returns (uint256 z) ``` Equivalent to `(x * y) / WAD` rounded down, but without overflow checks. @@ -201,10 +198,7 @@ Equivalent to `(x * WAD) / y` rounded down. ### rawDivWad(uint256,uint256) ```solidity -function rawDivWad(uint256 x, uint256 y) - internal - pure - returns (uint256 z) +function rawDivWad(uint256 x, uint256 y) internal pure returns (uint256 z) ``` Equivalent to `(x * WAD) / y` rounded down, but without overflow and divide by zero checks. @@ -256,6 +250,22 @@ Returns `exp(x)`, denominated in `WAD`. Credit to Remco Bloemen under MIT license: https://2π.com/22/exp-ln Note: This function is an approximation. Monotonically increasing. +### expWadFast(int256) + +```solidity +function expWadFast(int256 x) internal pure returns (int256 r) +``` + +Returns `exp(x)`, denominated in `WAD`. Cheaper than `expWad`. +Let `E = exp(x / 1e18) * 1e18` denote the exact, infinite-precision result. +Never overestimates: the result is greater than `E * (1 - 2.58e-22) - 1` +and not more than `E`, so it is `floor(E)` or `floor(E) - 1` whenever +`x <= 8265113944572514620` (results up to `~3885 * 1e18`). +`expWadFast(0) = 1e18` exactly. +Monotonically increasing, including across all `2**k` seams. +The bounds are certified by exact-rational interval proofs, reproducible via +https://github.com/ddallaire/wad-exponentials + ### lnWad(int256) ```solidity @@ -266,10 +276,26 @@ Returns `ln(x)`, denominated in `WAD`. Credit to Remco Bloemen under MIT license: https://2π.com/22/exp-ln Note: This function is an approximation. Monotonically increasing. +### lnWadFast(int256) + +```solidity +function lnWadFast(int256 x) internal pure returns (int256 r) +``` + +Returns `ln(x)`, denominated in `WAD`. Cheaper than `lnWad`. +Let `L = ln(x / 1e18) * 1e18` denote the exact, infinite-precision result. +Never overestimates: always returns `floor(L)` or `floor(L) - 1`. +`lnWadFast(1e18) = 0` exactly. Monotonically increasing. +The bounds are certified by exact-rational interval proofs, reproducible via +https://github.com/ddallaire/wad-exponentials + ### lambertW0Wad(int256) ```solidity -function lambertW0Wad(int256 x) internal pure returns (int256 w) +function lambertW0Wad(int256 x) + internal + pure + returns (int256 w) ``` Returns `W_0(x)`, denominated in `WAD`. @@ -759,7 +785,7 @@ function lerp(uint256 a, uint256 b, uint256 t, uint256 begin, uint256 end) Returns `a + (b - a) * (t - begin) / (end - begin)`, with `t` clamped between `begin` and `end` (inclusive). Agnostic to the order of (`a`, `b`) and (`end`, `begin`). -If `begins == end`, returns `t <= begin ? a : b`. +If `begin == end`, returns `t <= begin ? a : b`. ### lerp(int256,int256,int256,int256,int256) @@ -773,7 +799,7 @@ function lerp(int256 a, int256 b, int256 t, int256 begin, int256 end) Returns `a + (b - a) * (t - begin) / (end - begin)`. with `t` clamped between `begin` and `end` (inclusive). Agnostic to the order of (`a`, `b`) and (`end`, `begin`). -If `begins == end`, returns `t <= begin ? a : b`. +If `begin == end`, returns `t <= begin ? a : b`. ### isEven(uint256) diff --git a/src/utils/FixedPointMathLib.sol b/src/utils/FixedPointMathLib.sol index 37103d59f3..5b81a7dd81 100644 --- a/src/utils/FixedPointMathLib.sol +++ b/src/utils/FixedPointMathLib.sol @@ -271,6 +271,84 @@ library FixedPointMathLib { } } + /// @dev Returns `exp(x)`, denominated in `WAD`. Cheaper than `expWad`. + /// Let `E = exp(x / 1e18) * 1e18` denote the exact, infinite-precision result. + /// Never overestimates: the result is greater than `E * (1 - 2.58e-22) - 1` + /// and not more than `E`, so it is `floor(E)` or `floor(E) - 1` whenever + /// `x <= 8265113944572514620` (results up to `~3885 * 1e18`). + /// `expWadFast(0) = 1e18` exactly. + /// Monotonically increasing, including across all `2**k` seams. + /// The bounds are certified by exact-rational interval proofs, reproducible via + /// https://github.com/ddallaire/wad-exponentials + function expWadFast(int256 x) internal pure returns (int256 r) { + unchecked { + // Accept `-41446531673892822313 < x < 135305999368893231589` with a + // single unsigned comparison; sort out the two edges on the cold path. + if (uint256(x) + 41446531673892822312 >= 176752531042786053901) { + // When the true result is less than 1 wei we return zero. + // This happens when `x <= (log(1e-18) * 1e18) ~ -4.15e19`. + if (x <= -41446531673892822313) return r; + + /// @solidity memory-safe-assembly + assembly { + // When the result is greater than `(2**255 - 1) / 1e18` we can not + // represent it as an int. This happens when + // `x >= floor(log((2**255 - 1) / 1e18) * 1e18) ≈ 135`. + mstore(0x00, 0xa37bfec9) // `ExpOverflow()`. + revert(0x1c, 0x04) + } + } + + // Convert `x` from `10**18` fixed point to `2**96` fixed point. + x = (x << 78) / 5 ** 18; + + // Reduce to `x' in (-½ ln 2, ½ ln 2) * 2**96` with `exp(x) = 2**k * exp(x')`. + // `6196328019 = round(2**128 / (ln 2 * 2**96))`; `k` is in the range `[-60, 195]`. + int256 k = (x * 6196328019 + 2 ** 127) >> 128; + x = x - k * 54916777467707473351141471128; + + // `exp(x') = (E + x' * O) / (E - x' * O)`, a symmetric rational with + // `E`, `O` polynomials in `x'^2`. `E` is monic: its last Horner stage + // needs no `>> 96`, so with the constant term pre-shifted, `e` and `t` + // are in `2**192` basis. The coefficients are jointly fitted so that the + // certified relative error stays one-sided after the margin below, and + // the error at the `+½ ln 2` edge is negative, which makes every seam + // between adjacent `2**k` octaves step upward (monotonicity). + int256 u = (x * x) >> 96; + int256 e = (((u + 66584530426202717196765975783591) * u) >> 96) + + 5993331421273161380223160234961546; + e = e * u + (52742053377336245150083666490725880 << 96); + int256 o = + ((3328678398953600544402144014937 * u) >> 96) + 799080153247570479545590910204849; + o = ((o * u) >> 96) + 26371026688668122575032340971836885; + int256 t = x * o; + + /// @solidity memory-safe-assembly + assembly { + // Div in assembly because solidity adds a zero check despite the unchecked. + // The denominator is positive on the whole reduced domain. + r := sdiv(add(e, t), sar(96, sub(e, t))) + } + + // Multiply by `2**k * 1e18 / 2**96`, less a margin that keeps the result + // at or below `E` at every certified error extreme. `r < 1.5 * 2**96`, + // so the product cannot overflow, and the shift amount is never negative. + r = int256( + (uint256(r) * 633825300114114700748270193445868131307960658286) >> uint256(195 - k) + ); + + /// @solidity memory-safe-assembly + assembly { + // `exp(0) = 1` is the only exact integer result in the domain; the + // margin lands it one unit low, so add one back exactly there. + // The reduced `x` is zero iff the input was zero (the truncated + // base conversion never lands on `k * 54916777467707473351141471128` + // for any other input; checked exhaustively for every `k`). + r := add(iszero(x), r) + } + } + } + /// @dev Returns `ln(x)`, denominated in `WAD`. /// Credit to Remco Bloemen under MIT license: https://2π.com/22/exp-ln /// Note: This function is an approximation. Monotonically increasing. @@ -346,6 +424,85 @@ library FixedPointMathLib { } } + /// @dev Returns `ln(x)`, denominated in `WAD`. Cheaper than `lnWad`. + /// Let `L = ln(x / 1e18) * 1e18` denote the exact, infinite-precision result. + /// Never overestimates: always returns `floor(L)` or `floor(L) - 1`. + /// `lnWadFast(1e18) = 0` exactly. Monotonically increasing. + /// The bounds are certified by exact-rational interval proofs, reproducible via + /// https://github.com/ddallaire/wad-exponentials + function lnWadFast(int256 x) internal pure returns (int256 r) { + /// @solidity memory-safe-assembly + assembly { + // Compute `k = log2(x) - 96`, `r = 159 - k = 255 - log2(x) = 255 ^ log2(x)`. + r := shl(7, lt(0xffffffffffffffffffffffffffffffff, x)) + r := or(r, shl(6, lt(0xffffffffffffffff, shr(r, x)))) + r := or(r, shl(5, lt(0xffffffff, shr(r, x)))) + r := or(r, shl(4, lt(0xffff, shr(r, x)))) + r := or(r, shl(3, lt(0xff, shr(r, x)))) + // We place the check here for more optimal stack operations. + if iszero(sgt(x, 0)) { + mstore(0x00, 0x1615e638) // `LnWadUndefined()`. + revert(0x1c, 0x04) + } + // forgefmt: disable-next-item + r := xor(r, byte(and(0x1f, shr(shr(r, x), 0x8421084210842108cc6318c6db6d54be)), + 0xf8f9f9faf9fdfafbf9fdfcfdfafbfcfef9fafdfafcfcfbfefafafcfbffffffff)) + + // Reduce range of x to (1, 2) * 2**96 + // ln(2^k * x) = k * ln(2) + ln(x) + x := shr(159, shl(r, x)) + + // `s = (x - sqrt(2)) * 2**96 / (x + sqrt(2))`, so that + // `ln(x) = ln(2)/2 + 2 * atanh(s)`. Centering on `sqrt(2)` halves the + // fit domain: `|s| <= 3 - 2 * sqrt(2)`. + let s := + sdiv( + shl(96, sub(x, 112045541949572279837463876455)), + add(x, 112045541949572279837463876455) + ) + + // `2 * atanh(s) = s * A(w) / B(w)`, a (3, 3)-term odd rational in `w = s^2`. + let w := sar(96, mul(s, s)) + let a := + add( + sar(96, mul(sub(w, 1813347344949966953757847210329), w)), + 5824670411451500986303020460168 + ) + a := sub(sar(96, mul(a, w)), 4518264490991587979207438354337) + let b := + sub( + sar(96, mul(188151507788160136135094921663, w)), + 1676640319226537252003611223372 + ) + b := add(sar(96, mul(b, w)), 3665379287557676720634158507137) + b := sub(sar(96, mul(b, w)), 2259132245495793985525851698055) + + // `B` is bounded away from zero on the whole domain. + let p := sdiv(mul(s, a), b) + + // Add `(2k + 1) * ln(2)/2` and `ln(2**96 / 10**18)`, then convert to `WAD`, + // all in `5**18 * 2**192` basis. The additive constant is lowered by a + // certified margin (~0.0417 wei) so the accumulator never exceeds + // `L * 2**174`; downward errors total under 1 wei, so `sar(174, p)` + // lands on `floor(L)` or `floor(L) - 1`. + p := mul(302231454903657293676544000000000000000000, p) + p := add( + mul( + 8298788776342807110743642979096973734596910279609939088954046749604186, + sub(319, shl(1, r)) + ), + p + ) + p := add(600920179829731861735705478226805774338444159116519230494951146088873754, p) + r := sar(174, p) + + // `ln(1e18 / 1e18) = 0` is the only exact integer result in the domain; + // the margin lands it exactly on `-1`, so add one back precisely there. + // (`floor(L) = -1` is unreachable: no input has `L in [-1, 0)`.) + r := add(iszero(not(r)), r) + } + } + /// @dev Returns `W_0(x)`, denominated in `WAD`. /// See: https://en.wikipedia.org/wiki/Lambert_W_function /// a.k.a. Product log function. This is an approximation of the principal branch. diff --git a/src/utils/clz/FixedPointMathLib.sol b/src/utils/clz/FixedPointMathLib.sol index 5b19001c84..57c54c83c5 100644 --- a/src/utils/clz/FixedPointMathLib.sol +++ b/src/utils/clz/FixedPointMathLib.sol @@ -271,6 +271,84 @@ library FixedPointMathLib { } } + /// @dev Returns `exp(x)`, denominated in `WAD`. Cheaper than `expWad`. + /// Let `E = exp(x / 1e18) * 1e18` denote the exact, infinite-precision result. + /// Never overestimates: the result is greater than `E * (1 - 2.58e-22) - 1` + /// and not more than `E`, so it is `floor(E)` or `floor(E) - 1` whenever + /// `x <= 8265113944572514620` (results up to `~3885 * 1e18`). + /// `expWadFast(0) = 1e18` exactly. + /// Monotonically increasing, including across all `2**k` seams. + /// The bounds are certified by exact-rational interval proofs, reproducible via + /// https://github.com/ddallaire/wad-exponentials + function expWadFast(int256 x) internal pure returns (int256 r) { + unchecked { + // Accept `-41446531673892822313 < x < 135305999368893231589` with a + // single unsigned comparison; sort out the two edges on the cold path. + if (uint256(x) + 41446531673892822312 >= 176752531042786053901) { + // When the true result is less than 1 wei we return zero. + // This happens when `x <= (log(1e-18) * 1e18) ~ -4.15e19`. + if (x <= -41446531673892822313) return r; + + /// @solidity memory-safe-assembly + assembly { + // When the result is greater than `(2**255 - 1) / 1e18` we can not + // represent it as an int. This happens when + // `x >= floor(log((2**255 - 1) / 1e18) * 1e18) ≈ 135`. + mstore(0x00, 0xa37bfec9) // `ExpOverflow()`. + revert(0x1c, 0x04) + } + } + + // Convert `x` from `10**18` fixed point to `2**96` fixed point. + x = (x << 78) / 5 ** 18; + + // Reduce to `x' in (-½ ln 2, ½ ln 2) * 2**96` with `exp(x) = 2**k * exp(x')`. + // `6196328019 = round(2**128 / (ln 2 * 2**96))`; `k` is in the range `[-60, 195]`. + int256 k = (x * 6196328019 + 2 ** 127) >> 128; + x = x - k * 54916777467707473351141471128; + + // `exp(x') = (E + x' * O) / (E - x' * O)`, a symmetric rational with + // `E`, `O` polynomials in `x'^2`. `E` is monic: its last Horner stage + // needs no `>> 96`, so with the constant term pre-shifted, `e` and `t` + // are in `2**192` basis. The coefficients are jointly fitted so that the + // certified relative error stays one-sided after the margin below, and + // the error at the `+½ ln 2` edge is negative, which makes every seam + // between adjacent `2**k` octaves step upward (monotonicity). + int256 u = (x * x) >> 96; + int256 e = (((u + 66584530426202717196765975783591) * u) >> 96) + + 5993331421273161380223160234961546; + e = e * u + (52742053377336245150083666490725880 << 96); + int256 o = + ((3328678398953600544402144014937 * u) >> 96) + 799080153247570479545590910204849; + o = ((o * u) >> 96) + 26371026688668122575032340971836885; + int256 t = x * o; + + /// @solidity memory-safe-assembly + assembly { + // Div in assembly because solidity adds a zero check despite the unchecked. + // The denominator is positive on the whole reduced domain. + r := sdiv(add(e, t), sar(96, sub(e, t))) + } + + // Multiply by `2**k * 1e18 / 2**96`, less a margin that keeps the result + // at or below `E` at every certified error extreme. `r < 1.5 * 2**96`, + // so the product cannot overflow, and the shift amount is never negative. + r = int256( + (uint256(r) * 633825300114114700748270193445868131307960658286) >> uint256(195 - k) + ); + + /// @solidity memory-safe-assembly + assembly { + // `exp(0) = 1` is the only exact integer result in the domain; the + // margin lands it one unit low, so add one back exactly there. + // The reduced `x` is zero iff the input was zero (the truncated + // base conversion never lands on `k * 54916777467707473351141471128` + // for any other input; checked exhaustively for every `k`). + r := add(iszero(x), r) + } + } + } + /// @dev Returns `ln(x)`, denominated in `WAD`. /// Credit to Remco Bloemen under MIT license: https://2π.com/22/exp-ln /// Note: This function is an approximation. Monotonically increasing. @@ -338,6 +416,79 @@ library FixedPointMathLib { } } + /// @dev Returns `ln(x)`, denominated in `WAD`. Cheaper than `lnWad`. + /// Let `L = ln(x / 1e18) * 1e18` denote the exact, infinite-precision result. + /// Never overestimates: always returns `floor(L)` or `floor(L) - 1`. + /// `lnWadFast(1e18) = 0` exactly. Monotonically increasing. + /// The bounds are certified by exact-rational interval proofs, reproducible via + /// https://github.com/ddallaire/wad-exponentials + function lnWadFast(int256 x) internal pure returns (int256 r) { + /// @solidity memory-safe-assembly + assembly { + if iszero(sgt(x, 0)) { + mstore(0x00, 0x1615e638) // `LnWadUndefined()`. + revert(0x1c, 0x04) + } + + // Compute `k = log2(x) - 96`, `r = 159 - k = 255 - log2(x) = 255 - (255 - clz(x)) + // then k = clz(x)`. + r := clz(x) + + // Reduce range of x to (1, 2) * 2**96 + // ln(2^k * x) = k * ln(2) + ln(x) + x := shr(159, shl(r, x)) + + // `s = (x - sqrt(2)) * 2**96 / (x + sqrt(2))`, so that + // `ln(x) = ln(2)/2 + 2 * atanh(s)`. Centering on `sqrt(2)` halves the + // fit domain: `|s| <= 3 - 2 * sqrt(2)`. + let s := + sdiv( + shl(96, sub(x, 112045541949572279837463876455)), + add(x, 112045541949572279837463876455) + ) + + // `2 * atanh(s) = s * A(w) / B(w)`, a (3, 3)-term odd rational in `w = s^2`. + let w := sar(96, mul(s, s)) + let a := + add( + sar(96, mul(sub(w, 1813347344949966953757847210329), w)), + 5824670411451500986303020460168 + ) + a := sub(sar(96, mul(a, w)), 4518264490991587979207438354337) + let b := + sub( + sar(96, mul(188151507788160136135094921663, w)), + 1676640319226537252003611223372 + ) + b := add(sar(96, mul(b, w)), 3665379287557676720634158507137) + b := sub(sar(96, mul(b, w)), 2259132245495793985525851698055) + + // `B` is bounded away from zero on the whole domain. + let p := sdiv(mul(s, a), b) + + // Add `(2k + 1) * ln(2)/2` and `ln(2**96 / 10**18)`, then convert to `WAD`, + // all in `5**18 * 2**192` basis. The additive constant is lowered by a + // certified margin (~0.0417 wei) so the accumulator never exceeds + // `L * 2**174`; downward errors total under 1 wei, so `sar(174, p)` + // lands on `floor(L)` or `floor(L) - 1`. + p := mul(302231454903657293676544000000000000000000, p) + p := add( + mul( + 8298788776342807110743642979096973734596910279609939088954046749604186, + sub(319, shl(1, r)) + ), + p + ) + p := add(600920179829731861735705478226805774338444159116519230494951146088873754, p) + r := sar(174, p) + + // `ln(1e18 / 1e18) = 0` is the only exact integer result in the domain; + // the margin lands it exactly on `-1`, so add one back precisely there. + // (`floor(L) = -1` is unreachable: no input has `L in [-1, 0)`.) + r := add(iszero(not(r)), r) + } + } + /// @dev Returns `W_0(x)`, denominated in `WAD`. /// See: https://en.wikipedia.org/wiki/Lambert_W_function /// a.k.a. Product log function. This is an approximation of the principal branch. diff --git a/test/FixedPointMathLib.t.sol b/test/FixedPointMathLib.t.sol index edc37f066c..c522bc8bae 100644 --- a/test/FixedPointMathLib.t.sol +++ b/test/FixedPointMathLib.t.sol @@ -806,6 +806,279 @@ contract FixedPointMathLibTest is SoladyTest { return FixedPointMathLib.lnWad(x); } + // The `expWadFast` / `lnWadFast` constants, margins, and golden vectors below + // are derived and certified (exact-rational one-sided interval proofs plus + // exhaustive seam checks) by the scripts in + // https://github.com/ddallaire/wad-exponentials + + function expWadFast(int256 x) public pure returns (int256) { + return FixedPointMathLib.expWadFast(x); + } + + function lnWadFast(int256 x) public pure returns (int256) { + return FixedPointMathLib.lnWadFast(x); + } + + function testExpWadFast() public { + // `[x, expWadFast(x), floor(E)]` with `E = exp(x / 1e18) * 1e18` exact. + // Certified: `expWadFast(x) <= E`, and for `x <= 8265113944572514620` + // the result is `floor(E)` or `floor(E) - 1`. + // forgefmt: disable-next-item + int256[3][28] memory v = [ + [int256(0), int256(1000000000000000000), int256(1000000000000000000)], + [int256(1), int256(1000000000000000000), int256(1000000000000000001)], + [int256(-1), int256(999999999999999998), int256(999999999999999999)], + [int256(2), int256(1000000000000000001), int256(1000000000000000002)], + [int256(-2), int256(999999999999999997), int256(999999999999999998)], + [int256(1000000000000000000), int256(2718281828459045235), int256(2718281828459045235)], + [int256(-1000000000000000000), int256(367879441171442321), int256(367879441171442321)], + [int256(2000000000000000000), int256(7389056098930650227), int256(7389056098930650227)], + [int256(-2000000000000000000), int256(135335283236612691), int256(135335283236612691)], + [int256(3000000000000000000), int256(20085536923187667740), int256(20085536923187667740)], + [int256(-3000000000000000000), int256(49787068367863942), int256(49787068367863942)], + [int256(500000000000000000), int256(1648721270700128146), int256(1648721270700128146)], + [int256(-500000000000000000), int256(606530659712633423), int256(606530659712633423)], + [int256(300000000000000000), int256(1349858807576003103), int256(1349858807576003103)], + [int256(-300000000000000000), int256(740818220681717866), int256(740818220681717866)], + [int256(693147180559945309), int256(1999999999999999999), int256(1999999999999999999)], + [int256(10000000000000000000), int256(22026465794806716516952), int256(22026465794806716516957)], + [int256(50000000000000000000), int256(5184705528587072464086860598170708410702), int256(5184705528587072464087453322933485384827)], + [int256(100000000000000000000), int256(26881171418161354484121294484800788523547843748099059358033104), int256(26881171418161354484126255515800135873611118773741922415191608)], + [int256(135305999368893231588), int256(57896044618658097649807840469722558706031793413539064376035740825232986243708), int256(57896044618658097649816762928942336782129491980154662247847962410455084893091)], + [int256(-41446531673892822312), int256(1), int256(1)], + [int256(-41446531673892822311), int256(1), int256(1)], + [int256(8265113944572514620), int256(3885915731585811120099), int256(3885915731585811120100)], + [int256(8265113944572514619), int256(3885915731585811116213), int256(3885915731585811116214)], + [int256(346573590264282617), int256(1414213562350905984), int256(1414213562350905984)], + [int256(346573590264282618), int256(1414213562350905986), int256(1414213562350905986)], + [int256(1234567890123456789), int256(3436893084346008004), int256(3436893084346008004)], + [int256(-9876543210987654321), int256(51365531132686), int256(51365531132686)] + ]; + unchecked { + for (uint256 i; i != 28; ++i) { + int256 r = FixedPointMathLib.expWadFast(v[i][0]); + assertEq(r, v[i][1]); + assertLe(r, v[i][2]); // Never overestimates `floor(E)`. + if (v[i][0] <= 8265113944572514620) assertGe(r, v[i][2] - 1); + } + } + assertEq(FixedPointMathLib.expWadFast(-41446531673892822313), 0); + assertEq(FixedPointMathLib.expWadFast(type(int256).min), 0); + } + + function testLnWadFast() public { + // `[x, lnWadFast(x), floor(L)]` with `L = ln(x / 1e18) * 1e18` exact. + // Certified: the result is always `floor(L)` or `floor(L) - 1`. + // forgefmt: disable-next-item + int256[3][30] memory v = [ + [int256(1), int256(-41446531673892822313), int256(-41446531673892822313)], + [int256(2), int256(-40753384493332877003), int256(-40753384493332877003)], + [int256(3), int256(-40347919385224712621), int256(-40347919385224712621)], + [int256(1000), int256(-34538776394910685261), int256(-34538776394910685261)], + [int256(1000000000), int256(-20723265836946411157), int256(-20723265836946411157)], + [int256(1000000000000000), int256(-6907755278982137053), int256(-6907755278982137053)], + [int256(100000000000000000), int256(-2302585092994045685), int256(-2302585092994045685)], + [int256(300000000000000000), int256(-1203972804325935993), int256(-1203972804325935993)], + [int256(700000000000000000), int256(-356674943938732379), int256(-356674943938732379)], + [int256(999999999999999999), int256(-2), int256(-2)], + [int256(1000000000000000000), int256(0), int256(0)], + [int256(1000000000000000001), int256(0), int256(0)], + [int256(1000000000000000002), int256(1), int256(1)], + [int256(1414213562373095048), int256(346573590279972654), int256(346573590279972654)], + [int256(2718281828459045235), int256(999999999999999999), int256(999999999999999999)], + [int256(2000000000000000000), int256(693147180559945309), int256(693147180559945309)], + [int256(5000000000000000000), int256(1609437912434100374), int256(1609437912434100374)], + [int256(10000000000000000000), int256(2302585092994045683), int256(2302585092994045684)], + [int256(31415926535897932384), int256(3447314978843445858), int256(3447314978843445858)], + [int256(1000000000000000000000), int256(6907755278982137052), int256(6907755278982137052)], + [int256(1000000000000000000000000), int256(13815510557964274104), int256(13815510557964274104)], + [int256(1000000000000000000000000000), int256(20723265836946411156), int256(20723265836946411156)], + [int256(1000000000000000000000000000000000000), int256(41446531673892822312), int256(41446531673892822312)], + [int256(999999999999999999999), int256(6907755278982137052), int256(6907755278982137052)], + [int256(79228162514264337593543950336), int256(25095597659861927391), int256(25095597659861927391)], + [int256(340282366920938463463374607431768211456), int256(47276307437780177293), int256(47276307437780177293)], + [int256(1461501637330902918203684832716283019655932542976), int256(69457017215698427194), int256(69457017215698427194)], + [int256(1606938044258990275541962092341162602522202993782792835301376), int256(97182904438096239571), int256(97182904438096239571)], + [int256(28948022309329048855892746252171976963317496166410141009864396001978282409984), int256(134612852188333286279), int256(134612852188333286279)], + [int256(57896044618658097711785492504343953926634992332820282019728792003956564819967), int256(135305999368893231588), int256(135305999368893231589)] + ]; + unchecked { + for (uint256 i; i != 30; ++i) { + int256 r = FixedPointMathLib.lnWadFast(v[i][0]); + assertEq(r, v[i][1]); + assertTrue(r == v[i][2] || r == v[i][2] - 1); + } + } + } + + function testExpWadFastOverflowReverts() public { + vm.expectRevert(FixedPointMathLib.ExpOverflow.selector); + this.expWadFast(135305999368893231589); + vm.expectRevert(FixedPointMathLib.ExpOverflow.selector); + this.expWadFast(type(int256).max); + assertGt(this.expWadFast(135305999368893231588), 0); + } + + function testLnWadFastNegativeReverts() public { + vm.expectRevert(FixedPointMathLib.LnWadUndefined.selector); + this.lnWadFast(-1); + vm.expectRevert(FixedPointMathLib.LnWadUndefined.selector); + this.lnWadFast(-2 ** 255); + vm.expectRevert(FixedPointMathLib.LnWadUndefined.selector); + this.lnWadFast(0); + } + + function testExpWadFastMonotonicallyIncreasing(int256 a, int256 b) public { + a = _boundExpWadFastInput(a); + b = _boundExpWadFastInput(b); + if (a > b) (a, b) = (b, a); + assertLe(FixedPointMathLib.expWadFast(a), FixedPointMathLib.expWadFast(b)); + } + + function testExpWadFastMonotonicallyIncreasingAround(int256 t) public { + t = _boundExpWadFastInput(t); + unchecked { + for (int256 x = t - 2; x != t + 2; ++x) { + assertLe( + FixedPointMathLib.expWadFast(_boundExpWadFastInput(x)), + FixedPointMathLib.expWadFast(_boundExpWadFastInput(x + 1)) + ); + } + } + } + + function testExpWadFastSeamsMonotonic() public { + // Exhaustively checks the adjacent-input pair at every `2**k` seam of the + // range reduction, where the reduced argument jumps across the full fitted + // domain. These are the only monotonicity-critical points; the certified + // per-wei accumulator gain dominates everywhere else. + unchecked { + for (int256 k = -59; k <= 195; ++k) { + int256 x = _expWadFastSeamOf(k); + assertLe(FixedPointMathLib.expWadFast(x - 1), FixedPointMathLib.expWadFast(x)); + } + } + } + + function testLnWadFastOctaveSeamsMonotonic() public { + // Exhaustively checks the adjacent-input pair at every power-of-two + // mantissa boundary. + unchecked { + for (uint256 n = 1; n != 255; ++n) { + assertLe( + FixedPointMathLib.lnWadFast(int256(1 << n) - 1), + FixedPointMathLib.lnWadFast(int256(1 << n)) + ); + } + } + } + + function testLnWadFastMonotonicallyIncreasing(uint256 a, uint256 b) public { + a = _bound(a, 1, 2 ** 255 - 1); + b = _bound(b, 1, 2 ** 255 - 1); + if (a > b) (a, b) = (b, a); + assertLe(FixedPointMathLib.lnWadFast(int256(a)), FixedPointMathLib.lnWadFast(int256(b))); + } + + function testExpWadFastLnWadFastRoundTrip(int256 x) public { + // Both functions never overestimate, so the round trips can never exceed + // the identity: `lnWadFast(expWadFast(x)) <= x` and vice versa. + x = _boundExpWadFastInput(x); + int256 e = FixedPointMathLib.expWadFast(x); + if (e > 0) assertLe(FixedPointMathLib.lnWadFast(e), x); + } + + function testLnWadFastExpWadFastRoundTrip(uint256 x_) public { + int256 x = int256(_bound(x_, 1, 2 ** 255 - 1)); + int256 l = FixedPointMathLib.lnWadFast(x); + assertLe(FixedPointMathLib.expWadFast(l), x); + } + + function testLnWadFastDifferential(uint256 x_) public { + int256 x = int256(_bound(x_, 1, 2 ** 255 - 1)); + assertLe( + FixedPointMathLib.dist(FixedPointMathLib.lnWadFast(x), FixedPointMathLib.lnWad(x)), 2 + ); + } + + function testExpWadFastDifferential(int256 x) public { + x = _boundExpWadFastInput(x); + uint256 fast = uint256(FixedPointMathLib.expWadFast(x)); + uint256 orig = uint256(FixedPointMathLib.expWad(x)); + // `expWadFast` is certified within `2.58e-22` relative one-sided; + // `expWad`'s measured error is below `1e-17` relative. + assertLe(FixedPointMathLib.dist(fast, orig), orig / 1e17 + 2); + } + + function testExpWadGas() public { + unchecked { + uint256 acc; + for (int256 x = -40e18; x <= 130e18; x += 10e18) { + acc ^= uint256(FixedPointMathLib.expWad(x)); + } + assertTrue(acc != 0); + } + } + + function testExpWadFastGas() public { + unchecked { + uint256 acc; + for (int256 x = -40e18; x <= 130e18; x += 10e18) { + acc ^= uint256(FixedPointMathLib.expWadFast(x)); + } + assertTrue(acc != 0); + } + } + + function testLnWadGas() public { + unchecked { + uint256 acc; + for (int256 x = 1e9; x <= 1e63; x *= 1e6) { + acc ^= uint256(FixedPointMathLib.lnWad(x)); + } + assertTrue(acc != 0); + } + } + + function testLnWadFastGas() public { + unchecked { + uint256 acc; + for (int256 x = 1e9; x <= 1e63; x *= 1e6) { + acc ^= uint256(FixedPointMathLib.lnWadFast(x)); + } + assertTrue(acc != 0); + } + } + + function _boundExpWadFastInput(int256 x) private pure returns (int256) { + unchecked { + // Map into `[-41446531673892822312, 135305999368893231588]`. + return -41446531673892822312 + int256(uint256(x) % 176752531042786053901); + } + } + + function _expWadFastSeamOf(int256 k) private pure returns (int256 x) { + // Smallest `x` whose range reduction lands on octave `k`, by binary + // search over the exact `k(x)` map (monotone nondecreasing). + unchecked { + int256 lo = -41446531673892822312; + int256 hi = 135305999368893231588; + while (lo < hi) { + int256 mid = (lo + hi) >> 1; + if (_expWadFastKOf(mid) >= k) hi = mid; + else lo = mid + 1; + } + assert(_expWadFastKOf(lo) == k && _expWadFastKOf(lo - 1) == k - 1); + return lo; + } + } + + function _expWadFastKOf(int256 x) private pure returns (int256) { + unchecked { + return ((x << 78) / 5 ** 18 * 6196328019 + 2 ** 127) >> 128; + } + } + function testRPow() public { assertEq(FixedPointMathLib.rpow(0, 0, 0), 0); assertEq(FixedPointMathLib.rpow(1, 0, 0), 0); diff --git a/test/clz/FixedPointMathLib.t.sol b/test/clz/FixedPointMathLib.t.sol index 5eeca929ad..6b70a693b3 100644 --- a/test/clz/FixedPointMathLib.t.sol +++ b/test/clz/FixedPointMathLib.t.sol @@ -806,6 +806,279 @@ contract FixedPointMathLibWithCLZTest is SoladyTest { return FixedPointMathLib.lnWad(x); } + // The `expWadFast` / `lnWadFast` constants, margins, and golden vectors below + // are derived and certified (exact-rational one-sided interval proofs plus + // exhaustive seam checks) by the scripts in + // https://github.com/ddallaire/wad-exponentials + + function expWadFast(int256 x) public pure returns (int256) { + return FixedPointMathLib.expWadFast(x); + } + + function lnWadFast(int256 x) public pure returns (int256) { + return FixedPointMathLib.lnWadFast(x); + } + + function testExpWadFast() public { + // `[x, expWadFast(x), floor(E)]` with `E = exp(x / 1e18) * 1e18` exact. + // Certified: `expWadFast(x) <= E`, and for `x <= 8265113944572514620` + // the result is `floor(E)` or `floor(E) - 1`. + // forgefmt: disable-next-item + int256[3][28] memory v = [ + [int256(0), int256(1000000000000000000), int256(1000000000000000000)], + [int256(1), int256(1000000000000000000), int256(1000000000000000001)], + [int256(-1), int256(999999999999999998), int256(999999999999999999)], + [int256(2), int256(1000000000000000001), int256(1000000000000000002)], + [int256(-2), int256(999999999999999997), int256(999999999999999998)], + [int256(1000000000000000000), int256(2718281828459045235), int256(2718281828459045235)], + [int256(-1000000000000000000), int256(367879441171442321), int256(367879441171442321)], + [int256(2000000000000000000), int256(7389056098930650227), int256(7389056098930650227)], + [int256(-2000000000000000000), int256(135335283236612691), int256(135335283236612691)], + [int256(3000000000000000000), int256(20085536923187667740), int256(20085536923187667740)], + [int256(-3000000000000000000), int256(49787068367863942), int256(49787068367863942)], + [int256(500000000000000000), int256(1648721270700128146), int256(1648721270700128146)], + [int256(-500000000000000000), int256(606530659712633423), int256(606530659712633423)], + [int256(300000000000000000), int256(1349858807576003103), int256(1349858807576003103)], + [int256(-300000000000000000), int256(740818220681717866), int256(740818220681717866)], + [int256(693147180559945309), int256(1999999999999999999), int256(1999999999999999999)], + [int256(10000000000000000000), int256(22026465794806716516952), int256(22026465794806716516957)], + [int256(50000000000000000000), int256(5184705528587072464086860598170708410702), int256(5184705528587072464087453322933485384827)], + [int256(100000000000000000000), int256(26881171418161354484121294484800788523547843748099059358033104), int256(26881171418161354484126255515800135873611118773741922415191608)], + [int256(135305999368893231588), int256(57896044618658097649807840469722558706031793413539064376035740825232986243708), int256(57896044618658097649816762928942336782129491980154662247847962410455084893091)], + [int256(-41446531673892822312), int256(1), int256(1)], + [int256(-41446531673892822311), int256(1), int256(1)], + [int256(8265113944572514620), int256(3885915731585811120099), int256(3885915731585811120100)], + [int256(8265113944572514619), int256(3885915731585811116213), int256(3885915731585811116214)], + [int256(346573590264282617), int256(1414213562350905984), int256(1414213562350905984)], + [int256(346573590264282618), int256(1414213562350905986), int256(1414213562350905986)], + [int256(1234567890123456789), int256(3436893084346008004), int256(3436893084346008004)], + [int256(-9876543210987654321), int256(51365531132686), int256(51365531132686)] + ]; + unchecked { + for (uint256 i; i != 28; ++i) { + int256 r = FixedPointMathLib.expWadFast(v[i][0]); + assertEq(r, v[i][1]); + assertLe(r, v[i][2]); // Never overestimates `floor(E)`. + if (v[i][0] <= 8265113944572514620) assertGe(r, v[i][2] - 1); + } + } + assertEq(FixedPointMathLib.expWadFast(-41446531673892822313), 0); + assertEq(FixedPointMathLib.expWadFast(type(int256).min), 0); + } + + function testLnWadFast() public { + // `[x, lnWadFast(x), floor(L)]` with `L = ln(x / 1e18) * 1e18` exact. + // Certified: the result is always `floor(L)` or `floor(L) - 1`. + // forgefmt: disable-next-item + int256[3][30] memory v = [ + [int256(1), int256(-41446531673892822313), int256(-41446531673892822313)], + [int256(2), int256(-40753384493332877003), int256(-40753384493332877003)], + [int256(3), int256(-40347919385224712621), int256(-40347919385224712621)], + [int256(1000), int256(-34538776394910685261), int256(-34538776394910685261)], + [int256(1000000000), int256(-20723265836946411157), int256(-20723265836946411157)], + [int256(1000000000000000), int256(-6907755278982137053), int256(-6907755278982137053)], + [int256(100000000000000000), int256(-2302585092994045685), int256(-2302585092994045685)], + [int256(300000000000000000), int256(-1203972804325935993), int256(-1203972804325935993)], + [int256(700000000000000000), int256(-356674943938732379), int256(-356674943938732379)], + [int256(999999999999999999), int256(-2), int256(-2)], + [int256(1000000000000000000), int256(0), int256(0)], + [int256(1000000000000000001), int256(0), int256(0)], + [int256(1000000000000000002), int256(1), int256(1)], + [int256(1414213562373095048), int256(346573590279972654), int256(346573590279972654)], + [int256(2718281828459045235), int256(999999999999999999), int256(999999999999999999)], + [int256(2000000000000000000), int256(693147180559945309), int256(693147180559945309)], + [int256(5000000000000000000), int256(1609437912434100374), int256(1609437912434100374)], + [int256(10000000000000000000), int256(2302585092994045683), int256(2302585092994045684)], + [int256(31415926535897932384), int256(3447314978843445858), int256(3447314978843445858)], + [int256(1000000000000000000000), int256(6907755278982137052), int256(6907755278982137052)], + [int256(1000000000000000000000000), int256(13815510557964274104), int256(13815510557964274104)], + [int256(1000000000000000000000000000), int256(20723265836946411156), int256(20723265836946411156)], + [int256(1000000000000000000000000000000000000), int256(41446531673892822312), int256(41446531673892822312)], + [int256(999999999999999999999), int256(6907755278982137052), int256(6907755278982137052)], + [int256(79228162514264337593543950336), int256(25095597659861927391), int256(25095597659861927391)], + [int256(340282366920938463463374607431768211456), int256(47276307437780177293), int256(47276307437780177293)], + [int256(1461501637330902918203684832716283019655932542976), int256(69457017215698427194), int256(69457017215698427194)], + [int256(1606938044258990275541962092341162602522202993782792835301376), int256(97182904438096239571), int256(97182904438096239571)], + [int256(28948022309329048855892746252171976963317496166410141009864396001978282409984), int256(134612852188333286279), int256(134612852188333286279)], + [int256(57896044618658097711785492504343953926634992332820282019728792003956564819967), int256(135305999368893231588), int256(135305999368893231589)] + ]; + unchecked { + for (uint256 i; i != 30; ++i) { + int256 r = FixedPointMathLib.lnWadFast(v[i][0]); + assertEq(r, v[i][1]); + assertTrue(r == v[i][2] || r == v[i][2] - 1); + } + } + } + + function testExpWadFastOverflowReverts() public { + vm.expectRevert(FixedPointMathLib.ExpOverflow.selector); + this.expWadFast(135305999368893231589); + vm.expectRevert(FixedPointMathLib.ExpOverflow.selector); + this.expWadFast(type(int256).max); + assertGt(this.expWadFast(135305999368893231588), 0); + } + + function testLnWadFastNegativeReverts() public { + vm.expectRevert(FixedPointMathLib.LnWadUndefined.selector); + this.lnWadFast(-1); + vm.expectRevert(FixedPointMathLib.LnWadUndefined.selector); + this.lnWadFast(-2 ** 255); + vm.expectRevert(FixedPointMathLib.LnWadUndefined.selector); + this.lnWadFast(0); + } + + function testExpWadFastMonotonicallyIncreasing(int256 a, int256 b) public { + a = _boundExpWadFastInput(a); + b = _boundExpWadFastInput(b); + if (a > b) (a, b) = (b, a); + assertLe(FixedPointMathLib.expWadFast(a), FixedPointMathLib.expWadFast(b)); + } + + function testExpWadFastMonotonicallyIncreasingAround(int256 t) public { + t = _boundExpWadFastInput(t); + unchecked { + for (int256 x = t - 2; x != t + 2; ++x) { + assertLe( + FixedPointMathLib.expWadFast(_boundExpWadFastInput(x)), + FixedPointMathLib.expWadFast(_boundExpWadFastInput(x + 1)) + ); + } + } + } + + function testExpWadFastSeamsMonotonic() public { + // Exhaustively checks the adjacent-input pair at every `2**k` seam of the + // range reduction, where the reduced argument jumps across the full fitted + // domain. These are the only monotonicity-critical points; the certified + // per-wei accumulator gain dominates everywhere else. + unchecked { + for (int256 k = -59; k <= 195; ++k) { + int256 x = _expWadFastSeamOf(k); + assertLe(FixedPointMathLib.expWadFast(x - 1), FixedPointMathLib.expWadFast(x)); + } + } + } + + function testLnWadFastOctaveSeamsMonotonic() public { + // Exhaustively checks the adjacent-input pair at every power-of-two + // mantissa boundary. + unchecked { + for (uint256 n = 1; n != 255; ++n) { + assertLe( + FixedPointMathLib.lnWadFast(int256(1 << n) - 1), + FixedPointMathLib.lnWadFast(int256(1 << n)) + ); + } + } + } + + function testLnWadFastMonotonicallyIncreasing(uint256 a, uint256 b) public { + a = _bound(a, 1, 2 ** 255 - 1); + b = _bound(b, 1, 2 ** 255 - 1); + if (a > b) (a, b) = (b, a); + assertLe(FixedPointMathLib.lnWadFast(int256(a)), FixedPointMathLib.lnWadFast(int256(b))); + } + + function testExpWadFastLnWadFastRoundTrip(int256 x) public { + // Both functions never overestimate, so the round trips can never exceed + // the identity: `lnWadFast(expWadFast(x)) <= x` and vice versa. + x = _boundExpWadFastInput(x); + int256 e = FixedPointMathLib.expWadFast(x); + if (e > 0) assertLe(FixedPointMathLib.lnWadFast(e), x); + } + + function testLnWadFastExpWadFastRoundTrip(uint256 x_) public { + int256 x = int256(_bound(x_, 1, 2 ** 255 - 1)); + int256 l = FixedPointMathLib.lnWadFast(x); + assertLe(FixedPointMathLib.expWadFast(l), x); + } + + function testLnWadFastDifferential(uint256 x_) public { + int256 x = int256(_bound(x_, 1, 2 ** 255 - 1)); + assertLe( + FixedPointMathLib.dist(FixedPointMathLib.lnWadFast(x), FixedPointMathLib.lnWad(x)), 2 + ); + } + + function testExpWadFastDifferential(int256 x) public { + x = _boundExpWadFastInput(x); + uint256 fast = uint256(FixedPointMathLib.expWadFast(x)); + uint256 orig = uint256(FixedPointMathLib.expWad(x)); + // `expWadFast` is certified within `2.58e-22` relative one-sided; + // `expWad`'s measured error is below `1e-17` relative. + assertLe(FixedPointMathLib.dist(fast, orig), orig / 1e17 + 2); + } + + function testExpWadGas() public { + unchecked { + uint256 acc; + for (int256 x = -40e18; x <= 130e18; x += 10e18) { + acc ^= uint256(FixedPointMathLib.expWad(x)); + } + assertTrue(acc != 0); + } + } + + function testExpWadFastGas() public { + unchecked { + uint256 acc; + for (int256 x = -40e18; x <= 130e18; x += 10e18) { + acc ^= uint256(FixedPointMathLib.expWadFast(x)); + } + assertTrue(acc != 0); + } + } + + function testLnWadGas() public { + unchecked { + uint256 acc; + for (int256 x = 1e9; x <= 1e63; x *= 1e6) { + acc ^= uint256(FixedPointMathLib.lnWad(x)); + } + assertTrue(acc != 0); + } + } + + function testLnWadFastGas() public { + unchecked { + uint256 acc; + for (int256 x = 1e9; x <= 1e63; x *= 1e6) { + acc ^= uint256(FixedPointMathLib.lnWadFast(x)); + } + assertTrue(acc != 0); + } + } + + function _boundExpWadFastInput(int256 x) private pure returns (int256) { + unchecked { + // Map into `[-41446531673892822312, 135305999368893231588]`. + return -41446531673892822312 + int256(uint256(x) % 176752531042786053901); + } + } + + function _expWadFastSeamOf(int256 k) private pure returns (int256 x) { + // Smallest `x` whose range reduction lands on octave `k`, by binary + // search over the exact `k(x)` map (monotone nondecreasing). + unchecked { + int256 lo = -41446531673892822312; + int256 hi = 135305999368893231588; + while (lo < hi) { + int256 mid = (lo + hi) >> 1; + if (_expWadFastKOf(mid) >= k) hi = mid; + else lo = mid + 1; + } + assert(_expWadFastKOf(lo) == k && _expWadFastKOf(lo - 1) == k - 1); + return lo; + } + } + + function _expWadFastKOf(int256 x) private pure returns (int256) { + unchecked { + return ((x << 78) / 5 ** 18 * 6196328019 + 2 ** 127) >> 128; + } + } + function testRPow() public { assertEq(FixedPointMathLib.rpow(0, 0, 0), 0); assertEq(FixedPointMathLib.rpow(1, 0, 0), 0);