1 /* mpfr_root -- kth root.
2
3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 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
41 int
mpfr_root(mpfr_ptr y,mpfr_srcptr x,unsigned long k,mpfr_rnd_t rnd_mode)42 mpfr_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