Annotation of OpenXM_contrib/gmp/mpfr/get_str.c, Revision 1.1.1.2
1.1 maekawa 1: /* mpfr_get_str -- output a floating-point number to a string
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 <stdlib.h>
23: #include <string.h>
24: #include "gmp.h"
25: #include "gmp-impl.h"
26: #include "longlong.h"
27: #include "mpfr.h"
1.1.1.2 ! ohara 28: #include "mpfr-impl.h"
! 29:
! 30: #define ERR 5
1.1 maekawa 31:
32: /*
33: Convert op to a string in base 'base' with 'n' digits and writes the
34: mantissa in 'str', the exponent in 'expptr'.
35: The result is rounded wrt 'rnd_mode'.
36:
37: For op = 3.1416 we get str = "31416" and expptr=1.
38: */
1.1.1.2 ! ohara 39: char*
! 40: mpfr_get_str (char *str, mp_exp_t *expptr, int base, size_t n,
! 41: mpfr_srcptr op, mp_rnd_t rnd_mode)
1.1 maekawa 42: {
1.1.1.2 ! ohara 43: int neg;
! 44:
! 45: if (base < 2 || base > 36)
! 46: return NULL;
1.1 maekawa 47:
1.1.1.2 ! ohara 48: if (n == 0) /* determine n from precision of op */
! 49: {
! 50: n = __mp_bases[base].chars_per_bit_exactly * MPFR_PREC(op);
! 51: if (n < 2)
! 52: n = 2;
1.1 maekawa 53: }
54:
1.1.1.2 ! ohara 55: /* Do not use MPFR_PREC_MIN as this applies to base 2 only. Perhaps we
! 56: should allow n == 1 for directed rounding modes and odd bases. */
! 57: if (n < 2)
! 58: return NULL;
! 59:
! 60: if (MPFR_IS_NAN(op))
! 61: {
! 62: if (str == NULL)
! 63: str = (*__gmp_allocate_func)(4);
! 64: str[0] = 'N';
! 65: str[1] = 'a';
! 66: str[2] = 'N';
! 67: str[3] = '\0';
! 68: return str;
1.1 maekawa 69: }
70:
1.1.1.2 ! ohara 71: neg = MPFR_SIGN(op) < 0; /* 0 if positive, 1 if negative */
! 72:
! 73: if (MPFR_IS_INF(op))
! 74: {
! 75: char *str0;
! 76:
! 77: if (str == NULL)
! 78: str = (*__gmp_allocate_func)(neg + 4);
! 79: str0 = str;
! 80: if (neg)
! 81: *str++ = '-';
! 82: *str++ = 'I';
! 83: *str++ = 'n';
! 84: *str++ = 'f';
! 85: *str = '\0';
! 86: return str0; /* strlen(str0) = neg + 3 */
1.1 maekawa 87: }
88:
1.1.1.2 ! ohara 89: /* op is a floating-point number */
! 90:
! 91: if (MPFR_IS_ZERO(op))
! 92: {
! 93: char *str0;
! 94:
! 95: if (str == NULL)
! 96: str = (*__gmp_allocate_func)(neg + n + 1);
! 97: str0 = str;
! 98: if (neg)
! 99: *str++ = '-';
! 100: memset(str, '0', n);
! 101: str[n] = '\0';
! 102: *expptr = 0; /* a bit like frexp() in ISO C99 */
! 103: return str0; /* strlen(str0) = neg + n */
1.1 maekawa 104: }
1.1.1.2 ! ohara 105: else
! 106: {
! 107: int str_is_null;
! 108: int pow2;
! 109: mp_exp_t e, f;
! 110: mpfr_t a, b;
! 111: mpz_t bz;
! 112: char *str0;
! 113: size_t len;
! 114:
! 115: str_is_null = str == NULL;
! 116:
! 117: if (IS_POW2(base)) /* Is base a power of 2? */
! 118: {
! 119: count_leading_zeros (pow2, (mp_limb_t) base);
! 120: pow2 = BITS_PER_MP_LIMB - pow2 - 1; /* base = 2^pow2 */
! 121: }
! 122: else
! 123: pow2 = 0;
! 124:
! 125: /* first determines the exponent */
! 126: e = MPFR_EXP(op);
! 127:
! 128: /* the absolute value of op is between 1/2*2^e and 2^e */
! 129: /* the output exponent f is such that base^(f-1) <= |op| < base^f
! 130: i.e. f = 1 + floor(log(|op|)/log(base))
! 131: = 1 + floor((log(|m|)+e*log(2))/log(base)) */
! 132: /* f = 1 + (int)floor((log(d)/LOG2+(double)e)*LOG2/log((double)base)); */
! 133: /* when base = 2^pow2, then |op| < 2^(pow2*f)
! 134: i.e. e <= pow2*f and f = ceil(e/pow2) */
! 135: if (pow2)
! 136: f = e <= 0 ? e / pow2 : (e - 1) / pow2 + 1; /* int overflow avoided */
! 137: else
! 138: {
! 139: double d;
! 140:
! 141: d = mpfr_get_d3(op, 0, GMP_RNDN);
! 142: d = ((double) e + (double) _mpfr_floor_log2 (ABS(d)))
! 143: * __mp_bases[base].chars_per_bit_exactly;
! 144: MPFR_ASSERTN(d >= MP_EXP_T_MIN);
! 145: MPFR_ASSERTN(d <= MP_EXP_T_MAX);
! 146: /* warning: (mp_exp_t) d rounds towards 0 */
! 147: f = (mp_exp_t) d; /* f = floor(d) if d >= 0 and ceil(d) if d < 0 */
! 148: if (f <= d)
! 149: {
! 150: MPFR_ASSERTN(f < MP_EXP_T_MAX);
! 151: f++;
! 152: }
! 153: }
1.1 maekawa 154:
1.1.1.2 ! ohara 155: mpfr_init (a);
! 156: mpfr_init (b);
! 157: mpz_init (bz);
! 158:
! 159: str0 = str;
! 160:
! 161: do
! 162: {
! 163: mp_prec_t prec, q;
! 164:
! 165: /* now the first n digits of the mantissa are obtained from
! 166: rnd(op*base^(n-f)) */
! 167: if (pow2)
! 168: {
! 169: MPFR_ASSERTN(n <= MPFR_INTPREC_MAX / pow2);
! 170: prec = (mp_prec_t) n * pow2;
! 171: }
! 172: else
! 173: {
! 174: double d;
! 175:
! 176: d = (double) n / __mp_bases[base].chars_per_bit_exactly;
! 177: MPFR_ASSERTN(d <= MPFR_INTPREC_MAX - 1);
! 178: prec = (mp_prec_t) d + 1;
! 179: }
! 180:
! 181: MPFR_ASSERTN(prec <= MPFR_INTPREC_MAX - ERR);
! 182: /* one has to use at least q bits */
! 183: q = ((prec + (ERR-1)) / BITS_PER_MP_LIMB + 1) * BITS_PER_MP_LIMB;
! 184: mpfr_set_prec (a, q);
! 185: mpfr_set_prec (b, q);
! 186:
! 187: while (1)
! 188: {
! 189: mp_exp_unsigned_t p;
! 190: int div;
! 191:
! 192: if (f < 0)
! 193: {
! 194: p = (mp_exp_unsigned_t) n - f;
! 195: MPFR_ASSERTN(p > n);
! 196: div = 0;
! 197: }
! 198: else if (n >= f)
! 199: {
! 200: p = n - f;
! 201: div = 0;
! 202: }
! 203: else
! 204: {
! 205: p = f - n;
! 206: div = 1;
! 207: }
! 208:
! 209: if (pow2)
! 210: {
! 211: MPFR_ASSERTN(p <= ULONG_MAX / pow2);
! 212: MPFR_ASSERTN(p <= __mpfr_emax / pow2);
! 213: if (div)
! 214: mpfr_div_2ui (b, op, pow2*p, rnd_mode);
! 215: else
! 216: mpfr_mul_2ui (b, op, pow2*p, rnd_mode);
! 217: }
! 218: else
! 219: {
! 220: /* compute base^p with q bits */
! 221: mpfr_set_prec (b, q);
! 222: if (p == 0)
! 223: {
! 224: mpfr_set (b, op, rnd_mode);
! 225: mpfr_set_ui (a, 1, rnd_mode);
! 226: }
! 227: else
! 228: {
! 229: mp_rnd_t rnd1;
! 230:
! 231: mpfr_set_prec (a, q);
! 232: if (div)
! 233: {
! 234: /* if div, we divide by base^p, so we have to invert
! 235: the rounding mode to compute base^p */
! 236: switch (rnd_mode)
! 237: {
! 238: case GMP_RNDN: rnd1 = GMP_RNDN; break;
! 239: case GMP_RNDZ: rnd1 = GMP_RNDU; break;
! 240: case GMP_RNDU: rnd1 = GMP_RNDZ; break;
! 241: case GMP_RNDD: rnd1 = GMP_RNDU; break;
! 242: default: MPFR_ASSERTN(0);
! 243: }
! 244: }
! 245: else
! 246: {
! 247: rnd1 = rnd_mode;
! 248: }
! 249: mpfr_ui_pow_ui (a, base, p, rnd1);
! 250: if (div)
! 251: {
! 252: mpfr_set_ui (b, 1, rnd_mode);
! 253: mpfr_div (a, b, a, rnd_mode);
! 254: }
! 255: /* now a is an approximation to 1/base^(f-n) */
! 256: mpfr_mul (b, op, a, rnd_mode);
! 257: }
! 258: }
! 259:
! 260: if (neg)
! 261: MPFR_CHANGE_SIGN(b); /* put b positive */
! 262:
! 263: if (prec <= (MPFR_INTPREC_MAX - BITS_PER_MP_LIMB) / 2 &&
! 264: q > 2 * prec + BITS_PER_MP_LIMB)
! 265: {
! 266: /* if the intermediate precision exceeds twice that of the
! 267: input, a worst-case for the division cannot occur */
! 268: rnd_mode = GMP_RNDN;
! 269: break;
! 270: }
! 271: else if (pow2 ||
! 272: mpfr_can_round (b, q-ERR, rnd_mode, rnd_mode, prec))
! 273: break;
! 274:
! 275: MPFR_ASSERTN(q <= MPFR_INTPREC_MAX - BITS_PER_MP_LIMB);
! 276: q += BITS_PER_MP_LIMB;
! 277: }
! 278:
! 279: {
! 280: mp_rnd_t rnd = rnd_mode;
! 281:
! 282: if (neg)
! 283: switch (rnd_mode)
! 284: {
! 285: case GMP_RNDU: rnd = GMP_RNDZ; break;
! 286: case GMP_RNDD: rnd = GMP_RNDU; break;
! 287: }
! 288:
! 289: mpfr_round_prec (b, rnd, MPFR_EXP(b));
! 290: }
! 291:
! 292: prec = MPFR_EXP(b); /* may have changed due to rounding */
! 293:
! 294: {
! 295: mp_size_t n, e;
! 296: int sh;
! 297:
! 298: /* now the mantissa is the integer part of b */
! 299: n = 1 + (prec - 1) / BITS_PER_MP_LIMB;
! 300: _mpz_realloc (bz, n);
! 301: sh = prec % BITS_PER_MP_LIMB;
! 302: e = 1 + (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
! 303: MPFR_ASSERTN(e >= n);
! 304: e -= n;
! 305: if (sh != 0)
! 306: mpn_rshift (PTR(bz), MPFR_MANT(b) + e, n, BITS_PER_MP_LIMB - sh);
! 307: else
! 308: MPN_COPY(PTR(bz), MPFR_MANT(b) + e, n);
! 309: bz->_mp_size = n;
! 310: }
! 311:
! 312: /* computes the number of characters needed */
! 313: /* n+1 may not be enough for 100000... */
! 314: if (str0 == NULL)
! 315: str0 = (*__gmp_allocate_func) (neg + n + 2);
! 316:
! 317: str = str0; /* restore initial value in case we had to restart */
! 318:
! 319: if (neg)
! 320: *str++ = '-';
! 321:
! 322: mpz_get_str (str, base, bz); /* n digits of mantissa */
! 323: len = strlen(str);
! 324: }
! 325: while (len > n && (f++, 1));
! 326:
! 327: if (len == n - 1)
! 328: {
! 329: len++;
! 330: f--;
! 331: str[n-1]='0';
! 332: str[n]='\0';
! 333: }
! 334:
! 335: *expptr = f;
! 336: mpfr_clear (a);
! 337: mpfr_clear (b);
! 338: mpz_clear (bz);
! 339:
! 340: /* if the given string was null, ensure we return a block
! 341: which is exactly strlen(str0)+1 bytes long (useful for
! 342: __gmp_free_func and the C++ wrapper) */
! 343:
! 344: /* NOTE: len != n + 1 is always satisfied; either this condition
! 345: is useless or there is a bug somewhere */
! 346: if (str_is_null && len != n + 1)
! 347: str0 = (*__gmp_reallocate_func) (str0, neg + n + 2, neg + len + 1);
! 348:
! 349: return str0;
! 350: }
! 351: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to get_str.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr