[<prev] [next>] [day] [month] [year] [list]
Message-ID: <20250501200648.2878841-1-u@pkh.me>
Date: Thu, 1 May 2025 22:04:39 +0200
From: Clément Bœsch <u@....me>
To: musl@...ts.openwall.com
Cc: Clément Bœsch <u@....me>
Subject: [PATCH] math: fix comments about cbrt iterations being Halley and not Newton
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¹ 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² 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.
A similar patch has been sent and accepted³ in the FreeBSD repository
where this code originates from.
¹ https://en.wikipedia.org/wiki/Halley's_method
² https://www.wolframalpha.com/input?i=T-2%28T%5E3-x%293T%5E2%2F%282%283T%5E2%29%5E2-%28T%5E3-x%296T%29
³ https://github.com/freebsd/freebsd-src/pull/1684
Signed-off-by: Clément Bœsch <u@....me>
---
Note: I'm not subscribed to the mailing-list, so CCs are appreciated
---
src/math/cbrt.c | 4 ++--
src/math/cbrtf.c | 4 ++--
src/math/cbrtl.c | 2 +-
3 files changed, 5 insertions(+), 5 deletions(-)
diff --git a/src/math/cbrt.c b/src/math/cbrt.c
index 7599d3e3..42ac9520 100644
--- a/src/math/cbrt.c
+++ b/src/math/cbrt.c
@@ -84,7 +84,7 @@ double cbrt(double x)
* the result is larger in magnitude than cbrt(x) but not much more than
* 2 23-bit ulps larger). With rounding towards zero, the error bound
* would be ~5/6 instead of ~4/6. With a maximum error of 2 23-bit ulps
- * in the rounded t, the infinite-precision error in the Newton
+ * in the rounded t, the infinite-precision error in the Halley
* approximation barely affects third digit in the final error
* 0.667; the error in the rounded t can be up to about 3 23-bit ulps
* before the final error is larger than 0.667 ulps.
@@ -93,7 +93,7 @@ double cbrt(double x)
u.i = (u.i + 0x80000000) & 0xffffffffc0000000ULL;
t = u.f;
- /* one step Newton iteration to 53 bits with error < 0.667 ulps */
+ /* one step Halley iteration to 53 bits with error < 0.667 ulps */
s = t*t; /* t*t is exact */
r = x/s; /* error <= 0.5 ulps; |r| < |t| */
w = t+t; /* t+t is exact */
diff --git a/src/math/cbrtf.c b/src/math/cbrtf.c
index 89c2c865..639c0d3e 100644
--- a/src/math/cbrtf.c
+++ b/src/math/cbrtf.c
@@ -46,7 +46,7 @@ float cbrtf(float x)
u.i |= hx;
/*
- * First step Newton iteration (solving t*t-x/t == 0) to 16 bits. In
+ * First step Halley iteration (solving t*t-x/t == 0) to 16 bits. In
* double precision so that its terms can be arranged for efficiency
* without causing overflow or underflow.
*/
@@ -55,7 +55,7 @@ float cbrtf(float x)
T = T*((double_t)x+x+r)/(x+r+r);
/*
- * Second step Newton iteration to 47 bits. In double precision for
+ * Second step Halley iteration to 47 bits. In double precision for
* efficiency and accuracy.
*/
r = T*T*T;
diff --git a/src/math/cbrtl.c b/src/math/cbrtl.c
index ceff9136..02d23eb5 100644
--- a/src/math/cbrtl.c
+++ b/src/math/cbrtl.c
@@ -109,7 +109,7 @@ long double cbrtl(long double x)
#endif
/*
- * Final step Newton iteration to 64 or 113 bits with
+ * Final step Halley iteration to 64 or 113 bits with
* error < 0.667 ulps
*/
s = t*t; /* t*t is exact */
--
2.49.0
Powered by blists - more mailing lists
Confused about mailing lists and their use?
Read about mailing lists on Wikipedia
and check out these
guidelines on proper formatting of your messages.
