1/* mpfr_root -- kth root. 2 3Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#define MPFR_NEED_LONGLONG_H 24#include "mpfr-impl.h" 25 26 /* The computation of y = x^(1/k) is done as follows: 27 28 Let x = sign * m * 2^(k*e) where m is an integer 29 30 with 2^(k*(n-1)) <= m < 2^(k*n) where n = PREC(y) 31 32 and m = s^k + r where 0 <= r and m < (s+1)^k 33 34 we want that s has n bits i.e. s >= 2^(n-1), or m >= 2^(k*(n-1)) 35 i.e. m must have at least k*(n-1)+1 bits 36 37 then, not taking into account the sign, the result will be 38 x^(1/k) = s * 2^e or (s+1) * 2^e according to the rounding mode. 39 */ 40 41int 42mpfr_root (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode) 43{ 44 mpz_t m; 45 mpfr_exp_t e, r, sh; 46 mpfr_prec_t n, size_m, tmp; 47 int inexact, negative; 48 MPFR_SAVE_EXPO_DECL (expo); 49 50 MPFR_LOG_FUNC 51 (("x[%Pu]=%.*Rg k=%lu rnd=%d", 52 mpfr_get_prec (x), mpfr_log_prec, x, k, rnd_mode), 53 ("y[%Pu]=%.*Rg inexact=%d", 54 mpfr_get_prec (y), mpfr_log_prec, y, inexact)); 55 56 if (MPFR_UNLIKELY (k <= 1)) 57 { 58 if (k < 1) /* k==0 => y=x^(1/0)=x^(+Inf) */ 59#if 0 60 /* For 0 <= x < 1 => +0. 61 For x = 1 => 1. 62 For x > 1, => +Inf. 63 For x < 0 => NaN. 64 */ 65 { 66 if (MPFR_IS_NEG (x) && !MPFR_IS_ZERO (x)) 67 { 68 MPFR_SET_NAN (y); 69 MPFR_RET_NAN; 70 } 71 inexact = mpfr_cmp (x, __gmpfr_one); 72 if (inexact == 0) 73 return mpfr_set_ui (y, 1, rnd_mode); /* 1 may be Out of Range */ 74 else if (inexact < 0) 75 return mpfr_set_ui (y, 0, rnd_mode); /* 0+ */ 76 else 77 { 78 mpfr_set_inf (y, 1); 79 return 0; 80 } 81 } 82#endif 83 { 84 MPFR_SET_NAN (y); 85 MPFR_RET_NAN; 86 } 87 else /* y =x^(1/1)=x */ 88 return mpfr_set (y, x, rnd_mode); 89 } 90 91 /* Singular values */ 92 else if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 93 { 94 if (MPFR_IS_NAN (x)) 95 { 96 MPFR_SET_NAN (y); /* NaN^(1/k) = NaN */ 97 MPFR_RET_NAN; 98 } 99 else if (MPFR_IS_INF (x)) /* +Inf^(1/k) = +Inf 100 -Inf^(1/k) = -Inf if k odd 101 -Inf^(1/k) = NaN if k even */ 102 { 103 if (MPFR_IS_NEG(x) && (k % 2 == 0)) 104 { 105 MPFR_SET_NAN (y); 106 MPFR_RET_NAN; 107 } 108 MPFR_SET_INF (y); 109 MPFR_SET_SAME_SIGN (y, x); 110 MPFR_RET (0); 111 } 112 else /* x is necessarily 0: (+0)^(1/k) = +0 113 (-0)^(1/k) = -0 */ 114 { 115 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 116 MPFR_SET_ZERO (y); 117 MPFR_SET_SAME_SIGN (y, x); 118 MPFR_RET (0); 119 } 120 } 121 122 /* Returns NAN for x < 0 and k even */ 123 else if (MPFR_IS_NEG (x) && (k % 2 == 0)) 124 { 125 MPFR_SET_NAN (y); 126 MPFR_RET_NAN; 127 } 128 129 /* General case */ 130 MPFR_SAVE_EXPO_MARK (expo); 131 mpz_init (m); 132 133 e = mpfr_get_z_2exp (m, x); /* x = m * 2^e */ 134 if ((negative = MPFR_IS_NEG(x))) 135 mpz_neg (m, m); 136 r = e % (mpfr_exp_t) k; 137 if (r < 0) 138 r += k; /* now r = e (mod k) with 0 <= e < r */ 139 /* x = (m*2^r) * 2^(e-r) where e-r is a multiple of k */ 140 141 MPFR_MPZ_SIZEINBASE2 (size_m, m); 142 /* for rounding to nearest, we want the round bit to be in the root */ 143 n = MPFR_PREC (y) + (rnd_mode == MPFR_RNDN); 144 145 /* we now multiply m by 2^(r+k*sh) so that root(m,k) will give 146 exactly n bits: we want k*(n-1)+1 <= size_m + k*sh + r <= k*n 147 i.e. sh = floor ((kn-size_m-r)/k) */ 148 if ((mpfr_exp_t) size_m + r > k * (mpfr_exp_t) n) 149 sh = 0; /* we already have too many bits */ 150 else 151 sh = (k * (mpfr_exp_t) n - (mpfr_exp_t) size_m - r) / k; 152 sh = k * sh + r; 153 if (sh >= 0) 154 { 155 mpz_mul_2exp (m, m, sh); 156 e = e - sh; 157 } 158 else if (r > 0) 159 { 160 mpz_mul_2exp (m, m, r); 161 e = e - r; 162 } 163 164 /* invariant: x = m*2^e, with e divisible by k */ 165 166 /* we reuse the variable m to store the kth root, since it is not needed 167 any more: we just need to know if the root is exact */ 168 inexact = mpz_root (m, m, k) == 0; 169 170 MPFR_MPZ_SIZEINBASE2 (tmp, m); 171 sh = tmp - n; 172 if (sh > 0) /* we have to flush to 0 the last sh bits from m */ 173 { 174 inexact = inexact || ((mpfr_exp_t) mpz_scan1 (m, 0) < sh); 175 mpz_fdiv_q_2exp (m, m, sh); 176 e += k * sh; 177 } 178 179 if (inexact) 180 { 181 if (negative) 182 rnd_mode = MPFR_INVERT_RND (rnd_mode); 183 if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA 184 || (rnd_mode == MPFR_RNDN && mpz_tstbit (m, 0))) 185 inexact = 1, mpz_add_ui (m, m, 1); 186 else 187 inexact = -1; 188 } 189 190 /* either inexact is not zero, and the conversion is exact, i.e. inexact 191 is not changed; or inexact=0, and inexact is set only when 192 rnd_mode=MPFR_RNDN and bit (n+1) from m is 1 */ 193 inexact += mpfr_set_z (y, m, MPFR_RNDN); 194 MPFR_SET_EXP (y, MPFR_GET_EXP (y) + e / (mpfr_exp_t) k); 195 196 if (negative) 197 { 198 MPFR_CHANGE_SIGN (y); 199 inexact = -inexact; 200 } 201 202 mpz_clear (m); 203 MPFR_SAVE_EXPO_FREE (expo); 204 return mpfr_check_range (y, inexact, rnd_mode); 205} 206