From 4d289962d1f8909abb173d20cd9ea69aada5fecf Mon Sep 17 00:00:00 2001 From: vlefevre Date: Wed, 11 Mar 2020 10:18:02 +0000 Subject: [src/cbrt.c] Description of the algorithm: correction; added a TODO. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@13771 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/cbrt.c | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) (limited to 'src') diff --git a/src/cbrt.c b/src/cbrt.c index 58a35681a..0536df4fe 100644 --- a/src/cbrt.c +++ b/src/cbrt.c @@ -23,19 +23,28 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" -/* The computation of y = x^(1/3) is done as follows: +/* The computation of y = x^(1/3) is done as follows. - Let x = sign * m * 2^(3*e) where m is an integer >= 2^(3n-3) with - n = PREC(y), i.e. m has at least 3n-2 bits. + Let n = PREC(y), or PREC(y) + 1 if the rounding mode is MPFR_RNDN. + Let x = sign * m * 2^(3*e) where m is an integer >= 2^(3n-3), i.e. + m has at least 3n-2 bits. Let s be the integer cube root of m, i.e. the maximum integer such that m = s^3 + t with t >= 0. + TODO: Couldn't the size of m be fixed between 3n-2 and 3n? In the case + where the initial size of m is > 3n, if a discarded bit was non-zero, + this could be remembered for the sticky bit. Said otherwise, discard + 3k bits of the mpz_root argument instead of discarding k bits of its + result (integer cube root). + The constraint m >= 2^(3n-3) allows one to have sufficient precision for s: s >= 2^(n-1), i.e. s has at least n bits. FIXME: The description below is incorrect if s has more than n bits - (since n is the target precision). + (since n is the target precision). Also, I don't understand the + last case (nearest) since 1 was added to PREC(y) to make it like + the directed rounding modes. Then x^(1/3) = s * 2^e if t = 0 x^(1/3) = (s+1) * 2^e if round up -- cgit v1.2.1