contrib/mpfr/src/sqrt.c

   1 /* mpfr_sqrt -- square root of a floating-point number
   2
   3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 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 #include "mpfr-impl.h"
  24
  25 int
  26 mpfr_sqrt (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
  27 {
  28   mp_size_t rsize; /* number of limbs of r (plus 1 if exact limb multiple) */
  29   mp_size_t rrsize;
  30   mp_size_t usize; /* number of limbs of u */
  31   mp_size_t tsize; /* number of limbs of the sqrtrem remainder */
  32   mp_size_t k;
  33   mp_size_t l;
  34   mpfr_limb_ptr rp, rp0;
  35   mpfr_limb_ptr up;
  36   mpfr_limb_ptr sp;
  37   mp_limb_t sticky0; /* truncated part of input */
  38   mp_limb_t sticky1; /* truncated part of rp[0] */
  39   mp_limb_t sticky;
  40   int odd_exp;
  41   int sh; /* number of extra bits in rp[0] */
  42   int inexact; /* return ternary flag */
  43   mpfr_exp_t expr;
  44   MPFR_TMP_DECL(marker);
  45
  46   MPFR_LOG_FUNC
  47     (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (u), mpfr_log_prec, u, rnd_mode),
  48      ("y[%Pu]=%.*Rg inexact=%d",
  49       mpfr_get_prec (r), mpfr_log_prec, r, inexact));
  50
  51   if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(u)))
  52     {
  53       if (MPFR_IS_NAN(u))
  54         {
  55           MPFR_SET_NAN(r);
  56           MPFR_RET_NAN;
  57         }
  58       else if (MPFR_IS_ZERO(u))
  59         {
  60           /* 0+ or 0- */
  61           MPFR_SET_SAME_SIGN(r, u);
  62           MPFR_SET_ZERO(r);
  63           MPFR_RET(0); /* zero is exact */
  64         }
  65       else
  66         {
  67           MPFR_ASSERTD(MPFR_IS_INF(u));
  68           /* sqrt(-Inf) = NAN */
  69           if (MPFR_IS_NEG(u))
  70             {
  71               MPFR_SET_NAN(r);
  72               MPFR_RET_NAN;
  73             }
  74           MPFR_SET_POS(r);
  75           MPFR_SET_INF(r);
  76           MPFR_RET(0);
  77         }
  78     }
  79   if (MPFR_UNLIKELY(MPFR_IS_NEG(u)))
  80     {
  81       MPFR_SET_NAN(r);
  82       MPFR_RET_NAN;
  83     }
  84   MPFR_SET_POS(r);
  85
  86   MPFR_TMP_MARK (marker);
  87   MPFR_UNSIGNED_MINUS_MODULO(sh,MPFR_PREC(r));
  88   if (sh == 0 && rnd_mode == MPFR_RNDN)
  89     sh = GMP_NUMB_BITS; /* ugly case */
  90   rsize = MPFR_LIMB_SIZE(r) + (sh == GMP_NUMB_BITS);
  91   /* rsize is the number of limbs of r + 1 if exact limb multiple and rounding
  92      to nearest, this is the number of wanted limbs for the square root */
  93   rrsize = rsize + rsize;
  94   usize = MPFR_LIMB_SIZE(u); /* number of limbs of u */
  95   rp0 = MPFR_MANT(r);
  96   rp = (sh < GMP_NUMB_BITS) ? rp0 : MPFR_TMP_LIMBS_ALLOC (rsize);
  97   up = MPFR_MANT(u);
  98   sticky0 = MPFR_LIMB_ZERO; /* truncated part of input */
  99   sticky1 = MPFR_LIMB_ZERO; /* truncated part of rp[0] */
 100   odd_exp = (unsigned int) MPFR_GET_EXP (u) & 1;
 101   inexact = -1; /* return ternary flag */
 102
 103   sp = MPFR_TMP_LIMBS_ALLOC (rrsize);
 104
 105   /* copy the most significant limbs of u to {sp, rrsize} */
 106   if (MPFR_LIKELY(usize <= rrsize)) /* in case r and u have the same precision,
 107                                        we have indeed rrsize = 2 * usize */
 108     {
 109       k = rrsize - usize;
 110       if (MPFR_LIKELY(k))
 111         MPN_ZERO (sp, k);
 112       if (odd_exp)
 113         {
 114           if (MPFR_LIKELY(k))
 115             sp[k - 1] = mpn_rshift (sp + k, up, usize, 1);
 116           else
 117             sticky0 = mpn_rshift (sp, up, usize, 1);
 118         }
 119       else
 120         MPN_COPY (sp + rrsize - usize, up, usize);
 121     }
 122   else /* usize > rrsize: truncate the input */
 123     {
 124       k = usize - rrsize;
 125       if (odd_exp)
 126         sticky0 = mpn_rshift (sp, up + k, rrsize, 1);
 127       else
 128         MPN_COPY (sp, up + k, rrsize);
 129       l = k;
 130       while (sticky0 == MPFR_LIMB_ZERO && l != 0)
 131         sticky0 = up[--l];
 132     }
 133
 134   /* sticky0 is non-zero iff the truncated part of the input is non-zero */
 135
 136   /* mpn_rootrem with NULL 2nd argument is faster than mpn_sqrtrem, thus use
 137      it if available and if the user asked to use GMP internal functions */
 138 #if defined(WANT_GMP_INTERNALS) && defined(HAVE___GMPN_ROOTREM)
 139   tsize = __gmpn_rootrem (rp, NULL, sp, rrsize, 2);
 140 #else
 141   tsize = mpn_sqrtrem (rp, NULL, sp, rrsize);
 142 #endif
 143
 144   /* a return value of zero in mpn_sqrtrem indicates a perfect square */
 145   sticky = sticky0 || tsize != 0;
 146
 147   /* truncate low bits of rp[0] */
 148   sticky1 = rp[0] & ((sh < GMP_NUMB_BITS) ? MPFR_LIMB_MASK(sh)
 149                      : ~MPFR_LIMB_ZERO);
 150   rp[0] -= sticky1;
 151
 152   sticky = sticky || sticky1;
 153
 154   expr = (MPFR_GET_EXP(u) + odd_exp) / 2;  /* exact */
 155
 156   if (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDD || sticky == MPFR_LIMB_ZERO)
 157     {
 158       inexact = (sticky == MPFR_LIMB_ZERO) ? 0 : -1;
 159       goto truncate;
 160     }
 161   else if (rnd_mode == MPFR_RNDN)
 162     {
 163       /* if sh < GMP_NUMB_BITS, the round bit is bit (sh-1) of sticky1
 164                   and the sticky bit is formed by the low sh-1 bits from
 165                   sticky1, together with the sqrtrem remainder and sticky0. */
 166       if (sh < GMP_NUMB_BITS)
 167         {
 168           if (sticky1 & (MPFR_LIMB_ONE << (sh - 1)))
 169             { /* round bit is set */
 170               if (sticky1 == (MPFR_LIMB_ONE << (sh - 1)) && tsize == 0
 171                   && sticky0 == 0)
 172                 goto even_rule;
 173               else
 174                 goto add_one_ulp;
 175             }
 176           else /* round bit is zero */
 177             goto truncate; /* with the default inexact=-1 */
 178         }
 179       else /* sh = GMP_NUMB_BITS: the round bit is the most significant bit
 180               of rp[0], and the remaining GMP_NUMB_BITS-1 bits contribute to
 181               the sticky bit */
 182         {
 183           if (sticky1 & MPFR_LIMB_HIGHBIT)
 184             { /* round bit is set */
 185               if (sticky1 == MPFR_LIMB_HIGHBIT && tsize == 0 && sticky0 == 0)
 186                 goto even_rule;
 187               else
 188                 goto add_one_ulp;
 189             }
 190           else /* round bit is zero */
 191             goto truncate; /* with the default inexact=-1 */
 192         }
 193     }
 194   else /* rnd_mode=GMP_RDNU, necessarily sticky <> 0, thus add 1 ulp */
 195     goto add_one_ulp;
 196
 197  even_rule: /* has to set inexact */
 198   if (sh < GMP_NUMB_BITS)
 199     inexact = (rp[0] & (MPFR_LIMB_ONE << sh)) ? 1 : -1;
 200   else
 201     inexact = (rp[1] & MPFR_LIMB_ONE) ? 1 : -1;
 202   if (inexact == -1)
 203     goto truncate;
 204   /* else go through add_one_ulp */
 205
 206  add_one_ulp:
 207   inexact = 1; /* always here */
 208   if (sh == GMP_NUMB_BITS)
 209     {
 210       rp ++;
 211       rsize --;
 212       sh = 0;
 213     }
 214   if (mpn_add_1 (rp0, rp, rsize, MPFR_LIMB_ONE << sh))
 215     {
 216       expr ++;
 217       rp[rsize - 1] = MPFR_LIMB_HIGHBIT;
 218     }
 219   goto end;
 220
 221  truncate: /* inexact = 0 or -1 */
 222   if (sh == GMP_NUMB_BITS)
 223     MPN_COPY (rp0, rp + 1, rsize - 1);
 224
 225  end:
 226   MPFR_ASSERTN (expr >= MPFR_EMIN_MIN && expr <= MPFR_EMAX_MAX);
 227   MPFR_EXP (r) = expr;
 228   MPFR_TMP_FREE(marker);
 229
 230   return mpfr_check_range (r, inexact, rnd_mode);
 231 }