xref: /dragonfly/contrib/mpfr/src/root.c (revision ab6d115f)
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