Annotation of OpenXM_contrib/gmp/mpfr/sqrt.c, Revision 1.1.1.2
1.1 maekawa 1: /* mpfr_sqrt -- square root of a floating-point number
2:
1.1.1.2 ! ohara 3: Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
1.1 maekawa 4:
5: This file is part of the MPFR Library.
6:
7: The MPFR Library is free software; you can redistribute it and/or modify
1.1.1.2 ! ohara 8: it under the terms of the GNU Lesser General Public License as published by
! 9: the Free Software Foundation; either version 2.1 of the License, or (at your
1.1 maekawa 10: option) any later version.
11:
12: The MPFR Library is distributed in the hope that it will be useful, but
13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.1.1.2 ! ohara 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 15: License for more details.
16:
1.1.1.2 ! ohara 17: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 18: along with the MPFR Library; see the file COPYING.LIB. If not, write to
19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20: MA 02111-1307, USA. */
21:
22: #include "gmp.h"
23: #include "gmp-impl.h"
24: #include "mpfr.h"
1.1.1.2 ! ohara 25: #include "mpfr-impl.h"
1.1 maekawa 26:
27: /* #define DEBUG */
28:
29: int
1.1.1.2 ! ohara 30: mpfr_sqrt (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode)
1.1 maekawa 31: {
32: mp_ptr up, rp, tmp, remp;
33: mp_size_t usize, rrsize;
34: mp_size_t rsize;
1.1.1.2 ! ohara 35: mp_size_t err;
1.1 maekawa 36: mp_limb_t q_limb;
1.1.1.2 ! ohara 37: int odd_exp_u;
1.1 maekawa 38: long rw, nw, k;
1.1.1.2 ! ohara 39: int inexact = 0, t;
1.1 maekawa 40: unsigned long cc = 0;
1.1.1.2 ! ohara 41: int can_round = 0;
1.1 maekawa 42:
1.1.1.2 ! ohara 43: TMP_DECL(marker);
1.1 maekawa 44:
1.1.1.2 ! ohara 45: if (MPFR_IS_NAN(u))
! 46: {
! 47: MPFR_SET_NAN(r);
! 48: MPFR_RET_NAN;
! 49: }
! 50:
! 51: if (MPFR_SIGN(u) < 0)
! 52: {
! 53: if (MPFR_IS_INF(u) || MPFR_NOTZERO(u))
! 54: {
! 55: MPFR_SET_NAN(r);
! 56: MPFR_RET_NAN;
! 57: }
! 58: else
! 59: { /* sqrt(-0) = -0 */
! 60: MPFR_CLEAR_FLAGS(r);
! 61: MPFR_SET_ZERO(r);
! 62: MPFR_SET_NEG(r);
! 63: MPFR_RET(0);
! 64: }
! 65: }
! 66:
! 67: MPFR_CLEAR_NAN(r);
! 68: MPFR_SET_POS(r);
! 69:
! 70: if (MPFR_IS_INF(u))
! 71: {
! 72: MPFR_SET_INF(r);
! 73: MPFR_RET(0);
! 74: }
! 75:
! 76: MPFR_CLEAR_INF(r);
! 77:
! 78: if (MPFR_IS_ZERO(u))
! 79: {
! 80: MPFR_SET_ZERO(r);
! 81: MPFR_RET(0); /* zero is exact */
! 82: }
! 83:
! 84: up = MPFR_MANT(u);
! 85: usize = (MPFR_PREC(u) - 1)/BITS_PER_MP_LIMB + 1;
1.1 maekawa 86:
87: #ifdef DEBUG
1.1.1.2 ! ohara 88: printf("Entering square root : ");
! 89: for(k = usize - 1; k >= 0; k--) { printf("%lu ", up[k]); }
! 90: printf(".\n");
1.1 maekawa 91: #endif
92:
1.1.1.2 ! ohara 93: /* Compare the mantissas */
1.1 maekawa 94:
1.1.1.2 ! ohara 95: odd_exp_u = (unsigned int) MPFR_EXP(u) & 1;
! 96: MPFR_ASSERTN(MPFR_PREC(r) <= MPFR_INTPREC_MAX - 3);
! 97: rrsize = (MPFR_PREC(r) + 2 + odd_exp_u) / BITS_PER_MP_LIMB + 1;
! 98: MPFR_ASSERTN(rrsize <= MP_SIZE_T_MAX/2);
! 99: rsize = rrsize << 1;
! 100: /* One extra bit is needed in order to get the square root with enough
! 101: precision ; take one extra bit for rrsize in order to solve more
! 102: easily the problem of rounding to nearest.
! 103: Need to have 2*rrsize = rsize...
! 104: Take one extra bit if the exponent of u is odd since we shall have
! 105: to shift then.
! 106: */
! 107:
! 108: TMP_MARK(marker);
! 109: if (odd_exp_u) /* Shift u one bit to the right */
! 110: {
! 111: if (MPFR_PREC(u) & (BITS_PER_MP_LIMB - 1))
! 112: {
! 113: up = TMP_ALLOC(usize * BYTES_PER_MP_LIMB);
! 114: mpn_rshift(up, MPFR_MANT(u), usize, 1);
! 115: }
! 116: else
! 117: {
! 118: up = TMP_ALLOC((usize + 1) * BYTES_PER_MP_LIMB);
! 119: if (mpn_rshift(up + 1, MPFR_MANT(u), usize, 1))
! 120: up[0] = GMP_LIMB_HIGHBIT;
! 121: else
! 122: up[0] = 0;
! 123: usize++;
! 124: }
1.1 maekawa 125: }
126:
1.1.1.2 ! ohara 127: MPFR_EXP(r) = MPFR_EXP(u) != MPFR_EMAX_MAX
! 128: ? (MPFR_EXP(u) + odd_exp_u) / 2
! 129: : (MPFR_EMAX_MAX - 1) / 2 + 1;
! 130:
! 131: do
! 132: {
! 133:
! 134: err = rsize * BITS_PER_MP_LIMB;
! 135: if (rsize < usize)
! 136: err--;
! 137: if (err > rrsize * BITS_PER_MP_LIMB)
! 138: err = rrsize * BITS_PER_MP_LIMB;
! 139:
! 140: tmp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
! 141: rp = (mp_ptr) TMP_ALLOC (rrsize * BYTES_PER_MP_LIMB);
! 142: remp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
! 143:
! 144: if (usize >= rsize)
! 145: {
! 146: MPN_COPY (tmp, up + usize - rsize, rsize);
! 147: }
! 148: else
! 149: {
! 150: MPN_COPY (tmp + rsize - usize, up, usize);
! 151: MPN_ZERO (tmp, rsize - usize);
! 152: }
! 153:
! 154: /* Do the real job */
1.1 maekawa 155:
156: #ifdef DEBUG
1.1.1.2 ! ohara 157: printf("Taking the sqrt of : ");
! 158: for(k = rsize - 1; k >= 0; k--)
! 159: printf("+%lu*2^%lu",tmp[k],k*BITS_PER_MP_LIMB);
! 160: printf(".\n");
1.1 maekawa 161: #endif
162:
1.1.1.2 ! ohara 163: q_limb = mpn_sqrtrem (rp, remp, tmp, rsize);
1.1 maekawa 164:
165: #ifdef DEBUG
1.1.1.2 ! ohara 166: printf ("The result is : \n");
! 167: printf ("sqrt : ");
! 168: for (k = rrsize - 1; k >= 0; k--)
! 169: printf ("%lu ", rp[k]);
! 170: printf ("(inexact = %lu)\n", q_limb);
1.1 maekawa 171: #endif
172:
1.1.1.2 ! ohara 173: can_round = mpfr_can_round_raw(rp, rrsize, 1, err,
! 174: GMP_RNDZ, rnd_mode, MPFR_PREC(r));
1.1 maekawa 175:
1.1.1.2 ! ohara 176: /* If we used all the limbs of both the dividend and the divisor,
! 177: then we have the correct RNDZ rounding */
! 178:
! 179: if (!can_round && (rsize < 2*usize))
! 180: {
1.1 maekawa 181: #ifdef DEBUG
1.1.1.2 ! ohara 182: printf("Increasing the precision.\n");
! 183: #endif
! 184: }
! 185: }
! 186: while (!can_round && (rsize < 2*usize) && (rsize += 2) && (rrsize++));
! 187: #ifdef DEBUG
! 188: printf ("can_round = %d\n", can_round);
1.1 maekawa 189: #endif
190:
1.1.1.2 ! ohara 191: /* This part may be deplaced upper to avoid a few mpfr_can_round_raw */
! 192: /* when the square root is exact. It is however very unprobable that */
! 193: /* it would improve the behaviour of the present code on average. */
! 194:
! 195: if (!q_limb) /* the sqrtrem call was exact, possible exact square root */
! 196: {
! 197: /* if we have taken into account the whole of up */
! 198: for (k = usize - rsize - 1; k >= 0; k++)
! 199: if (up[k] != 0)
! 200: break;
! 201:
! 202: if (k < 0)
! 203: goto fin; /* exact square root ==> inexact = 0 */
! 204: }
! 205:
! 206: if (can_round)
! 207: {
! 208: cc = mpfr_round_raw (rp, rp, err, 0, MPFR_PREC(r), rnd_mode, &inexact);
! 209: if (inexact == 0) /* exact high part: inex flag depends on remainder */
! 210: inexact = -q_limb;
! 211: rrsize = (MPFR_PREC(r) - 1)/BITS_PER_MP_LIMB + 1;
! 212: }
! 213: else
! 214: {
! 215: /* Use the return value of sqrtrem to decide of the rounding */
! 216: /* Note that at this point the sqrt has been computed */
! 217: /* EXACTLY. */
! 218: switch (rnd_mode)
! 219: {
! 220: case GMP_RNDZ :
! 221: case GMP_RNDD :
! 222: inexact = -1; /* result is truncated */
! 223: break;
! 224:
! 225: case GMP_RNDN :
! 226: /* Not in the situation ...0 111111 */
! 227: rw = (MPFR_PREC(r) + 1) & (BITS_PER_MP_LIMB - 1);
! 228: if (rw != 0)
! 229: {
! 230: rw = BITS_PER_MP_LIMB - rw;
! 231: nw = 0;
! 232: }
! 233: else
! 234: nw = 1;
! 235: if ((rp[nw] >> rw) & 1 && /* Not 0111111111 */
! 236: (q_limb || /* Nonzero remainder */
! 237: (rw ? (rp[nw] >> (rw + 1)) & 1 :
! 238: (rp[nw] >> (BITS_PER_MP_LIMB - 1)) & 1))) /* or even r. */
! 239: {
! 240: cc = mpn_add_1 (rp + nw, rp + nw, rrsize, MP_LIMB_T_ONE << rw);
! 241: inexact = 1;
! 242: }
! 243: else
! 244: inexact = -1;
! 245: break;
1.1 maekawa 246:
1.1.1.2 ! ohara 247: case GMP_RNDU:
! 248: /* we should arrive here only when the result is inexact, i.e.
! 249: either q_limb > 0 (the remainder from mpn_sqrtrem is non-zero)
! 250: or up[0..usize-rsize-1] is non zero, thus we have to add one
! 251: ulp, and inexact = 1 */
! 252: inexact = 1;
! 253: t = MPFR_PREC(r) & (BITS_PER_MP_LIMB - 1);
! 254: rsize = (MPFR_PREC(r) - 1)/BITS_PER_MP_LIMB + 1;
! 255: cc = mpn_add_1 (rp + rrsize - rsize, rp + rrsize - rsize, rsize,
! 256: t != 0 ?
! 257: MP_LIMB_T_ONE << (BITS_PER_MP_LIMB - t) :
! 258: MP_LIMB_T_ONE);
! 259: }
1.1 maekawa 260: }
261:
1.1.1.2 ! ohara 262: if (cc)
! 263: {
! 264: /* Is a shift necessary here? Isn't the result 1.0000...? */
! 265: mpn_rshift (rp, rp, rrsize, 1);
! 266: rp[rrsize-1] |= GMP_LIMB_HIGHBIT;
! 267: MPFR_EXP(r)++;
! 268: }
1.1 maekawa 269:
270: fin:
1.1.1.2 ! ohara 271: rsize = rrsize;
! 272: rrsize = (MPFR_PREC(r) - 1)/BITS_PER_MP_LIMB + 1;
! 273: MPN_COPY(MPFR_MANT(r), rp + rsize - rrsize, rrsize);
! 274:
! 275: if (MPFR_PREC(r) & (BITS_PER_MP_LIMB - 1))
! 276: MPFR_MANT(r)[0] &= ~((MP_LIMB_T_ONE <<
! 277: (BITS_PER_MP_LIMB -
! 278: (MPFR_PREC(r) & (BITS_PER_MP_LIMB - 1)))) - 1);
! 279:
! 280: TMP_FREE(marker);
! 281: return inexact;
1.1 maekawa 282: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to sqrt.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr