diff options
| -rw-r--r-- | cbrt.c | 2 | ||||
| -rw-r--r-- | mpfr.h | 1 | ||||
| -rw-r--r-- | mpfr.texi | 9 | ||||
| -rw-r--r-- | root.c | 157 | ||||
| -rw-r--r-- | tests/tgeneric.c | 18 | ||||
| -rw-r--r-- | tests/troot.c | 224 |
6 files changed, 405 insertions, 6 deletions
@@ -65,7 +65,7 @@ mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) MPFR_RET (0); } /* case 0: cbrt(+/- 0) = +/- 0 */ - else /* x is necessaryly 0 */ + else /* x is necessarily 0 */ { MPFR_ASSERTD (MPFR_IS_ZERO (x)); MPFR_SET_ZERO (y); @@ -490,6 +490,7 @@ __MPFR_DECLSPEC int mpfr_hypot _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, __MPFR_DECLSPEC int mpfr_erf _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_erfc _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_cbrt _MPFR_PROTO ((mpfr_ptr,mpfr_srcptr,mpfr_rnd_t)); +__MPFR_DECLSPEC int mpfr_root _MPFR_PROTO ((mpfr_ptr,mpfr_srcptr,unsigned long,mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_gamma _MPFR_PROTO((mpfr_ptr,mpfr_srcptr,mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_zeta _MPFR_PROTO ((mpfr_ptr,mpfr_srcptr,mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_fac_ui _MPFR_PROTO ((mpfr_ptr, unsigned long int, @@ -1037,8 +1037,13 @@ Set @var{rop} to NaN if @var{op} is negative. @end deftypefun @deftypefun int mpfr_cbrt (mpfr_t @var{rop}, mpfr_t @var{op}, mp_rnd_t @var{rnd}) -Set @var{rop} to the cubic root (defined over the real numbers) of @var{op} -rounded in the direction @var{rnd}. +@deftypefunx int mpfr_root (mpfr_t @var{rop}, mpfr_t @var{op}, unsigned long int @var{k}, mp_rnd_t @var{rnd}) +Set @var{rop} to the cubic root (resp. the @var{k}th root) +of @var{op} rounded in the direction @var{rnd}. +An odd (resp. even) root of a negative number (including -Inf) +returns a negative number (resp. NaN). +The @var{k}th root of -0 is defined to be -0, +whatever the parity of @var{k} (this may change in the future). @end deftypefun @deftypefun int mpfr_pow (mpfr_t @var{rop}, mpfr_t @var{op1}, mpfr_t @var{op2}, mp_rnd_t @var{rnd}) @@ -0,0 +1,157 @@ +/* mpfr_root -- kth root. + +Copyright 2005 Free Software Foundation. +Contributed by the Spaces project, INRIA Lorraine. + +This file is part of the MPFR Library. + +The MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your +option) any later version. + +The MPFR Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#define MPFR_NEED_LONGLONG_H +#include "mpfr-impl.h" + + /* The computation of y = x^(1/k) is done as follows: + + 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 + + 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 + + then, not taking into account the sign, the result will be + x^(1/k) = s * 2^e or (s+1) * 2^e according to the rounding mode. + */ + +int +mpfr_root (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mp_rnd_t rnd_mode) +{ + mpz_t m; + mp_exp_t e, r, sh; + mp_prec_t n, size_m, tmp; + int inexact, negative; + MPFR_SAVE_EXPO_DECL (expo); + + /* special values */ + 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_NEG(x) && (k % 2 == 0)) + { + MPFR_SET_NAN (y); + 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); + } + } + + /* General case */ + MPFR_SAVE_EXPO_MARK (expo); + mpz_init (m); + + e = mpfr_get_z_exp (m, x); /* x = m * 2^e */ + if ((negative = MPFR_IS_NEG(x))) + mpz_neg (m, m); + r = e % (mp_exp_t) k; + if (r < 0) + r += k; /* now r = e (mod k) with 0 <= e < r */ + /* 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 == GMP_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 ((mp_exp_t) size_m + r > k * (mp_exp_t) n) + sh = 0; /* we already have too many bits */ + else + sh = (k * (mp_exp_t) n - (mp_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; + } + + /* invariant: x = m*2^e, with e divisible by k */ + + /* 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, k) == 0; + + 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 || ((mp_exp_t) mpz_scan1 (m, 0) < sh); + mpz_div_2exp (m, m, sh); + e += k * sh; + } + + if (inexact) + { + if (negative) + rnd_mode = MPFR_INVERT_RND (rnd_mode); + if (rnd_mode == GMP_RNDU + || (rnd_mode == GMP_RNDN && mpz_tstbit (m, 0))) + inexact = 1, mpz_add_ui (m, m, 1); + else + inexact = -1; + } + + /* 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=GMP_RNDN and bit (n+1) from m is 1 */ + inexact += mpfr_set_z (y, m, GMP_RNDN); + MPFR_SET_EXP (y, MPFR_GET_EXP (y) + e / (mp_exp_t) k); + + if (negative) + { + MPFR_CHANGE_SIGN (y); + inexact = -inexact; + } + + mpz_clear (m); + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_check_range (y, inexact, rnd_mode); +} diff --git a/tests/tgeneric.c b/tests/tgeneric.c index 846fb347e..0648440eb 100644 --- a/tests/tgeneric.c +++ b/tests/tgeneric.c @@ -20,9 +20,14 @@ the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ /* define TWO_ARGS for two-argument functions like mpfr_pow */ +/* define TWO_ARGS_UI for two-argument functions like mpfr_root */ static void -test_generic (mp_prec_t p0, mp_prec_t p1, unsigned int N) +test_generic (mp_prec_t p0, mp_prec_t p1, +#ifdef TWO_ARGS_UI + unsigned long u, +#endif +unsigned int N) { mp_prec_t prec, yprec; mpfr_t x, y, z, t; @@ -64,7 +69,7 @@ test_generic (mp_prec_t p0, mp_prec_t p1, unsigned int N) #endif rnd = (mp_rnd_t) RND_RAND (); mpfr_set_prec (y, yprec); -#ifdef TWO_ARGS +#if (defined(TWO_ARGS) || defined(TWO_ARGS_UI)) compare = TEST_FUNCTION (y, x, u, rnd); #else compare = TEST_FUNCTION (y, x, rnd); @@ -72,7 +77,7 @@ test_generic (mp_prec_t p0, mp_prec_t p1, unsigned int N) if (mpfr_can_round (y, yprec, rnd, rnd, prec)) { mpfr_set (t, y, rnd); -#ifdef TWO_ARGS +#if (defined(TWO_ARGS) || defined(TWO_ARGS_UI)) inexact = TEST_FUNCTION (z, x, u, rnd); #else inexact = TEST_FUNCTION (z, x, rnd); @@ -85,6 +90,9 @@ test_generic (mp_prec_t p0, mp_prec_t p1, unsigned int N) printf ("\nu="); mpfr_out_str (stdout, 2, prec, u, GMP_RNDN); #endif +#ifdef TWO_ARGS + printf ("\nu=%lu", u); +#endif printf (" prec=%u rnd_mode=%s\n", (unsigned) prec, mpfr_print_rnd_mode (rnd)); printf ("got "); @@ -115,6 +123,9 @@ test_generic (mp_prec_t p0, mp_prec_t p1, unsigned int N) #ifdef TWO_ARGS printf ("u="); mpfr_print_binary (u); puts (""); #endif +#ifdef TWO_ARGS_UI + printf ("u=%lu\n", u); +#endif printf ("y="); mpfr_print_binary (y); puts (""); printf ("t="); mpfr_print_binary (t); puts (""); exit (1); @@ -134,5 +145,6 @@ test_generic (mp_prec_t p0, mp_prec_t p1, unsigned int N) #undef RAND_FUNCTION #undef TWO_ARGS +#undef TWO_ARGS_UI #undef TEST_FUNCTION #undef test_generic diff --git a/tests/troot.c b/tests/troot.c new file mode 100644 index 000000000..525f74370 --- /dev/null +++ b/tests/troot.c @@ -0,0 +1,224 @@ +/* Test file for mpfr_root. + +Copyright 2005 Free Software Foundation, Inc. + +This file is part of the MPFR Library. + +The MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your +option) any later version. + +The MPFR Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include <stdlib.h> + +#include "mpfr-test.h" + +static void +special (void) +{ + mpfr_t x, y; + + mpfr_init (x); + mpfr_init (y); + + /* root(NaN) = NaN */ + mpfr_set_nan (x); + mpfr_root (y, x, 17, GMP_RNDN); + if (!mpfr_nan_p (y)) + { + printf ("Error: root(NaN,17) <> NaN\n"); + exit (1); + } + + /* root(+Inf) = +Inf */ + mpfr_set_inf (x, 1); + mpfr_root (y, x, 42, GMP_RNDN); + if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0) + { + printf ("Error: root(+Inf,42) <> +Inf\n"); + exit (1); + } + + /* root(-Inf, 17) = -Inf */ + mpfr_set_inf (x, -1); + mpfr_root (y, x, 17, GMP_RNDN); + if (!mpfr_inf_p (y) || mpfr_sgn (y) > 0) + { + printf ("Error: root(-Inf,17) <> -Inf\n"); + exit (1); + } + /* root(-Inf, 42) = NaN */ + mpfr_set_inf (x, -1); + mpfr_root (y, x, 42, GMP_RNDN); + if (!mpfr_nan_p (y)) + { + printf ("Error: root(-Inf,42) <> -Inf\n"); + exit (1); + } + + /* root(+/-0) = +/-0 */ + mpfr_set_ui (x, 0, GMP_RNDN); + mpfr_root (y, x, 17, GMP_RNDN); + if (mpfr_cmp_ui (y, 0) || mpfr_sgn (y) < 0) + { + printf ("Error: root(+0,17) <> +0\n"); + exit (1); + } + mpfr_neg (x, x, GMP_RNDN); + mpfr_root (y, x, 42, GMP_RNDN); + if (mpfr_cmp_ui (y, 0) || mpfr_sgn (y) > 0) + { + printf ("Error: root(-0,42) <> -0\n"); + exit (1); + } + + mpfr_set_prec (x, 53); + mpfr_set_str (x, "8.39005285514734966412e-01", 10, GMP_RNDN); + mpfr_root (x, x, 3, GMP_RNDN); + if (mpfr_cmp_str1 (x, "9.43166207799662426048e-01")) + { + printf ("Error in root3 (1)\n"); + printf ("expected 9.43166207799662426048e-01\n"); + printf ("got "); + mpfr_dump (x); + exit (1); + } + + mpfr_set_prec (x, 32); + mpfr_set_prec (y, 32); + mpfr_set_str_binary (x, "0.10000100001100101001001001011001"); + mpfr_root (x, x, 3, GMP_RNDN); + mpfr_set_str_binary (y, "0.11001101011000100111000111111001"); + if (mpfr_cmp (x, y)) + { + printf ("Error in root3 (2)\n"); + exit (1); + } + + mpfr_set_prec (x, 32); + mpfr_set_prec (y, 32); + mpfr_set_str_binary (x, "-0.1100001110110000010101011001011"); + mpfr_root (x, x, 3, GMP_RNDD); + mpfr_set_str_binary (y, "-0.11101010000100100101000101011001"); + if (mpfr_cmp (x, y)) + { + printf ("Error in root3 (3)\n"); + exit (1); + } + + mpfr_set_prec (x, 82); + mpfr_set_prec (y, 27); + mpfr_set_str_binary (x, "0.1010001111011101011011000111001011001101100011110110010011011011011010011001100101e-7"); + mpfr_root (y, x, 3, GMP_RNDD); + mpfr_set_str_binary (x, "0.101011110001110001000100011E-2"); + if (mpfr_cmp (x, y)) + { + printf ("Error in root3 (4)\n"); + exit (1); + } + + mpfr_set_prec (x, 204); + mpfr_set_prec (y, 38); + mpfr_set_str_binary (x, "0.101000000001101000000001100111111011111001110110100001111000100110100111001101100111110001110001011011010110010011100101111001111100001010010100111011101100000011011000101100010000000011000101001010001001E-5"); + mpfr_root (y, x, 3, GMP_RNDD); + mpfr_set_str_binary (x, "0.10001001111010011011101000010110110010E-1"); + if (mpfr_cmp (x, y)) + { + printf ("Error in root3 (5)\n"); + exit (1); + } + + mpfr_clear (x); + mpfr_clear (y); +} + +#define TEST_FUNCTION mpfr_root +#define TWO_ARGS_UI +#include "tgeneric.c" + +int +main (void) +{ + mpfr_t x; + int r; + mp_prec_t p; + unsigned long k; + + MPFR_TEST_USE_RANDS (); + tests_start_mpfr (); + + special (); + + mpfr_init (x); + + for (p = 2; p < 100; p++) + { + mpfr_set_prec (x, p); + for (r = 0; r < GMP_RND_MAX; r++) + { + mpfr_set_ui (x, 1, GMP_RNDN); + k = 2 + randlimb () % 4; /* 2 <= k <= 5 */ + mpfr_root (x, x, k, (mp_rnd_t) r); + if (mpfr_cmp_ui (x, 1)) + { + printf ("Error in mpfr_root(%lu) for x=1, rnd=%s\ngot ", + k, mpfr_print_rnd_mode ((mp_rnd_t) r)); + mpfr_out_str (stdout, 2, 0, x, GMP_RNDN); + printf ("\n"); + exit (1); + } + mpfr_set_si (x, -1, GMP_RNDN); + if (k % 2) + { + mpfr_root (x, x, k, (mp_rnd_t) r); + if (mpfr_cmp_si (x, -1)) + { + printf ("Error in mpfr_root(%lu) for x=-1, rnd=%s\ngot ", + k, mpfr_print_rnd_mode ((mp_rnd_t) r)); + mpfr_out_str (stdout, 2, 0, x, GMP_RNDN); + printf ("\n"); + exit (1); + } + } + + if (p >= 5) + { + int i; + for (i = -12; i <= 12; i++) + { + mpfr_set_ui (x, 27, GMP_RNDN); + mpfr_mul_2si (x, x, 3*i, GMP_RNDN); + mpfr_root (x, x, 3, GMP_RNDN); + if (mpfr_cmp_si_2exp (x, 3, i)) + { + printf ("Error in mpfr_root(3) for " + "x = 27.0 * 2^(%d), rnd=%s\ngot ", + 3*i, mpfr_print_rnd_mode ((mp_rnd_t) r)); + mpfr_out_str (stdout, 2, 0, x, GMP_RNDN); + printf ("\ninstead of 3 * 2^(%d)\n", i); + exit (1); + } + } + } + } + } + mpfr_clear (x); + + test_generic (2, 200, 3, 10); + test_generic (2, 200, 4, 10); + test_generic (2, 200, 5, 10); + + tests_end_mpfr (); + return 0; +} |