From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on inbox.vuxu.org X-Spam-Level: X-Spam-Status: No, score=-3.1 required=5.0 tests=DKIM_INVALID,DKIM_SIGNED, HEADER_FROM_DIFFERENT_DOMAINS,MAILING_LIST_MULTI,RCVD_IN_DNSWL_MED autolearn=ham autolearn_force=no version=3.4.4 Received: from second.openwall.net (second.openwall.net [193.110.157.125]) by inbox.vuxu.org (Postfix) with SMTP id 2FB702E05D for ; Thu, 1 May 2025 22:07:14 +0200 (CEST) Received: (qmail 27848 invoked by uid 550); 1 May 2025 20:07:10 -0000 Mailing-List: contact musl-help@lists.openwall.com; run by ezmlm Precedence: bulk List-Post: List-Help: List-Unsubscribe: List-Subscribe: List-ID: Reply-To: musl@lists.openwall.com x-ms-reactions: disallow Received: (qmail 27823 invoked from network); 1 May 2025 20:07:10 -0000 DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=pkh.me; s=selector1; t=1746130020; h=from:from:reply-to:subject:subject:date:date:message-id:message-id: to:to:cc:cc:mime-version:mime-version:content-type:content-type: content-transfer-encoding:content-transfer-encoding; bh=D8UCwlXWhZz4n6nUK4ZtSiJUobm8VVhI4J65WucL3xs=; b=Of6Odk39syn36aRgz+2nXKzMdQK4iLeN12or5T4jb8afaDQSI89xq3PAartdLhDG1QbcsN wG9TVV2kR7GQFFWJu+bmsrZ1uyPFOsJYtig2Vl2YK7pRbnRFYtOTC8JYPQMH2e/lhNvcNK Rlb7fqomna8m1S5byK8JfIX0SPORY5Q= From: =?UTF-8?q?Cl=C3=A9ment=20B=C5=93sch?= To: musl@lists.openwall.com Cc: =?UTF-8?q?Cl=C3=A9ment=20B=C5=93sch?= Date: Thu, 1 May 2025 22:04:39 +0200 Message-ID: <20250501200648.2878841-1-u@pkh.me> X-Mailer: git-send-email 2.41.0 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Subject: [musl] [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ₙ) = = 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 --- 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