Seeds the sqrt Babylonian/Newton iteration from a float estimate - #245
Merged
Conversation
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. |
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. |
Contributor
Author
|
See #247 |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
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/1is correctly rounded (a hard IEEE 754 guarantee, unlikethe 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 invariantrequires, 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 sameway. 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.