Skip to content

Reduce gas of the fixed-point log10/pow10/pow/sqrt without losing precision #306

Closed

Description

@thedavidmeister

Follow-up to #297, which made log10, pow10, pow and sqrt accurate to about 1e-41 (pow10 correctly rounded to 41 digits) at a gas cost. Measured at #297's 98e8226, main → PR: sqrt(2) 11,820 → 51,431; log10(2) 6,036 → 31,678; pow10(0.5) 5,120 → 17,782.

Reduce that cost with precision unchanged. Every bc-pinned value and threshold test #297 lands must still pass, and exact results stay exact.

Levers to measure, each separately and combined:

  • A cheaper full-precision 512-bit mulDiv.
  • exp10 range reduction by binary digits against constants 10^(2^-i) compiled into the code, instead of 8 halvings and 8 squarings.
  • Fixed-degree polynomials with precomputed coefficients, instead of dividing by n on every term.
  • Native arithmetic wherever the bounds allow.

Report gas and accuracy against #297's merged state in one table.

Depends on #297.

Activity

  1. thedavidmeister commented on Oct 5, 2026

    @thedavidmeister
    ContributorAuthor

    Prototype state, 2026-10-05

    Branch 2026-10-05-issue-306-pow-gas, pushed, no PR, built on #297 as of 98e8226. It has to be rebuilt on #297's merged state before a PR.

    Gas with the same bench (GasBenchTest, about 2k call overhead included):

    Case #297 @ 98e8226 f0ab6d2 0a8065d 999c755 (head)
    sqrt(2) 50,164 19,309 3,914 4,014
    sqrt(4) 51,802 19,211 — —
    pow 2^1.5 53,064 23,405 10,073 —
    pow 7^1.234 54,603 23,386 — 23,014
    pow 0.5^0.3 55,681 20,610 — 19,981
    log10(2) 30,431 11,716 — 10,886
    log10(123456.789) 30,928 12,297 — 11,831
    pow10(0.5) 19,068 3,461 — —
    pow10(2) 2,071 2,071 — —

    (— = not re-measured at that commit.)

    Kept levers:

    • f0ab6d2: constant-divisor mulDiv, and exp10 range reduction by binary digits against in-code constants.
    • 0a8065d: sqrt by integer square root, correctly rounded at 41 digits. pow with b = n+½ goes through it. Against bc at scale 700 over 1,800 rows, no row is less accurate.
    • 999c755: scaleUp cascade replaces maximizeFull in log10. Outputs unchanged.

    Dropped: the generic cheaper mulDiv (−100 gas, not worth it). The precomputed-polynomial lever alone saved almost nothing.

    Uncommitted: log10Reduce, table-free log10 by binary reduction against 2^256-scale constants in code. Measured: log10(2) 10,886 → 8,618, log10(123456.789) → 8,831. Unrounded log10 error goes from about 7e-48 to about 1e-58. The broad bc comparison has not run yet. With this change nothing reads the log tables, so the tables data contract and its deploy could be retired.

    Rounding: log10 and pow10 both round to 41 significant digits. The 41 was chosen to match pow10's 1e50 scale with 10 guard digits, not derived. Next: derive the guard size from the measured raw error (at least about 100× below one output step) and keep every remaining digit.

  2. thedavidmeister commented on Oct 5, 2026

    @thedavidmeister
    ContributorAuthor

    Table-free log10, 53cb651

    log10Reduce divides the coefficient by 10^(2^-i) for i in 1..16, using reciprocals at the 2^256 scale written in code. This replaces the log table seed and the exp10Fixed estimate. log10, pow10, pow and sqrt no longer read tables.

    Bench (GasBenchTest, about 2k call overhead included). Error is relative, from bc.

    Case main 2b19ed1 gas main err 999c755 gas 999c755 err 53cb651 gas 53cb651 err
    sqrt(2) 10,544 1e-3.3 4,014 1e-40.7 4,014 1e-40.7
    sqrt(4) 16,616 exact 4,053 exact 4,053 exact
    sqrt(9) 10,577 exact 4,047 exact 4,047 exact
    pow 2^1.5 13,520 1e-3.3 10,173 1e-40.7 10,173 1e-40.7
    pow 7^1.234 20,196 1e-3.7 23,014 1e-40.9 21,464 1e-40.9
    pow 0.5^0.3 16,511 1e-3.9 19,981 1e-41.2 18,312 1e-41.2
    log10(2) 4,789 1e-4.0 10,886 1e-49.6 8,618 1e-54.3
    log10(123456.789) 16,124 1e-5.6 11,831 1e-50.4 8,831 1e-55.4
    pow10(0.5) 3,897 1e-4.1 3,461 1e-40.8 3,461 1e-40.8
    pow10(2) 3,532 exact 2,071 exact 2,071 exact

    log10 itself is unrounded, which is why its error is far below the 41-digit rounding that pow and pow10 apply. On main, log10(2) is cheap because main returns the 4-digit table value.

    Broad comparison against bc (scale 200) over 1,500 seeded rows: 500 log10 inputs with exponents from -200 to 200 and coefficients up to int224, 500 pow10 inputs, and 500 pow inputs with b from 1e-6 to 1e3.

    999c755 worst 53cb651 worst rows less accurate rows equal
    log10 absolute 1e-49.7 1e-54.2 0 111
    log10 relative 1e-50.4 1e-54.6 0
    pow10 relative 1e-40.3 1e-40.3 0 500
    pow relative 1e-40.3 1e-40.3 0 500

    pow and pow10 come out identical because their 41-digit rounding hides the log10 improvement. The full suite passes (542 tests).

    Not done here: the tables data contract and its deploy can now be retired, and the digit count still waits on #297's proven bound and correct-rounding tie check. Once that lands it gets merged into this branch, and the digit count is set from the bound.

  3. thedavidmeister commented on Oct 6, 2026

    @thedavidmeister
    ContributorAuthor

    Techniques from PRBMath worth adopting

    From reading PaulRBerg/prb-math (src/Common.sol, src/ud60x18/Math.sol). Its README averages for UD60x18: sqrt 710, log10 4,571, exp 2,244, pow 6,635.

    1. sqrt: fixed unrolled Newton from a bit-check estimate. Start at about 2^(msb/2) from a few bit checks, run exactly 7 r = (r + x/r) >> 1 with no loop or convergence test, then one final compare for rounding. Our integer sqrt (about 4,014 on Gas: cheaper fixed-point log10/pow10/pow/sqrt without losing precision (#306) #309) should adopt this.
    2. Constant-divisor division via a precomputed inverse (mulDiv18). For a fixed divisor d = 2^k·odd, take the remainder with mulmod, shift out 2^k, and multiply by odd's inverse mod 2^256, computed offline. Apply to the fixed 1e50, guard and power-of-ten divisors in the series. Check what Gas: cheaper fixed-point log10/pow10/pow/sqrt without losing precision (#306) #309's constant-divisor lever already does.
    3. Branch-free msb in assembly: factor := shl(k, gt(x, mask)); x := shr(factor, x); r := or(r, factor) over 8 steps. Adapt it to decimal digit counting (the log10 integer part, the scaleUp cascade).
    4. Group skipping in binary decomposition: test a byte of the exponent before its 8 bits. Apply to exp10/log10 binary range reduction.

    Not applicable: PRBMath's bit-per-squaring log2 loop (slower and less precise than our series), and its 64-bit exp2 precision (cheap only because it carries about 19 digits; we need about 41, hence 512-bit mulDiv).

  4. thedavidmeister commented on Oct 6, 2026

    @thedavidmeister
    ContributorAuthor

    PRBMath techniques applied, 4898920

    Two of the four techniques were kept, each as its own commit. The other two cost gas and were dropped. Every bench output is byte-identical to 53cb651. A differential fuzz against a copy of the 53cb651 implementation also matched on every run (100k runs each over sqrt, log10, pow10, the fractional pow path, mul, scaleUp, exp10Fixed and log10Reduce). The full suite passes.

    Lever Commit Result
    1. sqrt: fixed unrolled Newton 7459101 Kept. Math.sqrt is replaced by one shift step from 2^125, then 7 unrolled Newton steps and a floor correction. sqrt −534, pow 2^1.5 −534.
    2. Constant divisor by precomputed inverse 09fed5b 1e50 was already done (mulDivFixed, f0ab6d2), and the guard division is a plain 256-bit div. The one remaining generic division by a power of ten was mul's mulDiv(.., 10^n). It now shifts out 2^n and multiplies by exp(inverse of 5, n). Kept: pow 7^1.234 −412, 2^1.5 −309, 0.5^0.3 −103.
    3. Branch-free msb digit count — Dropped. scaleUp used an 8-step msb, (msb·1233)>>12 (checked against bc for msb 0–255) and one exp(10, k). The cost went up: sqrt +8, log10 +12 / +72, pow +8 to +12. The 7-compare cascade is cheaper.
    4. Group skipping — Dropped. Per-byte guards in exp10Fixed save 139 on pow10(0.5) but cost +44 on any input with bits in both bytes (pow 7^1.234, 0.5^0.3). In log10Reduce they cost +44 on every log10 and pow case.

    Bench: an external call from a test. "Net" subtracts a no-op function with the same signature called the same way (861 for one argument, 936 for two). Net is the function's own cost, comparable to PRBMath's README averages.

    Case 53cb651 gross 53cb651 net 7459101 net 09fed5b net (head) PRBMath UD60x18 avg
    sqrt(2) 3,645 2,784 2,250 2,250 sqrt 710
    sqrt(4) 3,645 2,784 2,250 2,250
    sqrt(9) 3,636 2,775 2,247 2,247
    pow 2^1.5 9,494 8,558 8,024 7,715 pow 6,635
    pow 7^1.234 20,731 19,795 19,795 19,383
    pow 0.5^0.3 17,585 16,649 16,649 16,546
    log10(2) 8,172 7,311 7,311 7,311 log10 4,571
    log10(123456.789) 8,384 7,523 7,523 7,523
    pow10(0.5) 3,027 2,166 2,166 2,166 exp 2,244
    pow10(2) 1,649 788 788 788

    The gross numbers come from a different harness than the earlier tables in this thread (that one had about 2k overhead), so compare net to net.

    Against PRBMath, at about 41 significant digits versus its 18: pow10 matches its exp. sqrt is about 3.2×, log10 about 1.6×, and pow 1.2× to 2.9×. The rest of sqrt is the 512-bit Newton step and square check, plus the unpack/pack around it. pow and log10 are dominated by the atanh series and the 16 reduction multiplies.

    Also on the branch: 4898920 trims CLAUDE.md under the 4096-byte context cap. That cap failed static from 53cb651 onward.

  5. thedavidmeister commented on Oct 6, 2026

    @thedavidmeister
    ContributorAuthor

    Techniques from other Solidity math libraries

    Sources read (fresh clones, 2026-10-06): Solady FixedPointMathLib (expWad, lnWad, powWad, sqrt, sqrtWad, cbrt, log2, log10, fullMulDiv), ABDKMath64x64 and ABDKMathQuad (sqrt, log_2, exp_2/pow_2, divuu, sqrtu), OpenZeppelin Math (mulDiv, mulShr, sqrt, log2, log10), Uniswap v3/v4 FullMath, v4 TickMath and BitMath, Balancer v2 LogExpMath (exp, _ln, _ln_36, pow), and Solmate FixedPointMathLib (an older Solady). The PRBMath items from the comment above are not repeated.

    Where the gas goes now

    Measured on 4898920, which has the same code as 09fed5b. These are internal calls timed with gasleft() in a scratch forge test, with no external call or pack/unpack, so they are lower than the bench above.

    Piece Gas
    log10(2) total 7,021
    of which log10Ratio (atanh series and 2 generic mulDiv) 4,329
    of which log10Reduce 839 to 1,116
    exp10Fixed(0.30103…) / (0.99999…) / (0.5) 3,326 / 5,289 / 520
    one mulDivFixed call ~180
    one generic mulDiv 548
    one inline 2^256 mulhi step (as in log10Reduce) ~62
    float mul / add / withTargetExponent / intFrac 1,194 / 2,342 / 381 / 285
    powFraction(7, 0.234) 13,774

    The cost is in mulDivFixed, called about 7 times in the atanh loop and up to 26 times in exp10Fixed. Each call costs about 3× the 512-bit multiply it performs.

    Ranked techniques

    1. Unrolled, call-free series and Horner steps. Balancer LogExpMath._ln_36 / exp, Solady expWad / lnWad.
    Balancer writes the atanh and Taylor series as a fixed number of straight-line steps, with no loop test and no checked arithmetic:

    num = (num * z_squared) / ONE_36;
    seriesSum += num / 3;
    num = (num * z_squared) / ONE_36;
    seriesSum += num / 5;   // ... to /15, "8 Taylor terms are sufficient"

    Solady keeps the whole polynomial in one assembly expression. For us, the same mulDivFixed arithmetic goes inline in a single assembly block, with no Solidity internal call and no Yul function, unrolled to the worst-case term count. Measured in a scratch prototype:

    • atanh series: 2,892 → 1,582 for z ≈ 1.7e-5, and 3,752 → 1,578 for z = 5e-4 (the unreduced path within 1e-3 of a power of ten). That is −1,310 to −2,170 on every log10 and pow. 8 terms cover z ≤ 5e-4. Once a term reaches 0 it stays 0, so the extra terms add nothing.
    • exp10Fixed: 3,309 → 2,572 (0.30103), 5,269 → 3,940 (0.99999), 2,339 → 1,890 (1.23e-5). That is −450 to −1,330 on pow10 and pow.

    Wrapping the same arithmetic in a Yul function inside the block saved almost nothing (3,309 → 3,277), so the full expansion is what saves the gas.
    Precision: bit-identical. A differential fuzz (5,096 runs each) matched the current exp10Fixed and the current loop on every input.

    2. Binary 2^256 basis with a high-word multiply. Solady expWad / lnWad (Remco Bloemen), Uniswap TickMath.getSqrtPriceAtTick, ABDK exp_2 / pow_2.
    All three accumulate in a binary scale, so the fixed-point rescale is a shift or a high word rather than a division by 10^k:

    // Solady expWad: "Convert to (-42, 136) * 2**96 for more intermediate precision and a binary basis"
    x = (x << 78) / 5 ** 18;
    // TickMath / ABDK
    if (absTick & 0x2 != 0) price = (price * 0xfff97272373d413259a46990580e213a) >> 128;

    Our log10Reduce already does this, at about 62 gas a step against about 130 for an expanded mulDivFixed. Applying it to the exp10Fixed accumulator and polynomial, and to the atanh series, would need four things:

    • the result held at 2^252, with 10^(2^-i) constants at 2^254 and shl(2, mulhi(..));
    • the Taylor coefficients, and the atanh 1/(2k+1) as a Horner polynomial in z², at a binary scale, which also removes the per-term div;
    • one conversion in, from the 1e50 remainder;
    • one conversion out, (r * 1e50) >> 252, a mul512 and a shift.

    A gas-only stand-in, 9 Horner steps with expanded mulhi, cost 563 against 1,582 for the expanded decimal series.
    Estimate on top of #1, conversions included: atanh about −800 to −900, and exp10Fixed about −600 to −1,200. Not prototyped end to end.
    Precision: changes raw values. Each step truncates at about 2^-252 (1e-76) instead of 1e-50, so the raw error gets smaller. Outputs are no longer bit-identical, so the bc comparison has to be rerun. The #308 error bound has to be re-derived with the new per-step error. A 41-digit rounded result can differ from today only where today's result was within the old error of a tie.

    3. Keep pow's b·log10(a) in fixed point. Solady powWad, Balancer pow.

    return expWad((lnWad(x) * y) / int256(WAD));  // Solady
    logx_times_y = ((ln_36_x / ONE_18) * y_int256 + ((ln_36_x % ONE_18) * y_int256) / ONE_18);  // Balancer

    Neither library turns the log back into a general number between ln and exp. Today powFraction returns log10's float, runs a float mul (1,194), and then pow10 runs intFrac (285) and two withTargetExponent (2 × 381). log10's small-exponent branch already builds exponent * 1e66 + seed * 1e16 + correction at 1e-66. Multiplying that fixed-point value by b's coefficient and handing the integer and fraction straight to exp10Fixed would skip all of those steps.
    Estimate: pow −1,500 to −2,000. Not prototyped.
    Precision: changes raw values, probably for the better. The float mul truncates its coefficient, and this path would not. The fixed-point product needs at least about 51 fraction digits. It also needs a fallback to today's path when |b·log10(a)| is too large for the fixed-point word, or when log10 takes its large-exponent add branch. The #308 bound has to cover the new path.

    4. One division by a + b in log10Ratio, not two. Balancer _ln_36.
    Balancer computes z = (x−1)/(x+1) once and multiplies it through the series (num = z; num = num * z_squared …). We divide by sum twice, once for z and once more at the end to keep relative precision when difference is tiny. Both divisions are generic mulDiv at 548 each. Instead, scale difference into [1e74, 1e76) first, using the cascade that is already there. Then take one mulDiv(scaledDifference, 1e50, sum) as the relatively precise z, derive z² from it with the existing mulDivPow10, and multiply z into the series result.
    Estimate: −350 to −450 on log10 and pow (one mulDiv saved, one mulDivPow10 added). Not prototyped.
    Precision: changes raw values. The last-digit truncation moves, and the bound needs a line for it.

    5. sqrt rounding check from low words only. Adapted from Solady sqrtWad, which checks rounding with a single-word residue, z := sub(z, gt(999999999999999999, sub(mulmod(z, z, x), 1))).
    For us, |N − root²| ≤ 2·root + 1, which is far below 2^255. So sub(mul(coefficient, scale), mul(root, root)), read as signed, is the exact residual, and the two mul512 calls and the 512-bit compare can go:

    let d := sub(mul(coefficient, scale), mul(root, root))
    if slt(d, 0) { d := add(d, sub(shl(1, root), 1)) root := sub(root, 1) }
    root := add(root, gt(d, root))

    Measured: sqrt 1,860 → 1,545 (−315), and pow with b = n+½ by the same amount.
    Precision: bit-identical. A fuzz of 5,096 runs over coefficients up to int224 and exponents ±1e6 matched on every run.

    6. Cheap estimate plus an exact check only when ambiguous. Uniswap v4 TickMath.getTickAtSqrtPrice.

    int24 tickLow = int24((log_sqrt10001 - 3402992956809132418596140100660247210) >> 128);
    int24 tickHi = int24((log_sqrt10001 + 291339464771989622907027621153398088495) >> 128);
    tick = tickLow == tickHi ? tickLow : getSqrtPriceAtTick(tickHi) <= sqrtPriceX96 ? tickHi : tickLow;

    This is Ziv's strategy with a proven error interval. For #308's correct rounding, the tie check runs only when the 41-digit rounding point falls inside the proven error interval. The same pattern would let the common path carry fewer guard digits (fewer series terms and reduction steps), and recompute at full precision only near a tie.
    Saving: depends on #308's design. Nothing is measured.
    Precision: the contract holds only if the interval is proven. That is the #308 dependency.

    Read and rejected for our precision

    • Rational (Padé) approximations (Solady expWad (6,7)-term, lnWad (8,8)-term). They are cheap at 1e-18. At 1e-51 with r ≤ 2^-16·ln10, a [5/5] Padé needs about 10 multiplies and a 512-bit division, against 9 Horner multiplies now. Even/odd splitting with less range reduction does not beat the current cost either.
    • Bit-per-squaring log (ABDK log_2, TickMath getTickAtSqrtPrice, Balancer-style ladders). One bit per squaring would mean about 170 512-bit squarings at our precision.
    • 128-bit constants with one mul … >> 128 (ABDK exp_2, TickMath). 2^-128 is about 3e-39, below 41 digits plus guard digits. That is why we need 512-bit products.
    • Fewer reduction steps and more series terms (Balancer stops _ln at 2^-4 and exp at 2^-2). With the costs measured above, r steps plus n terms is flat across r = 8 to 16 for our target. No saving.
    • Branchless first step (v4 price := xor(ONE, mul(xor(ONE, C), and(absTick, 1)))). Our exp10Fixed first step already assigns a constant, and a branchless mulhi costs more than the JUMPI it saves.
    • Solady fullMulDiv overflow fast path (eq(div(z, x), y)) and its fused final Newton step. Tens of gas on the two generic mulDiv calls. 2024 07 02 decimal lossy #4 removes one of them anyway.
    • Solady / OZ log10 digit cascades and the OZ log2 nibble byte() table. Already covered by the dropped branch-free digit count above.
    • via_ir, which would inline the mulDivFixed calls. It changes all bytecode and the deploy pins in rain.math.float.deploy. 2024 06 06 log #1 gets the same gain locally.

    Published gas numbers for comparable functions

    Library Function Published gas Source
    PRBMath v2.0.1 UD60x18 sqrt / log10 / ln / log2 / exp / exp2 / pow avg 710 / 4,571 / 3,814 / 3,426 / 2,244 / 2,156 / 6,635 (max 846 / 8,695 / 6,902 / 6,825 / 2,742 / 2,652 / 10,637) old README (8ff5a3f). The current README is qualitative only and says there are no automated gas reports.
    PRBMath v2.0.1 SD59x18 sqrt / log10 / ln / exp / pow avg 716 / 4,337 / 4,724 / 2,263 / 8,518 same
    ABDKMath64x64 v3.0 sqrt / log2 / ln / exp / exp2 / pow (integer exponent) avg 699 / 7,024 / 7,126 / 2,687 / 2,746 / 1,792 PRBMath README (8ff5a3f), from paulrberg/abdk-gas-estimations
    Uniswap v4 getSqrtPriceAtTick / getTickAtSqrtPrice 72,349 / 195,022 per 100 calls in a loop (about 723 / 1,950 each) snapshots/TickMathTestTest.json
    Solmate sqrt fuzz μ 1,035 (test harness) .gas-snapshot
    Solady, OpenZeppelin, Balancer, ABDKMathQuad — none published

    Precision context: PRBMath and Balancer carry about 18 digits (Balancer's _ln_36 carries 36), ABDK 64.64 about 19, ABDKMathQuad 113 bits (about 34 digits), and we carry 41.

    Scratch prototypes and the profiling test are local only. No code changes, and no PR.

  6. thedavidmeister commented on Oct 6, 2026

    @thedavidmeister
    ContributorAuthor

    Constraint on the binary 2^256-scale lever (user, 2026-10-06): "we want decimal floats, so if we use binary scale there has to be some conversion".

    So the binary scale is internal only:

    • Decimal → binary on the way in: one constant mulDiv.
    • Series in binary: the per-step saving.
    • Binary → decimal on the way out: one constant mulDiv, then the existing round to 41 significant decimal digits.

    Conditions:

    • The two conversion errors are added to the proven error bound, with the analysis redone and written next to the code.
    • Exact results never go through the binary path. sqrt stays integer decimal. Powers of ten, integer powers and log10(10^k) keep their exact early returns before any conversion.
    • The documented contract (rounded to nearest at decimal digits, monotone within 2E with the one-ulp flip, exact anchors) is unchanged.
    • Take the lever only if gas net of both conversions is lower, measured per function.
  7. added 3 commits that reference this issue on Oct 6, 2026
  8. thedavidmeister commented on Oct 11, 2026

    @thedavidmeister
    ContributorAuthor

    Closing per the user's ruling. #297 is closed. Gas without precision loss landed in #352 (inlined log10Ratio/exp10Fixed series: log10 −1759, pow10 −494) and #432 (sqrt seed, 4 Newton steps instead of 7: −28 to −47 on every sqrt and half-integer pow). Outputs are identical to main.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions