summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorvlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4>2020-03-11 10:18:02 +0000
committervlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4>2020-03-11 10:18:02 +0000
commit4d289962d1f8909abb173d20cd9ea69aada5fecf (patch)
tree031ccfed078bd2dd261a260fbd5135885e8b7d4b
parent316737fac3d09c6b6b657732611ce91d937ef4b8 (diff)
downloadmpfr-4d289962d1f8909abb173d20cd9ea69aada5fecf.tar.gz
[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
-rw-r--r--src/cbrt.c17
1 files changed, 13 insertions, 4 deletions
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