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