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