diff -Naurd mpfr-3.1.3-a/PATCHES mpfr-3.1.3-b/PATCHES --- mpfr-3.1.3-a/PATCHES 2016-02-15 15:20:51.242696408 +0000 +++ mpfr-3.1.3-b/PATCHES 2016-02-15 15:20:51.306696441 +0000 @@ -0,0 +1 @@ +root diff -Naurd mpfr-3.1.3-a/VERSION mpfr-3.1.3-b/VERSION --- mpfr-3.1.3-a/VERSION 2016-02-15 15:20:16.922677881 +0000 +++ mpfr-3.1.3-b/VERSION 2016-02-15 15:20:51.306696441 +0000 @@ -1 +1 @@ -3.1.3-p11 +3.1.3-p12 diff -Naurd mpfr-3.1.3-a/src/mpfr.h mpfr-3.1.3-b/src/mpfr.h --- mpfr-3.1.3-a/src/mpfr.h 2016-02-15 15:20:16.922677881 +0000 +++ mpfr-3.1.3-b/src/mpfr.h 2016-02-15 15:20:51.302696439 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 3 #define MPFR_VERSION_MINOR 1 #define MPFR_VERSION_PATCHLEVEL 3 -#define MPFR_VERSION_STRING "3.1.3-p11" +#define MPFR_VERSION_STRING "3.1.3-p12" /* Macros dealing with MPFR VERSION */ #define MPFR_VERSION_NUM(a,b,c) (((a) << 16L) | ((b) << 8) | (c)) diff -Naurd mpfr-3.1.3-a/src/root.c mpfr-3.1.3-b/src/root.c --- mpfr-3.1.3-a/src/root.c 2015-06-19 19:55:10.000000000 +0000 +++ mpfr-3.1.3-b/src/root.c 2016-02-15 15:20:51.282696429 +0000 @@ -23,13 +23,15 @@ #define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" - /* The computation of y = x^(1/k) is done as follows: + /* The computation of y = x^(1/k) is done as follows, except for large + values of k, for which this would be inefficient or yield internal + integer overflows: Let x = sign * m * 2^(k*e) where m is an integer with 2^(k*(n-1)) <= m < 2^(k*n) where n = PREC(y) - and m = s^k + r where 0 <= r and m < (s+1)^k + and m = s^k + t where 0 <= t and m < (s+1)^k we want that s has n bits i.e. s >= 2^(n-1), or m >= 2^(k*(n-1)) i.e. m must have at least k*(n-1)+1 bits @@ -38,11 +40,15 @@ x^(1/k) = s * 2^e or (s+1) * 2^e according to the rounding mode. */ +static int +mpfr_root_aux (mpfr_ptr y, mpfr_srcptr x, unsigned long k, + mpfr_rnd_t rnd_mode); + int mpfr_root (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode) { mpz_t m; - mpfr_exp_t e, r, sh; + mpfr_exp_t e, r, sh, f; mpfr_prec_t n, size_m, tmp; int inexact, negative; MPFR_SAVE_EXPO_DECL (expo); @@ -55,50 +61,27 @@ if (MPFR_UNLIKELY (k <= 1)) { - if (k < 1) /* k==0 => y=x^(1/0)=x^(+Inf) */ -#if 0 - /* For 0 <= x < 1 => +0. - For x = 1 => 1. - For x > 1, => +Inf. - For x < 0 => NaN. - */ + if (k == 0) { - if (MPFR_IS_NEG (x) && !MPFR_IS_ZERO (x)) - { - MPFR_SET_NAN (y); - MPFR_RET_NAN; - } - inexact = mpfr_cmp (x, __gmpfr_one); - if (inexact == 0) - return mpfr_set_ui (y, 1, rnd_mode); /* 1 may be Out of Range */ - else if (inexact < 0) - return mpfr_set_ui (y, 0, rnd_mode); /* 0+ */ - else - { - mpfr_set_inf (y, 1); - return 0; - } + MPFR_SET_NAN (y); + MPFR_RET_NAN; } -#endif - { - MPFR_SET_NAN (y); - MPFR_RET_NAN; - } - else /* y =x^(1/1)=x */ + else /* y = x^(1/1) = x */ return mpfr_set (y, x, rnd_mode); } /* Singular values */ - else if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) + if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) { if (MPFR_IS_NAN (x)) { MPFR_SET_NAN (y); /* NaN^(1/k) = NaN */ MPFR_RET_NAN; } - else if (MPFR_IS_INF (x)) /* +Inf^(1/k) = +Inf - -Inf^(1/k) = -Inf if k odd - -Inf^(1/k) = NaN if k even */ + + if (MPFR_IS_INF (x)) /* +Inf^(1/k) = +Inf + -Inf^(1/k) = -Inf if k odd + -Inf^(1/k) = NaN if k even */ { if (MPFR_IS_NEG(x) && (k % 2 == 0)) { @@ -106,27 +89,31 @@ MPFR_RET_NAN; } MPFR_SET_INF (y); - MPFR_SET_SAME_SIGN (y, x); - MPFR_RET (0); } else /* x is necessarily 0: (+0)^(1/k) = +0 (-0)^(1/k) = -0 */ { MPFR_ASSERTD (MPFR_IS_ZERO (x)); MPFR_SET_ZERO (y); - MPFR_SET_SAME_SIGN (y, x); - MPFR_RET (0); } + MPFR_SET_SAME_SIGN (y, x); + MPFR_RET (0); } /* Returns NAN for x < 0 and k even */ - else if (MPFR_IS_NEG (x) && (k % 2 == 0)) + if (MPFR_UNLIKELY (MPFR_IS_NEG (x) && (k % 2 == 0))) { MPFR_SET_NAN (y); MPFR_RET_NAN; } /* General case */ + + /* For large k, use exp(log(x)/k). The threshold of 100 seems to be quite + good when the precision goes to infinity. */ + if (k > 100) + return mpfr_root_aux (y, x, k, rnd_mode); + MPFR_SAVE_EXPO_MARK (expo); mpz_init (m); @@ -135,31 +122,24 @@ mpz_neg (m, m); r = e % (mpfr_exp_t) k; if (r < 0) - r += k; /* now r = e (mod k) with 0 <= e < r */ + r += k; /* now r = e (mod k) with 0 <= r < k */ + MPFR_ASSERTD (0 <= r && r < k); /* x = (m*2^r) * 2^(e-r) where e-r is a multiple of k */ MPFR_MPZ_SIZEINBASE2 (size_m, m); /* for rounding to nearest, we want the round bit to be in the root */ n = MPFR_PREC (y) + (rnd_mode == MPFR_RNDN); - /* we now multiply m by 2^(r+k*sh) so that root(m,k) will give - exactly n bits: we want k*(n-1)+1 <= size_m + k*sh + r <= k*n - i.e. sh = floor ((kn-size_m-r)/k) */ - if ((mpfr_exp_t) size_m + r > k * (mpfr_exp_t) n) - sh = 0; /* we already have too many bits */ + /* we now multiply m by 2^sh so that root(m,k) will give + exactly n bits: we want k*(n-1)+1 <= size_m + sh <= k*n + i.e. sh = k*f + r with f = max(floor((k*n-size_m-r)/k),0) */ + if ((mpfr_exp_t) size_m + r >= k * (mpfr_exp_t) n) + f = 0; /* we already have too many bits */ else - sh = (k * (mpfr_exp_t) n - (mpfr_exp_t) size_m - r) / k; - sh = k * sh + r; - if (sh >= 0) - { - mpz_mul_2exp (m, m, sh); - e = e - sh; - } - else if (r > 0) - { - mpz_mul_2exp (m, m, r); - e = e - r; - } + f = (k * (mpfr_exp_t) n - (mpfr_exp_t) size_m - r) / k; + sh = k * f + r; + mpz_mul_2exp (m, m, sh); + e = e - sh; /* invariant: x = m*2^e, with e divisible by k */ @@ -203,3 +183,97 @@ MPFR_SAVE_EXPO_FREE (expo); return mpfr_check_range (y, inexact, rnd_mode); } + +/* Compute y <- x^(1/k) using exp(log(x)/k). + Assume all special cases have been eliminated before. + In the extended exponent range, overflows/underflows are not possible. + Assume x > 0, or x < 0 and k odd. +*/ +static int +mpfr_root_aux (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode) +{ + int inexact, exact_root = 0; + mpfr_prec_t w; /* working precision */ + mpfr_t absx, t; + MPFR_GROUP_DECL(group); + MPFR_TMP_DECL(marker); + MPFR_ZIV_DECL(loop); + MPFR_SAVE_EXPO_DECL (expo); + + MPFR_TMP_INIT_ABS (absx, x); + + MPFR_TMP_MARK(marker); + w = MPFR_PREC(y) + 10; + /* Take some guard bits to prepare for the 'expt' lost bits below. + If |x| < 2^k, then log|x| < k, thus taking log2(k) bits should be fine. */ + if (MPFR_GET_EXP(x) > 0) + w += MPFR_INT_CEIL_LOG2 (MPFR_GET_EXP(x)); + MPFR_GROUP_INIT_1(group, w, t); + MPFR_SAVE_EXPO_MARK (expo); + MPFR_ZIV_INIT (loop, w); + for (;;) + { + mpfr_exp_t expt; + unsigned int err; + + mpfr_log (t, absx, MPFR_RNDN); + /* t = log|x| * (1 + theta) with |theta| <= 2^(-w) */ + mpfr_div_ui (t, t, k, MPFR_RNDN); + expt = MPFR_GET_EXP (t); + /* t = log|x|/k * (1 + theta) + eps with |theta| <= 2^(-w) + and |eps| <= 1/2 ulp(t), thus the total error is bounded + by 1.5 * 2^(expt - w) */ + mpfr_exp (t, t, MPFR_RNDN); + /* t = |x|^(1/k) * exp(tau) * (1 + theta1) with + |tau| <= 1.5 * 2^(expt - w) and |theta1| <= 2^(-w). + For |tau| <= 0.5 we have |exp(tau)-1| < 4/3*tau, thus + for w >= expt + 2 we have: + t = |x|^(1/k) * (1 + 2^(expt+2)*theta2) * (1 + theta1) with + |theta1|, |theta2| <= 2^(-w). + If expt+2 > 0, as long as w >= 1, we have: + t = |x|^(1/k) * (1 + 2^(expt+3)*theta3) with |theta3| < 2^(-w). + For expt+2 = 0, we have: + t = |x|^(1/k) * (1 + 2^2*theta3) with |theta3| < 2^(-w). + Finally for expt+2 < 0 we have: + t = |x|^(1/k) * (1 + 2*theta3) with |theta3| < 2^(-w). + */ + err = (expt + 2 > 0) ? expt + 3 + : (expt + 2 == 0) ? 2 : 1; + /* now t = |x|^(1/k) * (1 + 2^(err-w)) thus the error is at most + 2^(EXP(t) - w + err) */ + if (MPFR_LIKELY (MPFR_CAN_ROUND(t, w - err, MPFR_PREC(y), rnd_mode))) + break; + + /* If we fail to round correctly, check for an exact result or a + midpoint result with MPFR_RNDN (regarded as hard-to-round in + all precisions in order to determine the ternary value). */ + { + mpfr_t z, zk; + + mpfr_init2 (z, MPFR_PREC(y) + (rnd_mode == MPFR_RNDN)); + mpfr_init2 (zk, MPFR_PREC(x)); + mpfr_set (z, t, MPFR_RNDN); + inexact = mpfr_pow_ui (zk, z, k, MPFR_RNDN); + exact_root = !inexact && mpfr_equal_p (zk, absx); + if (exact_root) /* z is the exact root, thus round z directly */ + inexact = mpfr_set4 (y, z, rnd_mode, MPFR_SIGN (x)); + mpfr_clear (zk); + mpfr_clear (z); + if (exact_root) + break; + } + + MPFR_ZIV_NEXT (loop, w); + MPFR_GROUP_REPREC_1(group, w, t); + } + MPFR_ZIV_FREE (loop); + + if (!exact_root) + inexact = mpfr_set4 (y, t, rnd_mode, MPFR_SIGN (x)); + + MPFR_GROUP_CLEAR(group); + MPFR_TMP_FREE(marker); + MPFR_SAVE_EXPO_FREE (expo); + + return mpfr_check_range (y, inexact, rnd_mode); +} diff -Naurd mpfr-3.1.3-a/src/version.c mpfr-3.1.3-b/src/version.c --- mpfr-3.1.3-a/src/version.c 2016-02-15 15:20:16.922677881 +0000 +++ mpfr-3.1.3-b/src/version.c 2016-02-15 15:20:51.306696441 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "3.1.3-p11"; + return "3.1.3-p12"; } diff -Naurd mpfr-3.1.3-a/tests/troot.c mpfr-3.1.3-b/tests/troot.c --- mpfr-3.1.3-a/tests/troot.c 2015-06-19 19:55:10.000000000 +0000 +++ mpfr-3.1.3-b/tests/troot.c 2016-02-15 15:20:51.282696429 +0000 @@ -25,6 +25,19 @@ #include "mpfr-test.h" +#define DEFN(N) \ + static int root##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) \ + { return mpfr_root (y, x, N, rnd); } \ + static int pow##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) \ + { return mpfr_pow_ui (y, x, N, rnd); } + +DEFN(2) +DEFN(3) +DEFN(4) +DEFN(5) +DEFN(17) +DEFN(120) + static void special (void) { @@ -52,7 +65,7 @@ exit (1); } - /* root(-Inf, 17) = -Inf */ + /* root(-Inf, 17) = -Inf */ mpfr_set_inf (x, -1); mpfr_root (y, x, 17, MPFR_RNDN); if (!mpfr_inf_p (y) || mpfr_sgn (y) > 0) @@ -69,7 +82,7 @@ exit (1); } - /* root(+/-0) = +/-0 */ + /* root(+/-0, k) = +/-0 for k > 0 */ mpfr_set_ui (x, 0, MPFR_RNDN); mpfr_root (y, x, 17, MPFR_RNDN); if (mpfr_cmp_ui (y, 0) || mpfr_sgn (y) < 0) @@ -190,64 +203,39 @@ i = mpfr_root (y, x, 1, MPFR_RNDN); if (mpfr_cmp_ui (x, 17) || i != 0) { - printf ("Error in root (17^(1/1))\n"); + printf ("Error in root for 17^(1/1)\n"); exit (1); } -#if 0 - /* Check for k == 0: - For 0 <= x < 1 => +0. - For x = 1 => 1. - For x > 1, => +Inf. - For x < 0 => NaN. */ - i = mpfr_root (y, x, 0, MPFR_RNDN); - if (!MPFR_IS_INF (y) || !MPFR_IS_POS (y) || i != 0) - { - printf ("Error in root 17^(1/0)\n"); - exit (1); - } - mpfr_set_ui (x, 1, MPFR_RNDN); - i = mpfr_root (y, x, 0, MPFR_RNDN); - if (mpfr_cmp_ui (y, 1) || i != 0) - { - printf ("Error in root 1^(1/0)\n"); - exit (1); - } mpfr_set_ui (x, 0, MPFR_RNDN); i = mpfr_root (y, x, 0, MPFR_RNDN); - if (!MPFR_IS_ZERO (y) || !MPFR_IS_POS (y) || i != 0) - { - printf ("Error in root 0+^(1/0)\n"); - exit (1); - } - MPFR_CHANGE_SIGN (x); - i = mpfr_root (y, x, 0, MPFR_RNDN); - if (!MPFR_IS_ZERO (y) || !MPFR_IS_POS (y) || i != 0) + if (!MPFR_IS_NAN (y) || i != 0) { - printf ("Error in root 0-^(1/0)\n"); + printf ("Error in root for (+0)^(1/0)\n"); exit (1); } - mpfr_set_ui_2exp (x, 17, -5, MPFR_RNDD); + mpfr_neg (x, x, MPFR_RNDN); i = mpfr_root (y, x, 0, MPFR_RNDN); - if (!MPFR_IS_ZERO (y) || !MPFR_IS_POS (y) || i != 0) + if (!MPFR_IS_NAN (y) || i != 0) { - printf ("Error in root (17/2^5)^(1/0)\n"); + printf ("Error in root for (-0)^(1/0)\n"); exit (1); } -#endif - mpfr_set_ui (x, 0, MPFR_RNDN); + + mpfr_set_ui (x, 1, MPFR_RNDN); i = mpfr_root (y, x, 0, MPFR_RNDN); if (!MPFR_IS_NAN (y) || i != 0) { - printf ("Error in root 0+^(1/0)\n"); + printf ("Error in root for 1^(1/0)\n"); exit (1); } + /* Check for k==2 */ mpfr_set_si (x, -17, MPFR_RNDD); i = mpfr_root (y, x, 2, MPFR_RNDN); if (!MPFR_IS_NAN (y) || i != 0) { - printf ("Error in root (-17)^(1/2)\n"); + printf ("Error in root for (-17)^(1/2)\n"); exit (1); } @@ -255,11 +243,168 @@ mpfr_clear (y); } +/* https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=812779 + * https://bugzilla.gnome.org/show_bug.cgi?id=756960 + * is a GNOME Calculator bug (mpfr_root applied on a negative integer, + * which is converted to an unsigned integer), but the strange result + * is also due to a bug in MPFR. + */ +static void +bigint (void) +{ + mpfr_t x, y; + + mpfr_inits2 (64, x, y, (mpfr_ptr) 0); + + mpfr_set_ui (x, 10, MPFR_RNDN); + if (sizeof (unsigned long) * CHAR_BIT == 64) + { + mpfr_root (x, x, ULONG_MAX, MPFR_RNDN); + mpfr_set_ui_2exp (y, 1, -63, MPFR_RNDN); + mpfr_add_ui (y, y, 1, MPFR_RNDN); + if (! mpfr_equal_p (x, y)) + { + printf ("Error in bigint for ULONG_MAX\n"); + printf ("Expected "); + mpfr_dump (y); + printf ("Got "); + mpfr_dump (x); + exit (1); + } + } + + mpfr_set_ui (x, 10, MPFR_RNDN); + mpfr_root (x, x, 1234567890, MPFR_RNDN); + mpfr_set_str_binary (y, + "1.00000000000000000000000000001000000000101011000101000110010001"); + if (! mpfr_equal_p (x, y)) + { + printf ("Error in bigint for 1234567890\n"); + printf ("Expected "); + mpfr_dump (y); + printf ("Got "); + mpfr_dump (x); + exit (1); + } + + mpfr_clears (x, y, (mpfr_ptr) 0); +} + #define TEST_FUNCTION mpfr_root #define INTEGER_TYPE unsigned long -#define INT_RAND_FUNCTION() (INTEGER_TYPE) (randlimb () % 3 +2) +#define INT_RAND_FUNCTION() \ + (INTEGER_TYPE) (randlimb () & 1 ? randlimb () : randlimb () % 3 + 2) #include "tgeneric_ui.c" +static void +exact_powers (unsigned long bmax, unsigned long kmax) +{ + long b, k; + mpz_t z; + mpfr_t x, y; + int inex, neg; + + mpz_init (z); + for (b = 2; b <= bmax; b++) + for (k = 1; k <= kmax; k++) + { + mpz_ui_pow_ui (z, b, k); + mpfr_init2 (x, mpz_sizeinbase (z, 2)); + mpfr_set_ui (x, b, MPFR_RNDN); + mpfr_pow_ui (x, x, k, MPFR_RNDN); + mpz_set_ui (z, b); + mpfr_init2 (y, mpz_sizeinbase (z, 2)); + for (neg = 0; neg <= 1; neg++) + { + inex = mpfr_root (y, x, k, MPFR_RNDN); + if (inex != 0) + { + printf ("Error in exact_powers, b=%ld, k=%ld\n", b, k); + printf ("Expected inex=0, got %d\n", inex); + exit (1); + } + if (neg && (k & 1) == 0) + { + if (!MPFR_IS_NAN (y)) + { + printf ("Error in exact_powers, b=%ld, k=%ld\n", b, k); + printf ("Expected y=NaN\n"); + printf ("Got "); + mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN); + printf ("\n"); + exit (1); + } + } + else if (MPFR_IS_NAN (y) || mpfr_cmp_si (y, b) != 0) + { + printf ("Error in exact_powers, b=%ld, k=%ld\n", b, k); + printf ("Expected y=%ld\n", b); + printf ("Got "); + mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN); + printf ("\n"); + exit (1); + } + mpfr_neg (x, x, MPFR_RNDN); + b = -b; + } + mpfr_clear (x); + mpfr_clear (y); + } + mpz_clear (z); +} + +/* Compare root(x,2^h) with pow(x,2^(-h)). */ +static void +cmp_pow (void) +{ + mpfr_t x, y1, y2; + int h; + + mpfr_inits2 (128, x, y1, y2, (mpfr_ptr) 0); + + for (h = 1; h < sizeof (unsigned long) * CHAR_BIT; h++) + { + unsigned long k = (unsigned long) 1 << h; + int i; + + for (i = 0; i < 10; i++) + { + mpfr_rnd_t rnd; + unsigned int flags1, flags2; + int inex1, inex2; + + tests_default_random (x, 0, __gmpfr_emin, __gmpfr_emax, 1); + rnd = RND_RAND (); + mpfr_set_ui_2exp (y1, 1, -h, MPFR_RNDN); + mpfr_clear_flags (); + inex1 = mpfr_pow (y1, x, y1, rnd); + flags1 = __gmpfr_flags; + mpfr_clear_flags (); + inex2 = mpfr_root (y2, x, k, rnd); + flags2 = __gmpfr_flags; + if (!(mpfr_equal_p (y1, y2) && SAME_SIGN (inex1, inex2) && + flags1 == flags2)) + { + printf ("Error in cmp_pow on h=%d, i=%d, rnd=%s\n", + h, i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); + printf ("x = "); + mpfr_dump (x); + printf ("pow = "); + mpfr_dump (y1); + printf ("with inex = %d, flags =", inex1); + flags_out (flags1); + printf ("root = "); + mpfr_dump (y2); + printf ("with inex = %d, flags =", inex2); + flags_out (flags2); + exit (1); + } + } + } + + mpfr_clears (x, y1, y2, (mpfr_ptr) 0); +} + int main (void) { @@ -270,7 +415,10 @@ tests_start_mpfr (); + exact_powers (3, 1000); special (); + bigint (); + cmp_pow (); mpfr_init (x); @@ -329,6 +477,13 @@ test_generic_ui (2, 200, 30); + bad_cases (root2, pow2, "mpfr_root[2]", 8, -256, 255, 4, 128, 800, 40); + bad_cases (root3, pow3, "mpfr_root[3]", 8, -256, 255, 4, 128, 800, 40); + bad_cases (root4, pow4, "mpfr_root[4]", 8, -256, 255, 4, 128, 800, 40); + bad_cases (root5, pow5, "mpfr_root[5]", 8, -256, 255, 4, 128, 800, 40); + bad_cases (root17, pow17, "mpfr_root[17]", 8, -256, 255, 4, 128, 800, 40); + bad_cases (root120, pow120, "mpfr_root[120]", 8, -256, 255, 4, 128, 800, 40); + tests_end_mpfr (); return 0; }