Skip to content

⚡️ Optimize solady.cbrt #1527

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

Merged
atarpara merged 5 commits into from
May 19, 2026
Merged

⚡️ Optimize solady.cbrt #1527

Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
7e8f57b
⚡️ Optimize cbrt
atarpara May 16, 2026
4e0e6eb
Typo
atarpara May 16, 2026
c8353b1
Fixed constant
atarpara May 17, 2026
9f046e0
Optimize more
atarpara May 17, 2026
3d4000b
New constant
atarpara May 19, 2026
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
14 changes: 7 additions & 7 deletions src/utils/clz/FixedPointMathLib.sol
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -806,14 +806,14 @@ library FixedPointMathLib {
function cbrt(uint256 x) internal pure returns (uint256 z) {
/// @solidity memory-safe-assembly
assembly {
// Initial guess z ≈ c · 2^q where b = ⌊log₂(x)⌋, q = ⌊b / 3⌋. The
// 8-bit fixed-point multipliers `c`: 144/128, 181/128, and 229/128
// Initial guess z ≈ c · 2^q where b = ⌊log₂(x) + 2⌋, q = ⌊b / 3⌋. The
// 8-bit fixed-point multipliers `c`: 89/128, 115/128, and 141/128
// are selected by `b mod 3` to balance each octave's worst-case
// final error. This gives >98 bits of precision after only 5
// Newton-Raphson iterations. The `or(..., 1)` keeps z ≥ 1 when the
// shifted estimate is 0.
let b := sub(255, clz(x))
z := or(shr(7, shl(div(b, 3), byte(add(mod(b, 3), 29), 0x90b5e5))), 1)
// final error. This gives >97 bits of precision after only 5

@duncancmt duncancmt May 17, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The comment is subtly wrong. The worst-case error analysis bound is 91 bits, not 97. This doesn't matter for implementation correctness, though

Suggested change
// final error. This gives >97 bits of precision after only 5
// final error. This gives >91 bits of precision after only 5

@atarpara atarpara May 18, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@duncancmt are u sure? As per my it is 97 bit accuracy

Best formula:  c_s = (89 + 26*s) / 128,   applied with s = mod(n+2, 3)

case n mod 3 = 0:
  n = 3k + 0
  z0           = 2^floor((3k+0+2)/3) = 2^(k+0) = 1 * 2^k
  x^(1/3)      in [1.0000, 1.2599) * 2^k
  ratio r      in (0.7937, 1.0000]
  error eps    in (-0.2063, +0.0000]
  worst |eps|  = 0.2063    (2.28 bits)

  with constant multiplier:
    s = (n+2) mod 3 = (3k+2) mod 3 = 2
    c[s]         = (89 + 26*2) / 128 = 141/128 = 1.101562
    new ratio r  in (0.8743, 1.1016]
    new eps      in (-0.1257, +0.1016]
    worst |eps|  = 0.1257   (2.99 bits)

case n mod 3 = 1:
  n = 3k + 1
  z0           = 2^floor((3k+1+2)/3) = 2^(k+1) = 2 * 2^k
  x^(1/3)      in [1.2599, 1.5874) * 2^k
  ratio r      in (1.2599, 1.5874]
  error eps    in (+0.2599, +0.5874]
  worst |eps|  = 0.5874    (0.77 bits)

  with constant multiplier:
    s = (n+2) mod 3 = (3k+3) mod 3 = 0
    c[s]         = (89 + 26*0) / 128 = 89/128 = 0.695312
    new ratio r  in (0.8760, 1.1037]
    new eps      in (-0.1240, +0.1037]
    worst |eps|  = 0.1240   (3.01 bits)

case n mod 3 = 2:
  n = 3k + 2
  z0           = 2^floor((3k+2+2)/3) = 2^(k+1) = 2 * 2^k
  x^(1/3)      in [1.5874, 2.0000) * 2^k
  ratio r      in (1.0000, 1.2599]
  error eps    in (+0.0000, +0.2599]
  worst |eps|  = 0.2599    (1.94 bits)

  with constant multiplier:
    s = (n+2) mod 3 = (3k+4) mod 3 = 1
    c[s]         = (89 + 26*1) / 128 = 115/128 = 0.898438
    new ratio r  in (0.8984, 1.1320]
    new eps      in (-0.1016, +0.1320]
    worst |eps|  = 0.1320   (2.92 bits)

Optimal integer multipliers (by s = mod(b, 3)):
  c[0] = 89/128 = 0.695312
  c[1] = 115/128 = 0.898438
  c[2] = 141/128 = 1.101562

Worst |eps_0| (with best c) = 0.1320  (2.92 bits)

Newton trace from worst eps_0 = +0.1320:
  step 0: eps =  +1.3196e-01, bits =    2.92
  step 1: eps =  +1.4786e-02, bits =    6.08
  step 2: eps =  +2.1439e-04, bits =   12.19
  step 3: eps =  +4.5948e-08, bits =   24.38
  step 4: eps =  +2.1112e-15, bits =   48.75
  step 5: eps =  +4.4573e-30, bits =   97.50  <-- target

@duncancmt duncancmt May 18, 2026 •
edited
Loading

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The sign of the relative error matters because there is an un-squared term in the relative error recurrence relation

def f(e):
    return e**2 * (3 + 2 * e) / (3 * (1 + e) ** 2)
def f5(e):
    return f(f(f(f(f(e)))))
def log_2(x, prec=Decimal('0.0001')):
    return (x.ln() / Decimal(2).ln()).quantize(prec)
es = [Decimal('-0.12396'), Decimal('+0.10374'), Decimal('-0.10156'), Decimal('+0.13196'), Decimal('-0\
.12569'), Decimal('+0.10156')]
max([log_2(f5(abs(e))) for e in es])
# Decimal('-97.5018')
max([log_2(f5(e)) for e in es])
# Decimal('-91.8558')

// Newton-Raphson iterations.

let b := sub(257, clz(x))
z := shr(7, shl(div(b, 3), add(89, mul(26, mod(b, 3)))))

// 5 Newton-Raphson iterations
z := div(add(add(div(x, mul(z, z)), z), z), 3)
Expand Down