Fix precision loss in Newton-Raphson inverse - #552
Open
tompng wants to merge 1 commit into
Open
Conversation
The residual precision `SIZET2NUM(SIZET2NUM(n / 2))` was a double conversion. It accidentally requested about n digits, and the intended n / 2 is not enough: the digit schedule `n = (prec >> i) + 2` leaves only 1-2 margin digits, so rounding errors compound over iterations and the inverse loses many digits. The correction loop in divmod_by_inv_mul then runs once per lost unit of the quotient and a single divmod can take minutes. With ROUND_DOWN or ROUND_FLOOR the same divergence happened even with the accidental n-digit residual, because the iteration converges from below and truncation adds error in the same direction at every step. - Use 4 margin digits for the schedule and n / 2 + 4 digits for the residual - Calculate with ROUND_HALF_UP in VpDivdNewton and restore the rounding mode together with the precision limit - Initial approximation from the leading 18 digits with 64-bit division (about 8 correct digits instead of 2), and skip iterations below that. The previous one could produce frac[0] == BASE for x like 1.001 - Pass VpDivdNewtonInner arguments as a struct Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
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 residual precision
SIZET2NUM(SIZET2NUM(n / 2))was a double conversion. It accidentally requested about n digits, and the intended n / 2 is not enough: the digit schedulen = (prec >> i) + 2leaves only 1-2 margin digits, so rounding errors compound over iterations and the inverse loses many digits. The correction loop in divmod_by_inv_mul then runs once per lost unit of the quotient and a single divmod can take minutes. With ROUND_DOWN or ROUND_FLOOR the same divergence happened even with the accidental n-digit residual, because the iteration converges from below and truncation adds error in the same direction at every step.Additional information:
Initial approximation was 2-digits precision because BASE_FIG was 4(uint64 unavailable) or 9(uint64 available) before. 2-digits is the half of BASE_FIG=4.
BASE_FIG is now always 9, so we can increase the initial approximation and the margin.
(Larger margin requires larger initial approximation precision)