Skip to content

⚡️ Cheaper expWad and lnWad via symmetric rationals #1546

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ddallaire wants to merge 1 commit into
base: main
Choose a base branch
from
Open

⚡️ Cheaper expWad and lnWad via symmetric rationals #1546

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 37 additions & 11 deletions docs/utils/fixedpointmathlib.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand All @@ -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`.
Expand Down Expand Up @@ -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)

Expand All @@ -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)

Expand Down
157 changes: 157 additions & 0 deletions src/utils/FixedPointMathLib.sol
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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.
Expand Down
Loading