Annotation of OpenXM_contrib/gmp/mpfr/set_d.c, Revision 1.1
1.1 ! maekawa 1: /* mpfr_set_d, mpfr_get_d -- convert a multiple precision floating-point number
! 2: from/to a machine double precision float
! 3:
! 4: Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
! 5:
! 6: This file is part of the MPFR Library.
! 7:
! 8: The MPFR Library is free software; you can redistribute it and/or modify
! 9: it under the terms of the GNU Library General Public License as published by
! 10: the Free Software Foundation; either version 2 of the License, or (at your
! 11: option) any later version.
! 12:
! 13: The 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 Library General Public
! 16: License for more details.
! 17:
! 18: You should have received a copy of the GNU Library General Public License
! 19: along with the MPFR Library; see the file COPYING.LIB. If not, write to
! 20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 21: MA 02111-1307, USA. */
! 22:
! 23: #if __GNUC__ /* gcc "patched" headers seem to omit isnan... */
! 24: extern int isnan(double);
! 25: #endif
! 26: #include <math.h> /* for isnan and NaN */
! 27:
! 28: #include "gmp.h"
! 29: #include "gmp-impl.h"
! 30: #include "longlong.h"
! 31: #include "mpfr.h"
! 32:
! 33: #define NaN sqrt(-1) /* ensures a machine-independent NaN */
! 34:
! 35: /* Included from gmp-2.0.2, patched to support denorms */
! 36:
! 37: #ifdef XDEBUG
! 38: #undef _GMP_IEEE_FLOATS
! 39: #endif
! 40:
! 41: #ifndef _GMP_IEEE_FLOATS
! 42: #define _GMP_IEEE_FLOATS 0
! 43: #endif
! 44:
! 45: int
! 46: #if __STDC__
! 47: __mpfr_extract_double (mp_ptr rp, double d, int e)
! 48: #else
! 49: __mpfr_extract_double (rp, d, e)
! 50: mp_ptr rp;
! 51: double d;
! 52: int e;
! 53: #endif
! 54: /* e=0 iff BITS_PER_MP_LIMB=32 and rp has only one limb */
! 55: {
! 56: long exp;
! 57: mp_limb_t manl;
! 58: #if BITS_PER_MP_LIMB == 32
! 59: mp_limb_t manh;
! 60: #endif
! 61:
! 62: /* BUGS
! 63:
! 64: 1. Should handle Inf and NaN in IEEE specific code.
! 65: 2. Handle Inf and NaN also in default code, to avoid hangs.
! 66: 3. Generalize to handle all BITS_PER_MP_LIMB >= 32.
! 67: 4. This lits is incomplete and misspelled.
! 68: */
! 69:
! 70: if (d == 0.0)
! 71: {
! 72: rp[0] = 0;
! 73: #if BITS_PER_MP_LIMB == 32
! 74: if (e) rp[1] = 0;
! 75: #endif
! 76: return 0;
! 77: }
! 78:
! 79: #if _GMP_IEEE_FLOATS
! 80: {
! 81: union ieee_double_extract x;
! 82: x.d = d;
! 83:
! 84: exp = x.s.exp;
! 85: if (exp)
! 86: {
! 87: #if BITS_PER_MP_LIMB == 64
! 88: manl = (((mp_limb_t) 1 << 63)
! 89: | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11));
! 90: #else
! 91: manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
! 92: manl = x.s.manl << 11;
! 93: #endif
! 94: }
! 95: else
! 96: {
! 97: #if BITS_PER_MP_LIMB == 64
! 98: manl = ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11);
! 99: #else
! 100: manh = (x.s.manh << 11) | (x.s.manl >> 21);
! 101: manl = x.s.manl << 11;
! 102: #endif
! 103: }
! 104: }
! 105: #else
! 106: {
! 107: /* Unknown (or known to be non-IEEE) double format. */
! 108: exp = 0;
! 109: if (d >= 1.0)
! 110: {
! 111: if (d * 0.5 == d)
! 112: abort ();
! 113:
! 114: while (d >= 32768.0)
! 115: {
! 116: d *= (1.0 / 65536.0);
! 117: exp += 16;
! 118: }
! 119: while (d >= 1.0)
! 120: {
! 121: d *= 0.5;
! 122: exp += 1;
! 123: }
! 124: }
! 125: else if (d < 0.5)
! 126: {
! 127: while (d < (1.0 / 65536.0))
! 128: {
! 129: d *= 65536.0;
! 130: exp -= 16;
! 131: }
! 132: while (d < 0.5)
! 133: {
! 134: d *= 2.0;
! 135: exp -= 1;
! 136: }
! 137: }
! 138:
! 139: d *= MP_BASE_AS_DOUBLE;
! 140: #if BITS_PER_MP_LIMB == 64
! 141: manl = d;
! 142: #else
! 143: manh = d;
! 144: manl = (d - manh) * MP_BASE_AS_DOUBLE;
! 145: #endif
! 146:
! 147: exp += 1022;
! 148: }
! 149: #endif
! 150:
! 151: if (exp) exp = (unsigned) exp - 1022; else exp = -1021;
! 152:
! 153: #if BITS_PER_MP_LIMB == 64
! 154: rp[0] = manl;
! 155: #else
! 156: if (e) {
! 157: rp[1] = manh;
! 158: rp[0] = manl;
! 159: }
! 160: else {
! 161: rp[0] = manh;
! 162: }
! 163: #endif
! 164:
! 165: return exp;
! 166: }
! 167:
! 168: /* End of part included from gmp-2.0.2 */
! 169: /* Part included from gmp temporary releases */
! 170: double
! 171: #if __STDC__
! 172: __mpfr_scale2 (double d, int exp)
! 173: #else
! 174: __mpfr_scale2 (d, exp)
! 175: double d;
! 176: int exp;
! 177: #endif
! 178: {
! 179: #if _GMP_IEEE_FLOATS
! 180: {
! 181: union ieee_double_extract x;
! 182: x.d = d;
! 183: exp += x.s.exp;
! 184: x.s.exp = exp;
! 185: if (exp >= 2047)
! 186: {
! 187: /* Return +-infinity */
! 188: x.s.exp = 2047;
! 189: x.s.manl = x.s.manh = 0;
! 190: }
! 191: else if (exp < 1)
! 192: {
! 193: x.s.exp = 1; /* smallest exponent (biased) */
! 194: /* Divide result by 2 until we have scaled it to the right IEEE
! 195: denormalized number, but stop if it becomes zero. */
! 196: while (exp < 1 && x.d != 0)
! 197: {
! 198: x.d *= 0.5;
! 199: exp++;
! 200: }
! 201: }
! 202: return x.d;
! 203: }
! 204: #else
! 205: {
! 206: double factor, r;
! 207:
! 208: factor = 2.0;
! 209: if (exp < 0)
! 210: {
! 211: factor = 0.5;
! 212: exp = -exp;
! 213: }
! 214: r = d;
! 215: if (exp != 0)
! 216: {
! 217: if ((exp & 1) != 0)
! 218: r *= factor;
! 219: exp >>= 1;
! 220: while (exp != 0)
! 221: {
! 222: factor *= factor;
! 223: if ((exp & 1) != 0)
! 224: r *= factor;
! 225: exp >>= 1;
! 226: }
! 227: }
! 228: return r;
! 229: }
! 230: #endif
! 231: }
! 232:
! 233:
! 234: /* End of part included from gmp */
! 235:
! 236: void
! 237: #if __STDC__
! 238: mpfr_set_d(mpfr_t r, double d, unsigned char rnd_mode)
! 239: #else
! 240: mpfr_set_d(r, d, rnd_mode)
! 241: mpfr_t r;
! 242: double d;
! 243: unsigned char rnd_mode;
! 244: #endif
! 245: {
! 246: int signd, sizer; unsigned int cnt;
! 247:
! 248: if (d == 0) { SET_ZERO(r); return; }
! 249: else if (isnan(d)) { SET_NAN(r); return; }
! 250:
! 251: signd = (d < 0) ? -1 : 1;
! 252: d = ABS (d);
! 253: sizer = (PREC(r)-1)/BITS_PER_MP_LIMB + 1;
! 254:
! 255: /* warning: __mpfr_extract_double requires at least two limbs */
! 256: if (sizer < MPFR_LIMBS_PER_DOUBLE)
! 257: EXP(r) = __mpfr_extract_double (MANT(r), d, 0);
! 258: else
! 259: EXP(r) = __mpfr_extract_double (MANT(r) + sizer - MPFR_LIMBS_PER_DOUBLE, d, 1);
! 260:
! 261: if (sizer > MPFR_LIMBS_PER_DOUBLE)
! 262: MPN_ZERO(MANT(r), sizer - MPFR_LIMBS_PER_DOUBLE);
! 263:
! 264: count_leading_zeros(cnt, MANT(r)[sizer-1]);
! 265: if (cnt) mpn_lshift(MANT(r), MANT(r), sizer, cnt);
! 266:
! 267: EXP(r) -= cnt;
! 268: if (SIZE(r)*signd<0) CHANGE_SIGN(r);
! 269:
! 270: mpfr_round(r, rnd_mode, PREC(r));
! 271: return;
! 272: }
! 273:
! 274: double
! 275: #if __STDC__
! 276: mpfr_get_d2(mpfr_srcptr src, long e)
! 277: #else
! 278: mpfr_get_d2(src, e)
! 279: mpfr_srcptr(src);
! 280: long e;
! 281: #endif
! 282: {
! 283: double res;
! 284: mp_size_t size, i, n_limbs_to_use;
! 285: mp_ptr qp;
! 286: int negative;
! 287:
! 288: if (FLAG_NAN(src)) {
! 289: #ifdef DEBUG
! 290: printf("recognized NaN\n");
! 291: #endif
! 292: return NaN; }
! 293: if (NOTZERO(src)==0) return 0.0;
! 294: size = 1+(PREC(src)-1)/BITS_PER_MP_LIMB;
! 295: qp = MANT(src);
! 296: negative = (SIGN(src)==-1);
! 297:
! 298: /* Warning: don't compute the abs(res) and set the sign afterwards,
! 299: otherwise the current machine rounding mode will not be taken
! 300: correctly into account. */
! 301: /* res = (negative) ? -(double)qp[size - 1] : qp[size - 1]; */
! 302: res = 0.0;
! 303: /* Warning: an arbitrary number of limbs may be required for an exact
! 304: rounding. The following code is correct but not optimal since one
! 305: may be able to decide without considering all limbs. */
! 306: /* n_limbs_to_use = MIN (MPFR_LIMBS_PER_DOUBLE, size); */
! 307: n_limbs_to_use = size;
! 308: /* Accumulate the limbs from less significant to most significant
! 309: otherwise due to rounding we may accumulate several ulps,
! 310: especially in rounding towards -/+infinity. */
! 311: for (i = n_limbs_to_use; i>=1; i--)
! 312: res = res / MP_BASE_AS_DOUBLE +
! 313: ((negative) ? -(double)qp[size - i] : qp[size - i]);
! 314: res = __mpfr_scale2 (res, e - BITS_PER_MP_LIMB);
! 315:
! 316: return res;
! 317: }
! 318:
! 319: double
! 320: #if __STDC__
! 321: mpfr_get_d(mpfr_srcptr src)
! 322: #else
! 323: mpfr_get_d(src)
! 324: mpfr_srcptr src;
! 325: #endif
! 326: {
! 327: return mpfr_get_d2(src, EXP(src));
! 328: }
! 329:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to set_d.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr