contrib/mpfr/root.c

   1 /* mpfr_root -- kth root.
   2
   3 Copyright 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao 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 2.1 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.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 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, mp_rnd_t rnd_mode)
  43 {
  44   mpz_t m;
  45   mp_exp_t e, r, sh;
  46   mp_prec_t n, size_m, tmp;
  47   int inexact, negative;
  48   MPFR_SAVE_EXPO_DECL (expo);
  49
  50   if (MPFR_UNLIKELY (k <= 1))
  51     {
  52       if (k < 1) /* k==0 => y=x^(1/0)=x^(+Inf) */
  53 #if 0
  54         /* For 0 <= x < 1 => +0.
  55            For x = 1      => 1.
  56            For x > 1,     => +Inf.
  57            For x < 0      => NaN.
  58         */
  59         {
  60           if (MPFR_IS_NEG (x) && !MPFR_IS_ZERO (x))
  61             {
  62               MPFR_SET_NAN (y);
  63               MPFR_RET_NAN;
  64             }
  65           inexact = mpfr_cmp (x, __gmpfr_one);
  66           if (inexact == 0)
  67             return mpfr_set_ui (y, 1, rnd_mode); /* 1 may be Out of Range */
  68           else if (inexact < 0)
  69             return mpfr_set_ui (y, 0, rnd_mode); /* 0+ */
  70           else
  71             {
  72               mpfr_set_inf (y, 1);
  73               return 0;
  74             }
  75         }
  76 #endif
  77       {
  78         MPFR_SET_NAN (y);
  79         MPFR_RET_NAN;
  80       }
  81       else /* y =x^(1/1)=x */
  82         return mpfr_set (y, x, rnd_mode);
  83     }
  84
  85   /* Singular values */
  86   else if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
  87     {
  88       if (MPFR_IS_NAN (x))
  89         {
  90           MPFR_SET_NAN (y); /* NaN^(1/k) = NaN */
  91           MPFR_RET_NAN;
  92         }
  93       else if (MPFR_IS_INF (x)) /* +Inf^(1/k) = +Inf
  94                                    -Inf^(1/k) = -Inf if k odd
  95                                    -Inf^(1/k) = NaN if k even */
  96         {
  97           if (MPFR_IS_NEG(x) && (k % 2 == 0))
  98             {
  99               MPFR_SET_NAN (y);
 100               MPFR_RET_NAN;
 101             }
 102           MPFR_SET_INF (y);
 103           MPFR_SET_SAME_SIGN (y, x);
 104           MPFR_RET (0);
 105         }
 106       else /* x is necessarily 0: (+0)^(1/k) = +0
 107                                   (-0)^(1/k) = -0 */
 108         {
 109           MPFR_ASSERTD (MPFR_IS_ZERO (x));
 110           MPFR_SET_ZERO (y);
 111           MPFR_SET_SAME_SIGN (y, x);
 112           MPFR_RET (0);
 113         }
 114     }
 115
 116   /* Returns NAN for x < 0 and k even */
 117   else if (MPFR_IS_NEG (x) && (k % 2 == 0))
 118     {
 119       MPFR_SET_NAN (y);
 120       MPFR_RET_NAN;
 121     }
 122
 123   /* General case */
 124   MPFR_SAVE_EXPO_MARK (expo);
 125   mpz_init (m);
 126
 127   e = mpfr_get_z_exp (m, x);                /* x = m * 2^e */
 128   if ((negative = MPFR_IS_NEG(x)))
 129     mpz_neg (m, m);
 130   r = e % (mp_exp_t) k;
 131   if (r < 0)
 132     r += k; /* now r = e (mod k) with 0 <= e < r */
 133   /* x = (m*2^r) * 2^(e-r) where e-r is a multiple of k */
 134
 135   MPFR_MPZ_SIZEINBASE2 (size_m, m);
 136   /* for rounding to nearest, we want the round bit to be in the root */
 137   n = MPFR_PREC (y) + (rnd_mode == GMP_RNDN);
 138
 139   /* we now multiply m by 2^(r+k*sh) so that root(m,k) will give
 140      exactly n bits: we want k*(n-1)+1 <= size_m + k*sh + r <= k*n
 141      i.e. sh = floor ((kn-size_m-r)/k) */
 142   if ((mp_exp_t) size_m + r > k * (mp_exp_t) n)
 143     sh = 0; /* we already have too many bits */
 144   else
 145     sh = (k * (mp_exp_t) n - (mp_exp_t) size_m - r) / k;
 146   sh = k * sh + r;
 147   if (sh >= 0)
 148     {
 149       mpz_mul_2exp (m, m, sh);
 150       e = e - sh;
 151     }
 152   else if (r > 0)
 153     {
 154       mpz_mul_2exp (m, m, r);
 155       e = e - r;
 156     }
 157
 158   /* invariant: x = m*2^e, with e divisible by k */
 159
 160   /* we reuse the variable m to store the kth root, since it is not needed
 161      any more: we just need to know if the root is exact */
 162   inexact = mpz_root (m, m, k) == 0;
 163
 164   MPFR_MPZ_SIZEINBASE2 (tmp, m);
 165   sh = tmp - n;
 166   if (sh > 0) /* we have to flush to 0 the last sh bits from m */
 167     {
 168       inexact = inexact || ((mp_exp_t) mpz_scan1 (m, 0) < sh);
 169       mpz_div_2exp (m, m, sh);
 170       e += k * sh;
 171     }
 172
 173   if (inexact)
 174     {
 175       if (negative)
 176         rnd_mode = MPFR_INVERT_RND (rnd_mode);
 177       if (rnd_mode == GMP_RNDU
 178           || (rnd_mode == GMP_RNDN && mpz_tstbit (m, 0)))
 179         inexact = 1, mpz_add_ui (m, m, 1);
 180       else
 181         inexact = -1;
 182     }
 183
 184   /* either inexact is not zero, and the conversion is exact, i.e. inexact
 185      is not changed; or inexact=0, and inexact is set only when
 186      rnd_mode=GMP_RNDN and bit (n+1) from m is 1 */
 187   inexact += mpfr_set_z (y, m, GMP_RNDN);
 188   MPFR_SET_EXP (y, MPFR_GET_EXP (y) + e / (mp_exp_t) k);
 189
 190   if (negative)
 191     {
 192       MPFR_CHANGE_SIGN (y);
 193       inexact = -inexact;
 194     }
 195
 196   mpz_clear (m);
 197   MPFR_SAVE_EXPO_FREE (expo);
 198   return mpfr_check_range (y, inexact, rnd_mode);
 199 }