| /* mpfr_root -- kth root. |
| |
| Copyright 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc. |
| Contributed by the Arenaire and Cacao projects, INRIA. |
| |
| This file is part of the GNU MPFR Library. |
| |
| The GNU 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 GNU 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 GNU MPFR Library; see the file COPYING.LIB. If not, write to |
| the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, |
| MA 02110-1301, 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); |
| |
| 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 (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; |
| } |
| } |
| #endif |
| { |
| MPFR_SET_NAN (y); |
| MPFR_RET_NAN; |
| } |
| 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_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); |
| } |
| } |
| |
| /* Returns NAN for x < 0 and k even */ |
| else if (MPFR_IS_NEG (x) && (k % 2 == 0)) |
| { |
| MPFR_SET_NAN (y); |
| MPFR_RET_NAN; |
| } |
| |
| /* 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 kth 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); |
| } |