Annotation of OpenXM_contrib/gmp/mpfr/get_str.c, Revision 1.1
1.1 ! maekawa 1: /* mpfr_get_str -- output a floating-point number to a string
! 2:
! 3: Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
! 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
! 8: it under the terms of the GNU Library General Public License as published by
! 9: the Free Software Foundation; either version 2 of the License, or (at your
! 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
! 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
! 15: License for more details.
! 16:
! 17: You should have received a copy of the GNU Library General Public License
! 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 <math.h>
! 23: #include <stdio.h>
! 24: #include <stdlib.h>
! 25: #include <string.h>
! 26: #include "gmp.h"
! 27: #include "gmp-impl.h"
! 28: #include "longlong.h"
! 29: #include "mpfr.h"
! 30:
! 31: /*
! 32: Convert op to a string in base 'base' with 'n' digits and writes the
! 33: mantissa in 'str', the exponent in 'expptr'.
! 34: The result is rounded wrt 'rnd_mode'.
! 35:
! 36: For op = 3.1416 we get str = "31416" and expptr=1.
! 37: */
! 38: #if __STDC__
! 39: char *mpfr_get_str(char *str, mp_exp_t *expptr, int base, size_t n,
! 40: mpfr_srcptr op, unsigned char rnd_mode)
! 41: #else
! 42: char *mpfr_get_str(str, expptr, base, n, op, rnd_mode)
! 43: char *str;
! 44: mp_exp_t *expptr;
! 45: int base;
! 46: size_t n;
! 47: mpfr_srcptr op;
! 48: unsigned char rnd_mode;
! 49: #endif
! 50: {
! 51: double d; long e, q, div, p, err, prec, sh; mpfr_t a, b; mpz_t bz;
! 52: char *str0; unsigned char rnd1; int f, pow2, ok=0, neg;
! 53:
! 54: if (base<2 || 36<base) {
! 55: fprintf(stderr, "Error: too small or too large base in mpfr_get_str: %d\n",
! 56: base);
! 57: exit(1);
! 58: }
! 59:
! 60: neg = (SIGN(op)<0) ? 1 : 0;
! 61:
! 62: if (!NOTZERO(op)) {
! 63: if (str==NULL) str0=str=(*_mp_allocate_func)(neg + n + 2);
! 64: if (SIGN(op)<0) *str++ = '-';
! 65: for (f=0;f<n;f++) *str++ = '0';
! 66: *expptr = 1;
! 67: return str0;
! 68: }
! 69:
! 70: count_leading_zeros(pow2, (mp_limb_t)base);
! 71: pow2 = BITS_PER_MP_LIMB - pow2 - 1;
! 72: if (base != (1<<pow2)) pow2=0;
! 73: /* if pow2 <> 0, then base = 2^pow2 */
! 74:
! 75: /* first determines the exponent */
! 76: e = EXP(op);
! 77: d = fabs(mpfr_get_d2(op, 0));
! 78: /* the absolute value of op is between 1/2*2^e and 2^e */
! 79: /* the output exponent f is such that base^(f-1) <= |op| < base^f
! 80: i.e. f = 1 + floor(log(|op|)/log(base))
! 81: = 1 + floor((log(|m|)+e*log(2))/log(base)) */
! 82: f = 1 + (int) floor((log(d)+(double)e*log(2.0))/log((double)base));
! 83: if (n==0) {
! 84: /* performs exact rounding, i.e. returns y such that for rnd_mode=RNDN
! 85: for example, we have:
! 86: y*base^(f-n) <= x*2^(e-p) < (x+1)*2^(e-p) <= (y+1)*base^(f-n)
! 87: which implies 2^(EXP(op)-PREC(op)) <= base^(f-n)
! 88: */
! 89: n = f + (int) ceil(((double)PREC(op)-e)*log(2.0)/log((double)base));
! 90: }
! 91: /* now the first n digits of the mantissa are obtained from
! 92: rnd(op*base^(n-f)) */
! 93: prec = (long) ceil((double)n*log((double)base)/log(2.0));
! 94: err = 5;
! 95: q = prec+err;
! 96: /* one has to use at least q bits */
! 97: q = (((q-1)/BITS_PER_MP_LIMB)+1)*BITS_PER_MP_LIMB;
! 98: mpfr_init2(a,q); mpfr_init2(b,q);
! 99:
! 100: do {
! 101: p = n-f; if ((div=(p<0))) p=-p;
! 102: rnd1 = rnd_mode;
! 103: if (div) {
! 104: /* if div we divide by base^p so we have to invert the rounding mode */
! 105: switch (rnd1) {
! 106: case GMP_RNDN: rnd1=GMP_RNDN; break;
! 107: case GMP_RNDZ: rnd1=GMP_RNDU; break;
! 108: case GMP_RNDU: rnd1=GMP_RNDZ; break;
! 109: case GMP_RNDD: rnd1=GMP_RNDZ; break;
! 110: }
! 111: }
! 112:
! 113: if (pow2) {
! 114: if (div) mpfr_div_2exp(b, op, pow2*p, rnd_mode);
! 115: else mpfr_mul_2exp(b, op, pow2*p, rnd_mode);
! 116: }
! 117: else {
! 118: /* compute base^p with q bits and rounding towards zero */
! 119: mpfr_set_prec(b, q);
! 120: if (p==0) { mpfr_set(b, op, rnd_mode); mpfr_set_ui(a, 1, rnd_mode); }
! 121: else {
! 122: mpfr_set_prec(a, q);
! 123: mpfr_ui_pow_ui(a, base, p, rnd1);
! 124: if (div) {
! 125: mpfr_set_ui(b, 1, rnd_mode);
! 126: mpfr_div(a, b, a, rnd_mode);
! 127: }
! 128: /* now a is an approximation by default of 1/base^(f-n) */
! 129: mpfr_mul(b, op, a, rnd_mode);
! 130: }
! 131: }
! 132: if (neg) CHANGE_SIGN(b); /* put b positive */
! 133:
! 134: if (q>2*prec+BITS_PER_MP_LIMB) {
! 135: /* happens when just in the middle between two digits */
! 136: n--; q-=BITS_PER_MP_LIMB;
! 137: if (n==0) {
! 138: fprintf(stderr, "cannot determine leading digit\n"); exit(1);
! 139: }
! 140: }
! 141: ok = pow2 || mpfr_can_round(b, q-err, rnd_mode, rnd_mode, prec);
! 142:
! 143: if (ok) {
! 144: if (pow2) {
! 145: sh = e-PREC(op) + pow2*(n-f); /* error at most 2^e */
! 146: ok = mpfr_can_round(b, EXP(b)-sh-1, rnd_mode, rnd_mode, n*pow2);
! 147: }
! 148: else {
! 149: /* check that value is the same at distance 2^(e-PREC(op))/base^(f-n)
! 150: in opposite from rounding direction */
! 151: if (e>=PREC(op)) mpfr_mul_2exp(a, a, e-PREC(op), rnd_mode);
! 152: else mpfr_div_2exp(a, a, PREC(op)-e, rnd_mode);
! 153: if (rnd_mode==GMP_RNDN) {
! 154: mpfr_div_2exp(a, a, 2, rnd_mode);
! 155: mpfr_sub(b, b, a, rnd_mode); /* b - a/2 */
! 156: mpfr_mul_2exp(a, a, 2, rnd_mode);
! 157: mpfr_add(a, b, a, rnd_mode); /* b + a/2 */
! 158: }
! 159: else if ((rnd_mode==GMP_RNDU && neg==0) || (rnd_mode==GMP_RNDD && neg))
! 160: mpfr_sub(a, b, a, rnd_mode);
! 161: else mpfr_add(a, b, a, rnd_mode);
! 162: /* check that a and b are rounded similarly */
! 163: prec=EXP(b);
! 164: if (EXP(a) != prec) ok=0;
! 165: else {
! 166: mpfr_round(b, rnd_mode, prec);
! 167: mpfr_round(a, rnd_mode, prec);
! 168: if (mpfr_cmp(a, b)) ok=0;
! 169: }
! 170: }
! 171: if (ok==0) { /* n is too large */
! 172: n--;
! 173: if (n==0) {
! 174: fprintf(stderr, "cannot determine leading digit\n"); exit(1);
! 175: }
! 176: q -= BITS_PER_MP_LIMB;
! 177: }
! 178: }
! 179: } while (ok==0 && (q+=BITS_PER_MP_LIMB) );
! 180: if (neg)
! 181: switch (rnd_mode) {
! 182: case GMP_RNDU: rnd_mode=GMP_RNDZ; break;
! 183: case GMP_RNDD: rnd_mode=GMP_RNDU; break;
! 184: }
! 185:
! 186: prec=EXP(b); /* may have changed due to rounding */
! 187:
! 188: /* now the mantissa is the integer part of b */
! 189: mpz_init(bz); q=1+(prec-1)/BITS_PER_MP_LIMB;
! 190: _mpz_realloc(bz, q);
! 191: sh = prec%BITS_PER_MP_LIMB;
! 192: e = 1 + (PREC(b)-1)/BITS_PER_MP_LIMB-q;
! 193: if (sh) mpn_rshift(PTR(bz), MANT(b)+e, q, BITS_PER_MP_LIMB-sh);
! 194: else MPN_COPY(PTR(bz), MANT(b)+e, q);
! 195: bz->_mp_size=q;
! 196:
! 197: /* computes the number of characters needed */
! 198: q = neg + n + 2; /* n+1 may not be enough for 100000... */
! 199: if (str==NULL) str0=str=(*_mp_allocate_func)(q);
! 200: if (neg) *str++='-';
! 201: mpz_get_str(str, base, bz); /* n digits of mantissa */
! 202: if (strlen(str)==n+1) f++; /* possible due to rounding */
! 203: *expptr = f;
! 204: mpfr_clear(a); mpfr_clear(b); mpz_clear(bz);
! 205: return str0;
! 206: }
! 207:
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