[BACK]Return to get_str.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gmp / mpfr

Annotation of OpenXM_contrib/gmp/mpfr/get_str.c, Revision 1.1.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>