Skip to content

Seeds the sqrt Babylonian/Newton iteration from a float estimate - #245

Merged
ericmj merged 1 commit into
ericmj:mainfrom
tomciopp:sqrt-float-seed
Aug 26, 2026
Merged

Seeds the sqrt Babylonian/Newton iteration from a float estimate#245
ericmj merged 1 commit into
ericmj:mainfrom
tomciopp:sqrt-float-seed

Conversation

@tomciopp

Copy link
Copy Markdown
Contributor

The power of ten seed overshoots the true root by up to an order of
magnitude, so the loop takes eight or more bignum divisions to converge.
A double holds the scaled operand with at most 2^-52 relative error and
:math.sqrt/1 is correctly rounded (a hard IEEE 754 guarantee, unlike
the transcendental functions), so its float square root sits within a
few ULPs of the true root. Widened by 1e-12, (orders of magnitude beyond
that error bound) and stepped past truncation, the seed provably starts
at or above floor(√operand), which the loop's from above invariant
requires, while overshooting so little that convergence takes two or
three divisions. Since each Newton step roughly doubles the number of
correct digits, starting from the double's ~16 correct digits instead of
zero is where the entire saving comes from.

This float seeded integer Newton pattern is the established approach in
arbitrary precision libraries. The GNU Multiple Precision Arithmetic
Library (GMP) seeds its integer square root (mpn_sqrtrem) the same
way. That's cited as algorithmic precedent only: the BEAM has its own
bignum implementation and no GMP dependency. Correctness rests on
IEEE 754's correctly-rounded sqrt plus the classical from above
convergence invariant of the integer Newton iteration.

Operands at or above 10^300 fall back to the power of ten seed, since
converting them to a double would overflow. The scaled operand has
roughly 2 × precision digits, so the fallback only triggers for
precisions in the hundreds.

Results are identical: the loop converges to exactly floor(√operand)
from any seed at or above it. New tests pin inexact roots at the default
precision 34, including the odd exponent scaling path.

bench.exs sqrt job (Apple M5, OTP 29, MIX_ENV=prod):
1680.87 μs → 745.48 μs per pass (2.25x faster). Memory on the job rises
458.69 KB → 535.56 KB (~48 B per call for the float boxes the seed
computation allocates). All other jobs unchanged beyond run to run
noise; their memory measurements are identical.

@ericmj
ericmj merged commit 5d363f9 into ericmj:main Aug 26, 2026
2 checks passed
@ericmj

ericmj commented Aug 26, 2026

Copy link
Copy Markdown
Owner

@tomciopp Thanks for all your fixes and performance improvements! Are you planning any more PRs? Otherwise I will make a new release in the coming days.

@tomciopp

Copy link
Copy Markdown
Contributor Author

@ericmj I don't have anything queued, I will do one more pass and send an update later today whether I find something or not.

@tomciopp

Copy link
Copy Markdown
Contributor Author

See #247

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

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants