1 /* mpfr_root -- kth root.
3 Copyright 2005-2015 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 /* The computation of y = x^(1/k) is done as follows:
28 Let x = sign * m * 2^(k*e) where m is an integer
30 with 2^(k*(n-1)) <= m < 2^(k*n) where n = PREC(y)
32 and m = s^k + r where 0 <= r and m < (s+1)^k
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
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.
42 mpfr_root (mpfr_ptr y
, mpfr_srcptr x
, unsigned long k
, mpfr_rnd_t rnd_mode
)
46 mpfr_prec_t n
, size_m
, tmp
;
47 int inexact
, negative
;
48 MPFR_SAVE_EXPO_DECL (expo
);
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
));
56 if (MPFR_UNLIKELY (k
<= 1))
58 if (k
< 1) /* k==0 => y=x^(1/0)=x^(+Inf) */
60 /* For 0 <= x < 1 => +0.
66 if (MPFR_IS_NEG (x
) && !MPFR_IS_ZERO (x
))
71 inexact
= mpfr_cmp (x
, __gmpfr_one
);
73 return mpfr_set_ui (y
, 1, rnd_mode
); /* 1 may be Out of Range */
75 return mpfr_set_ui (y
, 0, rnd_mode
); /* 0+ */
87 else /* y =x^(1/1)=x */
88 return mpfr_set (y
, x
, rnd_mode
);
92 else if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
96 MPFR_SET_NAN (y
); /* NaN^(1/k) = NaN */
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 */
103 if (MPFR_IS_NEG(x
) && (k
% 2 == 0))
109 MPFR_SET_SAME_SIGN (y
, x
);
112 else /* x is necessarily 0: (+0)^(1/k) = +0
115 MPFR_ASSERTD (MPFR_IS_ZERO (x
));
117 MPFR_SET_SAME_SIGN (y
, x
);
122 /* Returns NAN for x < 0 and k even */
123 else if (MPFR_IS_NEG (x
) && (k
% 2 == 0))
130 MPFR_SAVE_EXPO_MARK (expo
);
133 e
= mpfr_get_z_2exp (m
, x
); /* x = m * 2^e */
134 if ((negative
= MPFR_IS_NEG(x
)))
136 r
= e
% (mpfr_exp_t
) k
;
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 */
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
);
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 */
151 sh
= (k
* (mpfr_exp_t
) n
- (mpfr_exp_t
) size_m
- r
) / k
;
155 mpz_mul_2exp (m
, m
, sh
);
160 mpz_mul_2exp (m
, m
, r
);
164 /* invariant: x = m*2^e, with e divisible by k */
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;
170 MPFR_MPZ_SIZEINBASE2 (tmp
, m
);
172 if (sh
> 0) /* we have to flush to 0 the last sh bits from m */
174 inexact
= inexact
|| ((mpfr_exp_t
) mpz_scan1 (m
, 0) < sh
);
175 mpz_fdiv_q_2exp (m
, m
, sh
);
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);
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
);
198 MPFR_CHANGE_SIGN (y
);
203 MPFR_SAVE_EXPO_FREE (expo
);
204 return mpfr_check_range (y
, inexact
, rnd_mode
);