contrib/mpfr/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 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 #include "mpfr-impl.h"
  24
  25 int
  26 mpfr_sqrt (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode)
  27 {
  28   mp_size_t rsize; /* number of limbs of r */
  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   mp_ptr rp;
  35   mp_ptr up;
  36   mp_ptr sp;
  37   mp_ptr tp;
  38   mp_limb_t sticky0; /* truncated part of input */
  39   mp_limb_t sticky1; /* truncated part of rp[0] */
  40   mp_limb_t sticky;
  41   int odd_exp;
  42   int sh; /* number of extra bits in rp[0] */
  43   int inexact; /* return ternary flag */
  44   mp_exp_t expr;
  45   MPFR_TMP_DECL(marker);
  46
  47   MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", u, u, rnd_mode),
  48                  ("y[%#R]=%R inexact=%d", r, r, inexact));
  49
  50   if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(u)))
  51     {
  52       if (MPFR_IS_NAN(u))
  53         {
  54           MPFR_SET_NAN(r);
  55           MPFR_RET_NAN;
  56         }
  57       else if (MPFR_IS_ZERO(u))
  58         {
  59           /* 0+ or 0- */
  60           MPFR_SET_SAME_SIGN(r, u);
  61           MPFR_SET_ZERO(r);
  62           MPFR_RET(0); /* zero is exact */
  63         }
  64       else
  65         {
  66           MPFR_ASSERTD(MPFR_IS_INF(u));
  67           /* sqrt(-Inf) = NAN */
  68           if (MPFR_IS_NEG(u))
  69             {
  70               MPFR_SET_NAN(r);
  71               MPFR_RET_NAN;
  72             }
  73           MPFR_SET_POS(r);
  74           MPFR_SET_INF(r);
  75           MPFR_RET(0);
  76         }
  77     }
  78   if (MPFR_UNLIKELY(MPFR_IS_NEG(u)))
  79     {
  80       MPFR_SET_NAN(r);
  81       MPFR_RET_NAN;
  82     }
  83   MPFR_CLEAR_FLAGS(r);
  84   MPFR_SET_POS(r);
  85
  86   rsize = MPFR_LIMB_SIZE(r); /* number of limbs of r */
  87   rrsize = rsize + rsize;
  88   usize = MPFR_LIMB_SIZE(u); /* number of limbs of u */
  89   rp = MPFR_MANT(r);
  90   up = MPFR_MANT(u);
  91   sticky0 = MPFR_LIMB_ZERO; /* truncated part of input */
  92   sticky1 = MPFR_LIMB_ZERO; /* truncated part of rp[0] */
  93   odd_exp = (unsigned int) MPFR_GET_EXP (u) & 1;
  94   inexact = -1; /* return ternary flag */
  95
  96   MPFR_TMP_MARK (marker);
  97   sp = (mp_limb_t *) MPFR_TMP_ALLOC (rrsize * sizeof (mp_limb_t));
  98
  99   /* copy the most significant limbs of u to {sp, rrsize} */
 100   if (MPFR_LIKELY(usize <= rrsize)) /* in case r and u have the same precision,
 101                                        we have indeed rrsize = 2 * usize */
 102     {
 103       k = rrsize - usize;
 104       if (MPFR_LIKELY(k))
 105         MPN_ZERO (sp, k);
 106       if (odd_exp)
 107         {
 108           if (MPFR_LIKELY(k))
 109             sp[k - 1] = mpn_rshift (sp + k, up, usize, 1);
 110           else
 111             sticky0 = mpn_rshift (sp, up, usize, 1);
 112         }
 113       else
 114         MPN_COPY (sp + rrsize - usize, up, usize);
 115     }
 116   else /* usize > rrsize: truncate the input */
 117     {
 118       k = usize - rrsize;
 119       if (odd_exp)
 120         sticky0 = mpn_rshift (sp, up + k, rrsize, 1);
 121       else
 122         MPN_COPY (sp, up + k, rrsize);
 123       l = k;
 124       while (sticky0 == MPFR_LIMB_ZERO && l != 0)
 125         sticky0 = up[--l];
 126     }
 127
 128   /* sticky0 is non-zero iff the truncated part of the input is non-zero */
 129
 130   tsize = mpn_sqrtrem (rp, tp = sp, sp, rrsize);
 131
 132   l = tsize;
 133   sticky = sticky0;
 134   while (sticky == MPFR_LIMB_ZERO && l != 0)
 135     sticky = tp[--l];
 136
 137   /* truncated low bits of rp[0] */
 138   MPFR_UNSIGNED_MINUS_MODULO(sh,MPFR_PREC(r));
 139   sticky1 = rp[0] & MPFR_LIMB_MASK(sh);
 140   rp[0] -= sticky1;
 141
 142   sticky = sticky || sticky1;
 143
 144   expr = (MPFR_GET_EXP(u) + odd_exp) / 2;  /* exact */
 145
 146   if (rnd_mode == GMP_RNDZ || rnd_mode == GMP_RNDD || sticky == MPFR_LIMB_ZERO)
 147     {
 148       inexact = (sticky == MPFR_LIMB_ZERO) ? 0 : -1;
 149       goto truncate;
 150     }
 151   else if (rnd_mode == GMP_RNDN)
 152     {
 153       /* if sh>0, the round bit is bit (sh-1) of sticky1
 154                   and the sticky bit is formed by the low sh-1 bits from
 155                   sticky1, together with {tp, tsize} and sticky0. */
 156       if (sh)
 157         {
 158           if (sticky1 & (MPFR_LIMB_ONE << (sh - 1)))
 159             { /* round bit is set */
 160               if (sticky1 == (MPFR_LIMB_ONE << (sh - 1)) && tsize == 0
 161                   && sticky0 == 0)
 162                 goto even_rule;
 163               else
 164                 goto add_one_ulp;
 165             }
 166           else /* round bit is zero */
 167             goto truncate; /* with the default inexact=-1 */
 168         }
 169       else
 170         {
 171           /* if sh=0, we have to compare {tp, tsize} with {rp, rsize}:
 172                 if {tp, tsize} < {rp, rsize}: truncate
 173                 if {tp, tsize} > {rp, rsize}: round up
 174                 if {tp, tsize} = {rp, rsize}: compare the truncated part of the
 175                                               input to 1/4
 176                    if < 1/4: truncate
 177                    if > 1/4: round up
 178                    if = 1/4: even rounding rule
 179              Set inexact = -1 if truncate
 180                  inexact = 1 if add one ulp
 181                  inexact = 0 if even rounding rule
 182           */
 183           if (tsize < rsize)
 184             inexact = -1;
 185           else if (tsize > rsize) /* FIXME: may happen? */
 186             inexact = 1;
 187           else /* tsize = rsize */
 188             {
 189               int cmp;
 190
 191               cmp = mpn_cmp (tp, rp, rsize);
 192               if (cmp > 0)
 193                 inexact = 1;
 194               else if (cmp < 0 || sticky0 == MPFR_LIMB_ZERO)
 195                 inexact = -1;
 196               /* now tricky case {tp, tsize} = {rp, rsize} */
 197               /* in case usize <= rrsize, the only case where sticky0 <> 0
 198                  is when the exponent of u is odd and usize = rrsize (k=0),
 199                  but in that case the truncated part is exactly 1/2, thus
 200                  we have to round up.
 201                  If the exponent of u is odd, and up[k] is odd, the truncated
 202                  part is >= 1/2, so we round up too. */
 203               else if (usize <= rrsize || (odd_exp && (up[k] & MPFR_LIMB_ONE)))
 204                 inexact = 1;
 205               else
 206                 {
 207                   /* now usize > rrsize:
 208                      (a) if the exponent of u is even, the 1/4 bit is the
 209                      2nd most significant bit of up[k-1];
 210                      (b) if the exponent of u is odd, the 1/4 bit is the
 211                      1st most significant bit of up[k-1]; */
 212                   sticky1 = MPFR_LIMB_ONE << (BITS_PER_MP_LIMB - 2 + odd_exp);
 213                   if (up[k - 1] < sticky1)
 214                     inexact = -1;
 215                   else if (up[k - 1] > sticky1)
 216                     inexact = 1;
 217                   else
 218                     {
 219                       /* up[k - 1] == sticky1: consider low k-1 limbs */
 220                       while (--k > 0 && up[k - 1] == MPFR_LIMB_ZERO)
 221                         ;
 222                       inexact = (k != 0);
 223                     }
 224                 } /* end of case {tp, tsize} = {rp, rsize} */
 225             } /* end of case tsize = rsize */
 226           if (inexact == -1)
 227             goto truncate;
 228           else if (inexact == 1)
 229             goto add_one_ulp;
 230           /* else go through even_rule */
 231         }
 232     }
 233   else /* rnd_mode=GMP_RDNU, necessarily sticky <> 0, thus add 1 ulp */
 234     goto add_one_ulp;
 235
 236  even_rule: /* has to set inexact */
 237   inexact = (rp[0] & (MPFR_LIMB_ONE << sh)) ? 1 : -1;
 238   if (inexact == -1)
 239     goto truncate;
 240   /* else go through add_one_ulp */
 241
 242  add_one_ulp:
 243   inexact = 1; /* always here */
 244   if (mpn_add_1 (rp, rp, rsize, MPFR_LIMB_ONE << sh))
 245     {
 246       expr ++;
 247       rp[rsize - 1] = MPFR_LIMB_HIGHBIT;
 248     }
 249
 250  truncate: /* inexact = 0 or -1 */
 251
 252   MPFR_ASSERTN (expr >= MPFR_EMIN_MIN && expr <= MPFR_EMAX_MAX);
 253   MPFR_EXP (r) = expr;
 254
 255   MPFR_TMP_FREE(marker);
 256   return mpfr_check_range (r, inexact, rnd_mode);
 257 }