contrib/mpfr/src/root.c

   1 /* mpfr_root -- kth root.
   2
   3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   4 Contributed by the Arenaire 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
  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 }