HomeFreeBSD

msun: fix cbrt iterations from Newton to Halley method

Description

msun: fix cbrt iterations from Newton to Halley method

Since we're inverting a cube, we have:

f(Tₙ)=Tₙ³-x    (1)

Its first and second derivatives are:

f'(Tₙ)=3Tₙ²    (2)
f"(Tₙ)=6Tₙ     (3)

Halley iteration[1] uses:

Tₙ₊₁=Tₙ-2f(Tₙ)f'(Tₙ)/(2f'(Tₙ)²-f(Tₙ)f"(Tₙ))    (4)

Replacing the terms of (4) using (1), (2) and (3):

Tₙ₊₁ = Tₙ-2f(Tₙ)f'(Tₙ)/(2f'(Tₙ)²-f(Tₙ)f"(Tₙ))

= Tₙ-2(Tₙ³-x)3Tₙ²/(2(3Tₙ²)²-(Tₙ³-x)6Tₙ)
= <snip, see WolframAlpha[2] alternate forms>
= Tₙ(2x+Tₙ³)/(x+2Tₙ³)

This formula corresponds to the exact expression used in the code.

Newton formula is Tₙ-f(Tₙ)/f'(Tₙ) which would have simplified to
(2Tₙ³+x)/(3Tₙ²) instead.

[1] https://en.wikipedia.org/wiki/Halley's_method
[2] https://www.wolframalpha.com/input?i=T-2%28T%5E3-x%293T%5E2%2F%282%283T%5E2%29%5E2-%28T%5E3-x%296T%29

Note: UTF8 in commit message due to the heavy math being hard to
recreate w/o it. -- imp

Signed-off-by: Clément Bœsch <u@pkh.me>
Reviewed by: imp
Pull Request: https://github.com/freebsd/freebsd-src/pull/1684

Details

Provenance
Clément Bœsch <u@pkh.me>Authored on May 1 2025, 5:19 PM
impCommitted on May 5 2025, 4:52 AM
Parents
rG5b9660caff69: Fix incorrect version introduced in manual pages
Branches
Unknown
Tags
Unknown