⚡️ Cheaper expWad and lnWad via symmetric rationals

[ddallaire]

Description

Cheaper expWad and lnWad using symmetric rational approximations, and a round-to-nearest 2⁹⁶-precision final step for lambertW0Wad.

expWad: a diagonal rational approximant of exp satisfies q(v) = p(-v), so splitting the numerator into even and odd parts p = E(v²) + v·O(v²) makes the denominator E − v·O, i.e. free. The core drops from 8 muls / 11 constants to 5 / 5 (−27% gas).

lnWad: ln(x) = ln2/2 + 2·atanh(s) with s = (x − √2)/(x + √2). atanh is odd, so the rational only has odd terms, and the √2 centering costs nothing (the tail's k·ln2 just becomes (2k+1)·ln2/2) while shrinking the fit domain enough for a (3,3) rational in s². 12 muls / 14 constants become 7 / 7 (−13% gas). Floor rounding convention unchanged, none of the existing ln test vectors move, and lnWad(1e18) is now exactly 0.

lambertW0Wad: the final Iacono–Boyd step inherited lnWad's last-wei rounding noise, which is what the hand-tuned nudges were calibrated against. It now inlines the ln core at 2⁹⁶ precision and rounds once, putting W within one wei of correctly rounded in that branch (98% exactly rounded, from 46%), and slightly cheaper overall since the cheaper exp/ln inside the iterations pay for it.

Measured accuracy is unchanged against a 140-digit reference: exp at the same 9.68e-18 pipeline floor with the same worst cases, ln within 1 wei. Every constant is reproducible from a script, unlike the previous ones (see the note in the test file). Details, derivations, benchmarks, and proofs here.

Also updated the exact test vectors that moved; several known values are now the correctly rounded ones (exp(-3), the omega constant, W(2**71 - 1), W(2**80 - 1)), and the rest reshuffle by a wei in both directions within the unchanged error envelope, net toward truth (details in the repo above).

Checklist

Ensure you completed all of the steps below before submitting your pull request:

• Ran forge fmt?
• Ran forge test?

[atarpara]

@ddallaire Thanks for this PR! @duncancmt and I worked on the same thing (lnWad) using clz a few months ago: #1537. I’ll review your approach as well.

I remember that at the time, we had an issue with the 18-decimal precision — we couldn’t get the true floor version of lnWad. @duncancmt then came up with a solution using ray precision.

[ddallaire]

@atarpara I'll be waiting for comments 👍 Note that one CI test run is failing as the test outgrew the configured gas limit of 100M, please tell me if you want to bump the gas limit or no, I can always split the test contract on its own for labertW0Wad.

[duncancmt]

My PR (and un-PR'd complementary work on exp in the same repo as referenced by the ln PR) are somewhat orthogonal to @ddallaire 's work here. @ddallaire 's work targets lambertW0 while mine doesn't bother. My work extends ln and exp to work in Ray precision (rather than Wad) because these values are almost always rates or rate-like; by Maker convention, these should be Ray, not Wad. As a consequence of the larger fixnum format, my implementation is less gas efficient.

Independent of the fixnum format, my implementation carefully quantifies the error in the approximation and compensates for/bounds it. Fundamentally, because ln and exp are both[0] transcendental, due to the table-maker's dilemma, without significant knowledge of or constraint on the domain, it's not possible to produce a exactly-computed correctly-rounded result in bounded computation (gas). Therefore, my implementation does the best that is possible in bounded computation: at most 1 wei of error and only error in one direction.

@ddallaire , does your implementation have this rounding behavior? If not, that would be a nice improvement. I'm happy to work with you on that and try to formally verify that property using the same stack that I did with my implementation.

---

[0] Although sqrt and cbrt are transcendental, because their inverses are not, it is possible to produce an exactly-computed correctly-rounded result in bounded computation by running the inverse (squaring or cubing) and comparing the result.

[ddallaire]

@duncancmt Thanks for the context on previous work. The goal of this implementation was to optimize the gas while being equally or more accurate than the previous version. If the current efforts are towards accuracy, you can disregard this PR. Error stats for these, feel free to validate:

• lnWad: error in [−1.04, +0.04] wei. ~99% of inputs land below the true value (the final wad conversion is a floor), ~1% dip just past −1 wei, and a small tail lands up to +0.04 wei above.
• expWad x > 0: relative error in [−2.4e−18, +1.7e−18], both directions, so the wei error scales with the result, as it must for a relative fit.
• expWad x ≤ 0: absolute error in [−1.9, +1.1] wei.

The 6,6 version from my repo which also cuts down on gas reduces the error by a lot.

[duncancmt]

Well, whether or not the stated error is adequate is a bit putting the cart before the horse. I think the major question is whether to keep the existing expWad/lnWad functions or whether to go to Ray-basis (at the expense of gas) while keeping Wad-basis wrappers around for backwards-compatibility. Answering that question rests entirely with @atarpara .

Independent of the answer to that question, I think it's very likely that your implementation could be coaxed into having floor-or-1-wei-less behavior with a bit of tweaking to the constants and a subtraction. My professional opinion is that having one-sided rounding error makes the design of systems with safe/conservative/optionality-adverse rounding behavior considerably easier, but this too could be a point for discussion.

Certainly either implementation would be an improvement on the current one! 👍

[ddallaire]

I'm not going to chime in on the wad vs ray discussion as I have no horse in this.

If the decision ends up going towards supporting both, I'd be happy to tweak both exp and ln function to make the error one sided as my PR doesn't supported that as is, my (6,6) exp implementation can be one-sided by removing a small bias, at the cost of 8 gas over the (5,5), unsure about the ln function as it's centered on lnWad(1e18) == 0

Ref:
https://github.com/ddallaire/wad-exponentials/blob/d1056df258aaea6ac4b0ed7dfdf740b9b311d691/src/ExpWad.sol#L93

[atarpara]

@ddallaire I agree with the design direction suggested by @duncancmt , but I don't think Solady needs to move to Ray precision here.

I'd prefer to keep the existing Wad/18-decimal basis and optimize around that. For Solady, the important property isn't higher internal precision, but rather that the approximation error is bounded and, more importantly, that the rounding direction is guaranteed.

So I'd be happy with a solution that stays at Wad precision as long as we can guarantee that the result is always rounded in one direction (e.g. true floor for lnWad), with the error bounded to at most 1 wei.

I think that property is more important for Solady than squeezing out additional precision, especially since it keeps the behavior predictable while keeping the implementation simple and efficient.

And yes, I'd definitely like to formally verify that rounding-direction guarantee if possible. Maybe @duncancmt can help us with that.

[duncancmt]

You can check my full implementation (again, at Ray-scale, but the general techniques may be applicable) in the branch in this PR 0xProject/0x-settler#599 . This is part of a larger PR stack, so I recommend pulling just the branch involved and examining the code/proof.

I'm surprised that the 6-coeff version of exp is required to get the desired one-sided rounding error. Perhaps the technique from my branch that involves doing a raw search over the low bits of the coeffs can be applied here? You may also find it possible to obtain smaller code size by adopting the mixed fixed-point Horner technique and sharing the constant coefficient between both polynomials (applicable to both ln and exp).

[Vectorized]

Let’s support both.

Just postfix the faster functions with “Fast”, while keeping the original functions in case some applications require the original precision.

FixedPointMathLib is like a museum on the different numerical techniques.

lambertW0Wad is really niche, so we can leave it using the original lnWad and expWad.

In the comments, also give a link to the 0x-settler implementations in case viewers want to use the ray versions.