From nobody Mon May 5 05:06:02 2025 X-Original-To: dev-commits-src-all@mlmmj.nyi.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2610:1c1:1:606c::19:1]) by mlmmj.nyi.freebsd.org (Postfix) with ESMTP id 4ZrTy31N0Gz5tsZR; Mon, 05 May 2025 05:06:03 +0000 (UTC) (envelope-from git@FreeBSD.org) Received: from mxrelay.nyi.freebsd.org (mxrelay.nyi.freebsd.org [IPv6:2610:1c1:1:606c::19:3]) (using TLSv1.3 with cipher TLS_AES_256_GCM_SHA384 (256/256 bits) key-exchange X25519 server-signature RSA-PSS (4096 bits) server-digest SHA256 client-signature RSA-PSS (4096 bits) client-digest SHA256) (Client CN "mxrelay.nyi.freebsd.org", Issuer "R11" (verified OK)) by mx1.freebsd.org (Postfix) with ESMTPS id 4ZrTy25ryWz3Frx; Mon, 05 May 2025 05:06:02 +0000 (UTC) (envelope-from git@FreeBSD.org) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=freebsd.org; s=dkim; t=1746421562; h=from:from:reply-to:subject:subject:date:date:message-id:message-id: to:to:cc:mime-version:mime-version:content-type:content-type: content-transfer-encoding:content-transfer-encoding; bh=D0VURfU0QQESbjTJ0OhUUh9pjiXzobIvu9vmg0S5J4I=; b=eHDJaKxiZmYz9YrMfu95EGXhcvLUlMLhFOoXZZ1YJCBjGLLI2lTpKRTNBerv+J2oBxRgk0 WLUGtloD/N1blKb2nkdxFv4i8vczUIP7zoTm4sITLMJG7KXpfa+sErDgDjtcl8U8Ezh8u1 CxkVjDK+u1AaSQ9cvAjo8xsyJ1d5VUb5dSohWqVXnlPYwSZblzmVpwNNbgp9enoA2Lbiwd XPvR97UBNUrPdLAOCN6vZUhZ/KeHqP4l96SkUdjElGAgQvcdtgsMInV+U7tQrDauyEswxk QGu/8RR2NfYFJH+oS4gRNxMTJTxbeDrrVQv37qxbkAFYyu0LxrI+eHdgKYbNoA== ARC-Message-Signature: i=1; a=rsa-sha256; c=relaxed/relaxed; d=freebsd.org; s=dkim; t=1746421562; h=from:from:reply-to:subject:subject:date:date:message-id:message-id: to:to:cc:mime-version:mime-version:content-type:content-type: content-transfer-encoding:content-transfer-encoding; bh=D0VURfU0QQESbjTJ0OhUUh9pjiXzobIvu9vmg0S5J4I=; b=Tu17n+PIACxtwyYnLsE95oynb61t7rBEEYh4vLvHrWd9ocnJMKYDY2T45d6MWTqfoWsCvP hcabKUJzRq9B1OtEQYx322B4dIx5qvkc7E1lUu+JrzAbngdEXK4cevsbWU1fjupTPQDsXR ZGyoPa0P5RvEHmNgcailPgGSfk6EQ//rdFnfO4MTd0D/BIaXZ6cq6dI3iZbkS8oJxKjKrz rlx7y66LVSlscj7COk46uMpTA53RDUdY9VL7EmV70fL2klCibE/1iaACGkt1XIbns3g1Qx ZLQgDfyJohiGvf6EJSJf/u+4eKbbghIznqTVhu+wvfnYLpiUE5rAE0siWP4EHA== ARC-Seal: i=1; s=dkim; d=freebsd.org; t=1746421562; a=rsa-sha256; cv=none; b=XxUOs5qaua6gIFo/uvaLxYqoc3DI2+TW43t2vNU4krXGVDH4KZ/YmzXhQ+o7TAnyx71I/4 pGUNRiYvba4c+xPt8OXFAY8b3jBvT5Ozz9j56Whr7VOtNgY7AKnvRNfSTTHMeaUrLKcQjh XVZXTZFeVd0zN4nYmuiJvMYA+uCpWQ+ZBiieqxXmjn4pSf1EoLa4Hnpx1RH9fkBEpqQsFl aeNAl0p5eUhxfkIyXwsz0rGggX7SEMRUx0wfDW+0nOCfh8OKTB6tU05pokx36WTmqBH2IC 5n81BC+siq3VKsYYvr7DF4y7SHHrx/EVpGwmIR5l4xSJLVKEAZOpfou58Y+uFw== ARC-Authentication-Results: i=1; mx1.freebsd.org; none Received: from gitrepo.freebsd.org (gitrepo.freebsd.org [IPv6:2610:1c1:1:6068::e6a:5]) (using TLSv1.3 with cipher TLS_AES_256_GCM_SHA384 (256/256 bits) key-exchange X25519 server-signature RSA-PSS (4096 bits) server-digest SHA256) (Client did not present a certificate) by mxrelay.nyi.freebsd.org (Postfix) with ESMTPS id 4ZrTy25RlKz10Kt; Mon, 05 May 2025 05:06:02 +0000 (UTC) (envelope-from git@FreeBSD.org) Received: from gitrepo.freebsd.org ([127.0.1.44]) by gitrepo.freebsd.org (8.18.1/8.18.1) with ESMTP id 5455628W013589; Mon, 5 May 2025 05:06:02 GMT (envelope-from git@gitrepo.freebsd.org) Received: (from git@localhost) by gitrepo.freebsd.org (8.18.1/8.18.1/Submit) id 545562LH013587; Mon, 5 May 2025 05:06:02 GMT (envelope-from git) Date: Mon, 5 May 2025 05:06:02 GMT Message-Id: <202505050506.545562LH013587@gitrepo.freebsd.org> To: src-committers@FreeBSD.org, dev-commits-src-all@FreeBSD.org, dev-commits-src-main@FreeBSD.org From: Warner Losh Subject: git: f887d0215fb4 - main - msun: fix cbrt iterations from Newton to Halley method List-Id: Commit messages for all branches of the src repository List-Archive: https://lists.freebsd.org/archives/dev-commits-src-all List-Help: List-Post: List-Subscribe: List-Unsubscribe: X-BeenThere: dev-commits-src-all@freebsd.org Sender: owner-dev-commits-src-all@FreeBSD.org MIME-Version: 1.0 Content-Type: text/plain; charset=utf-8 Content-Transfer-Encoding: 8bit X-Git-Committer: imp X-Git-Repository: src X-Git-Refname: refs/heads/main X-Git-Reftype: branch X-Git-Commit: f887d0215fb48e682acccf4cb95f3794974e1a9d Auto-Submitted: auto-generated The branch main has been updated by imp: URL: https://cgit.FreeBSD.org/src/commit/?id=f887d0215fb48e682acccf4cb95f3794974e1a9d commit f887d0215fb48e682acccf4cb95f3794974e1a9d Author: Clément Bœsch AuthorDate: 2025-05-01 17:19:36 +0000 Commit: Warner Losh CommitDate: 2025-05-05 04:52:49 +0000 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ₙ) = = 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 Reviewed by: imp Pull Request: https://github.com/freebsd/freebsd-src/pull/1684 --- lib/msun/src/s_cbrt.c | 4 ++-- lib/msun/src/s_cbrtf.c | 4 ++-- lib/msun/src/s_cbrtl.c | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/lib/msun/src/s_cbrt.c b/lib/msun/src/s_cbrt.c index 6bf84243adcd..568a36545216 100644 --- a/lib/msun/src/s_cbrt.c +++ b/lib/msun/src/s_cbrt.c @@ -90,7 +90,7 @@ 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. @@ -99,7 +99,7 @@ cbrt(double x) u.bits=(u.bits+0x80000000)&0xffffffffc0000000ULL; t=u.value; - /* 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/lib/msun/src/s_cbrtf.c b/lib/msun/src/s_cbrtf.c index a225d3edb982..c69e0fa5be12 100644 --- a/lib/msun/src/s_cbrtf.c +++ b/lib/msun/src/s_cbrtf.c @@ -50,7 +50,7 @@ cbrtf(float x) SET_FLOAT_WORD(t,sign|(hx/3+B1)); /* - * 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. */ @@ -59,7 +59,7 @@ cbrtf(float x) T=T*((double)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/lib/msun/src/s_cbrtl.c b/lib/msun/src/s_cbrtl.c index f1950e2d4cef..ff527cc5e5e7 100644 --- a/lib/msun/src/s_cbrtl.c +++ b/lib/msun/src/s_cbrtl.c @@ -126,7 +126,7 @@ 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 */