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 } 232