summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/cbrt.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/cbrt.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/cbrt.c104
1 files changed, 64 insertions, 40 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/cbrt.c b/Build/source/libs/mpfr/mpfr-src/src/cbrt.c
index 50301761348..85f6cd2fb91 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/cbrt.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/cbrt.c
@@ -1,6 +1,6 @@
/* mpfr_cbrt -- cube root function.
-Copyright 2002-2019 Free Software Foundation, Inc.
+Copyright 2002-2020 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
@@ -23,31 +23,35 @@ 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
+ Let n = PREC(y), or PREC(y) + 1 if the rounding mode is MPFR_RNDN.
+ We seek to compute an integer cube root in precision n and the
+ associated inexact bit (non-zero iff the remainder is non-zero).
- with 2^(3n-3) <= m < 2^(3n) where n = PREC(y)
+ Let us write x, possibly truncated, under the form sign * m * 2^(3*e)
+ where m is an integer such that 2^(3n-3) <= m < 2^(3n), i.e. m has
+ between 3n-2 and 3n bits.
- and m = s^3 + r where 0 <= r and m < (s+1)^3
+ Let s be the integer cube root of m, i.e. the maximum integer such that
+ m = s^3 + t with t >= 0. Thus 2^(n-1) <= s < 2^n, i.e. s has n bits.
- we want that s has n bits i.e. s >= 2^(n-1), or m >= 2^(3n-3)
- i.e. m must have at least 3n-2 bits
+ Then |x|^(1/3) = s * 2^e or (s+1) * 2^e depending on the rounding mode,
+ the sign, and whether s is "inexact" (i.e. t > 0 or the truncation of x
+ was not equal to x).
- then x^(1/3) = s * 2^e if r=0
- x^(1/3) = (s+1) * 2^e if round up
- x^(1/3) = (s-1) * 2^e if round down
- x^(1/3) = s * 2^e if nearest and r < 3/2*s^2+3/4*s+1/8
- (s+1) * 2^e otherwise
- */
+ Note: The truncation of x was allowed because any breakpoint has n bits
+ and its cube has at most 3n bits. Thus the truncation of x cannot yield
+ a cube root below RNDZ(x^(1/3)) in precision n. [TODO: add details.]
+*/
int
mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
mpz_t m;
- mpfr_exp_t e, r, sh;
- mpfr_prec_t n, size_m, tmp;
- int inexact, negative;
+ mpfr_exp_t e, d, sh;
+ mpfr_prec_t n, size_m;
+ int inexact, inexact2, negative, r;
MPFR_SAVE_EXPO_DECL (expo);
MPFR_LOG_FUNC (
@@ -89,40 +93,55 @@ mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
r = e % 3;
if (r < 0)
r += 3;
+ MPFR_ASSERTD (r >= 0 && r < 3 && (e - r) % 3 == 0);
+
/* x = (m*2^r) * 2^(e-r) = (m*2^r) * 2^(3*q) */
+ MPFR_LOG_MSG (("e=%" MPFR_EXP_FSPEC "d r=%d\n", (mpfr_eexp_t) e, r));
+
MPFR_MPZ_SIZEINBASE2 (size_m, m);
n = MPFR_PREC (y) + (rnd_mode == MPFR_RNDN);
- /* we want 3*n-2 <= size_m + 3*sh + r <= 3*n
- i.e. 3*sh + size_m + r <= 3*n */
- sh = (3 * (mpfr_exp_t) n - (mpfr_exp_t) size_m - r) / 3;
- sh = 3 * sh + r;
- if (sh >= 0)
+ /* We will need to multiply m by 2^(r'), truncated if r' < 0, and
+ subtract r' from e, so that m has between 3n-2 and 3n bits and
+ e becomes a multiple of 3.
+ Since r = e % 3, we write r' = 3 * sh + r.
+ We want 3 * n - 2 <= size_m + 3 * sh + r <= 3 * n.
+ Let d = 3 * n - size_m - r. Thus we want 0 <= d - 3 * sh <= 2,
+ i.e. sh = floor(d/3). */
+ d = 3 * (mpfr_exp_t) n - (mpfr_exp_t) size_m - r;
+ sh = d >= 0 ? d / 3 : - ((2 - d) / 3); /* floor(d/3) */
+ r += 3 * sh; /* denoted r' above */
+
+ e -= r;
+ MPFR_ASSERTD (e % 3 == 0);
+ e /= 3;
+
+ inexact = 0;
+
+ if (r > 0)
{
- mpz_mul_2exp (m, m, sh);
- e = e - sh;
+ mpz_mul_2exp (m, m, r);
}
- else if (r > 0)
+ else if (r < 0)
{
- mpz_mul_2exp (m, m, r);
- e = e - r;
+ r = -r;
+ inexact = mpz_scan1 (m, 0) < r;
+ mpz_fdiv_q_2exp (m, m, r);
}
- /* invariant: x = m*2^e, with e divisible by 3 */
-
/* we reuse the variable m to store the cube root, since it is not needed
any more: we just need to know if the root is exact */
- inexact = mpz_root (m, m, 3) == 0;
+ inexact = ! mpz_root (m, m, 3) || inexact;
- MPFR_MPZ_SIZEINBASE2 (tmp, m);
- sh = tmp - n;
- if (sh > 0) /* we have to flush to 0 the last sh bits from m */
- {
- inexact = inexact || ((mpfr_exp_t) mpz_scan1 (m, 0) < sh);
- mpz_fdiv_q_2exp (m, m, sh);
- e += 3 * sh;
- }
+#if MPFR_WANT_ASSERT > 0
+ {
+ mpfr_prec_t tmp;
+
+ MPFR_MPZ_SIZEINBASE2 (tmp, m);
+ MPFR_ASSERTN (tmp == n);
+ }
+#endif
if (inexact)
{
@@ -130,7 +149,10 @@ mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
rnd_mode = MPFR_INVERT_RND (rnd_mode);
if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA
|| (rnd_mode == MPFR_RNDN && mpz_tstbit (m, 0)))
- inexact = 1, mpz_add_ui (m, m, 1);
+ {
+ inexact = 1;
+ mpz_add_ui (m, m, 1);
+ }
else
inexact = -1;
}
@@ -138,8 +160,10 @@ mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
/* either inexact is not zero, and the conversion is exact, i.e. inexact
is not changed; or inexact=0, and inexact is set only when
rnd_mode=MPFR_RNDN and bit (n+1) from m is 1 */
- inexact += mpfr_set_z (y, m, MPFR_RNDN);
- MPFR_SET_EXP (y, MPFR_GET_EXP (y) + e / 3);
+ inexact2 = mpfr_set_z (y, m, MPFR_RNDN);
+ MPFR_ASSERTD (inexact == 0 || inexact2 == 0);
+ inexact += inexact2;
+ MPFR_SET_EXP (y, MPFR_GET_EXP (y) + e);
if (negative)
{