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

Diff for /OpenXM_contrib/gmp/mpfr/Attic/get_str.c between version 1.1.1.1 and 1.1.1.2

version 1.1.1.1, 2000/09/09 14:12:19 version 1.1.1.2, 2003/08/25 16:06:07
Line 1 
Line 1 
 /* mpfr_get_str -- output a floating-point number to a string  /* mpfr_get_str -- output a floating-point number to a string
   
 Copyright (C) 1999 PolKA project, Inria Lorraine and Loria  Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
   
 This file is part of the MPFR Library.  This file is part of the MPFR Library.
   
 The MPFR Library is free software; you can redistribute it and/or modify  The MPFR Library is free software; you can redistribute it and/or modify
 it under the terms of the GNU Library General Public License as published by  it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation; either version 2 of the License, or (at your  the Free Software Foundation; either version 2.1 of the License, or (at your
 option) any later version.  option) any later version.
   
 The MPFR Library is distributed in the hope that it will be useful, but  The MPFR Library is distributed in the hope that it will be useful, but
 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 License for more details.  License for more details.
   
 You should have received a copy of the GNU Library General Public License  You should have received a copy of the GNU Lesser General Public License
 along with the MPFR Library; see the file COPYING.LIB.  If not, write to  along with the MPFR Library; see the file COPYING.LIB.  If not, write to
 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,  the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 MA 02111-1307, USA. */  MA 02111-1307, USA. */
   
 #include <math.h>  
 #include <stdio.h>  
 #include <stdlib.h>  #include <stdlib.h>
 #include <string.h>  #include <string.h>
 #include "gmp.h"  #include "gmp.h"
 #include "gmp-impl.h"  #include "gmp-impl.h"
 #include "longlong.h"  #include "longlong.h"
 #include "mpfr.h"  #include "mpfr.h"
   #include "mpfr-impl.h"
   
   #define ERR 5
   
 /*  /*
   Convert op to a string in base 'base' with 'n' digits and writes the    Convert op to a string in base 'base' with 'n' digits and writes the
   mantissa in 'str', the exponent in 'expptr'.    mantissa in 'str', the exponent in 'expptr'.
Line 35  MA 02111-1307, USA. */
Line 36  MA 02111-1307, USA. */
   
   For op = 3.1416 we get str = "31416" and expptr=1.    For op = 3.1416 we get str = "31416" and expptr=1.
  */   */
 #if __STDC__  char*
 char *mpfr_get_str(char *str, mp_exp_t *expptr, int base, size_t n,  mpfr_get_str (char *str, mp_exp_t *expptr, int base, size_t n,
                   mpfr_srcptr op, unsigned char rnd_mode)                mpfr_srcptr op, mp_rnd_t rnd_mode)
 #else  
 char *mpfr_get_str(str, expptr, base, n, op, rnd_mode)  
      char *str;  
      mp_exp_t *expptr;  
      int base;  
      size_t n;  
      mpfr_srcptr op;  
      unsigned char rnd_mode;  
 #endif  
 {  {
   double d; long e, q, div, p, err, prec, sh; mpfr_t a, b; mpz_t bz;    int neg;
   char *str0; unsigned char rnd1; int f, pow2, ok=0, neg;  
   
   if (base<2 || 36<base) {    if (base < 2 || base > 36)
     fprintf(stderr, "Error: too small or too large base in mpfr_get_str: %d\n",      return NULL;
             base);  
     exit(1);  
   }  
   
   neg = (SIGN(op)<0) ? 1 : 0;    if (n == 0) /* determine n from precision of op */
       {
         n = __mp_bases[base].chars_per_bit_exactly * MPFR_PREC(op);
         if (n < 2)
           n = 2;
       }
   
   if (!NOTZERO(op)) {    /* Do not use MPFR_PREC_MIN as this applies to base 2 only. Perhaps we
     if (str==NULL) str0=str=(*_mp_allocate_func)(neg + n + 2);       should allow n == 1 for directed rounding modes and odd bases. */
     if (SIGN(op)<0) *str++ = '-';    if (n < 2)
     for (f=0;f<n;f++) *str++ = '0';      return NULL;
     *expptr = 1;  
     return str0;  
   }  
   
   count_leading_zeros(pow2, (mp_limb_t)base);    if (MPFR_IS_NAN(op))
   pow2 = BITS_PER_MP_LIMB - pow2 - 1;      {
   if (base != (1<<pow2)) pow2=0;        if (str == NULL)
   /* if pow2 <> 0, then base = 2^pow2 */          str = (*__gmp_allocate_func)(4);
         str[0] = 'N';
         str[1] = 'a';
         str[2] = 'N';
         str[3] = '\0';
         return str;
       }
   
   /* first determines the exponent */    neg = MPFR_SIGN(op) < 0; /* 0 if positive, 1 if negative */
   e = EXP(op);  
   d = fabs(mpfr_get_d2(op, 0));  
   /* the absolute value of op is between 1/2*2^e and 2^e */  
   /* the output exponent f is such that base^(f-1) <= |op| < base^f  
      i.e. f = 1 + floor(log(|op|)/log(base))  
      = 1 + floor((log(|m|)+e*log(2))/log(base)) */  
   f = 1 + (int) floor((log(d)+(double)e*log(2.0))/log((double)base));  
   if (n==0) {  
     /* performs exact rounding, i.e. returns y such that for rnd_mode=RNDN  
        for example, we have:  
        y*base^(f-n) <= x*2^(e-p) < (x+1)*2^(e-p) <= (y+1)*base^(f-n)  
        which implies 2^(EXP(op)-PREC(op)) <= base^(f-n)  
      */  
     n = f + (int) ceil(((double)PREC(op)-e)*log(2.0)/log((double)base));  
   }  
   /* now the first n digits of the mantissa are obtained from  
      rnd(op*base^(n-f)) */  
   prec = (long) ceil((double)n*log((double)base)/log(2.0));  
   err = 5;  
   q = prec+err;  
   /* one has to use at least q bits */  
   q = (((q-1)/BITS_PER_MP_LIMB)+1)*BITS_PER_MP_LIMB;  
   mpfr_init2(a,q); mpfr_init2(b,q);  
   
   do {    if (MPFR_IS_INF(op))
     p = n-f; if ((div=(p<0))) p=-p;      {
     rnd1 = rnd_mode;        char *str0;
     if (div) {  
       /* if div we divide by base^p so we have to invert the rounding mode */        if (str == NULL)
       switch (rnd1) {          str = (*__gmp_allocate_func)(neg + 4);
       case GMP_RNDN: rnd1=GMP_RNDN; break;        str0 = str;
       case GMP_RNDZ: rnd1=GMP_RNDU; break;        if (neg)
       case GMP_RNDU: rnd1=GMP_RNDZ; break;          *str++ = '-';
       case GMP_RNDD: rnd1=GMP_RNDZ; break;        *str++ = 'I';
       }        *str++ = 'n';
         *str++ = 'f';
         *str = '\0';
         return str0; /* strlen(str0) = neg + 3 */
     }      }
   
     if (pow2) {    /* op is a floating-point number */
       if (div) mpfr_div_2exp(b, op, pow2*p, rnd_mode);  
       else mpfr_mul_2exp(b, op, pow2*p, rnd_mode);    if (MPFR_IS_ZERO(op))
     }      {
     else {        char *str0;
        /* compute base^p with q bits and rounding towards zero */  
        mpfr_set_prec(b, q);        if (str == NULL)
        if (p==0) { mpfr_set(b, op, rnd_mode); mpfr_set_ui(a, 1, rnd_mode); }          str = (*__gmp_allocate_func)(neg + n + 1);
        else {        str0 = str;
          mpfr_set_prec(a, q);        if (neg)
          mpfr_ui_pow_ui(a, base, p, rnd1);          *str++ = '-';
          if (div) {        memset(str, '0', n);
            mpfr_set_ui(b, 1, rnd_mode);        str[n] = '\0';
            mpfr_div(a, b, a, rnd_mode);        *expptr = 0; /* a bit like frexp() in ISO C99 */
          }        return str0; /* strlen(str0) = neg + n */
          /* now a is an approximation by default of 1/base^(f-n) */  
          mpfr_mul(b, op, a, rnd_mode);  
        }  
     }      }
     if (neg) CHANGE_SIGN(b); /* put b positive */    else
       {
         int str_is_null;
         int pow2;
         mp_exp_t e, f;
         mpfr_t a, b;
         mpz_t bz;
         char *str0;
         size_t len;
   
     if (q>2*prec+BITS_PER_MP_LIMB) {        str_is_null = str == NULL;
       /* happens when just in the middle between two digits */  
       n--; q-=BITS_PER_MP_LIMB;        if (IS_POW2(base)) /* Is base a power of 2? */
       if (n==0) {          {
           fprintf(stderr, "cannot determine leading digit\n"); exit(1);            count_leading_zeros (pow2, (mp_limb_t) base);
             pow2 = BITS_PER_MP_LIMB - pow2 - 1; /* base = 2^pow2 */
         }          }
     }        else
     ok = pow2 || mpfr_can_round(b, q-err, rnd_mode, rnd_mode, prec);          pow2 = 0;
   
     if (ok) {        /* first determines the exponent */
       if (pow2) {        e = MPFR_EXP(op);
         sh = e-PREC(op) + pow2*(n-f); /* error at most 2^e */  
         ok = mpfr_can_round(b, EXP(b)-sh-1, rnd_mode, rnd_mode, n*pow2);  
       }  
       else {  
          /* check that value is the same at distance 2^(e-PREC(op))/base^(f-n)  
           in opposite from rounding direction */  
          if (e>=PREC(op)) mpfr_mul_2exp(a, a, e-PREC(op), rnd_mode);  
          else mpfr_div_2exp(a, a, PREC(op)-e, rnd_mode);  
          if (rnd_mode==GMP_RNDN) {  
            mpfr_div_2exp(a, a, 2, rnd_mode);  
            mpfr_sub(b, b, a, rnd_mode); /* b - a/2 */  
            mpfr_mul_2exp(a, a, 2, rnd_mode);  
            mpfr_add(a, b, a, rnd_mode); /* b + a/2 */  
          }  
          else if ((rnd_mode==GMP_RNDU && neg==0) || (rnd_mode==GMP_RNDD && neg))  
            mpfr_sub(a, b, a, rnd_mode);  
          else mpfr_add(a, b, a, rnd_mode);  
          /* check that a and b are rounded similarly */  
          prec=EXP(b);  
          if (EXP(a) != prec) ok=0;  
          else {  
            mpfr_round(b, rnd_mode, prec);  
            mpfr_round(a, rnd_mode, prec);  
            if (mpfr_cmp(a, b)) ok=0;  
          }  
        }  
       if (ok==0) { /* n is too large */  
         n--;  
         if (n==0) {  
           fprintf(stderr, "cannot determine leading digit\n"); exit(1);  
         }  
         q -= BITS_PER_MP_LIMB;  
       }  
     }  
   } while (ok==0 && (q+=BITS_PER_MP_LIMB) );  
   if (neg)  
     switch (rnd_mode) {  
     case GMP_RNDU: rnd_mode=GMP_RNDZ; break;  
     case GMP_RNDD: rnd_mode=GMP_RNDU; break;  
   }  
   
   prec=EXP(b); /* may have changed due to rounding */        /* the absolute value of op is between 1/2*2^e and 2^e */
         /* the output exponent f is such that base^(f-1) <= |op| < base^f
            i.e. f = 1 + floor(log(|op|)/log(base))
            = 1 + floor((log(|m|)+e*log(2))/log(base)) */
         /* f = 1 + (int)floor((log(d)/LOG2+(double)e)*LOG2/log((double)base)); */
         /* when base = 2^pow2, then |op| < 2^(pow2*f)
            i.e. e <= pow2*f and f = ceil(e/pow2) */
         if (pow2)
           f = e <= 0 ? e / pow2 : (e - 1) / pow2 + 1; /* int overflow avoided */
         else
           {
             double d;
   
   /* now the mantissa is the integer part of b */            d = mpfr_get_d3(op, 0, GMP_RNDN);
   mpz_init(bz); q=1+(prec-1)/BITS_PER_MP_LIMB;            d = ((double) e + (double) _mpfr_floor_log2 (ABS(d)))
   _mpz_realloc(bz, q);              * __mp_bases[base].chars_per_bit_exactly;
   sh = prec%BITS_PER_MP_LIMB;            MPFR_ASSERTN(d >= MP_EXP_T_MIN);
   e = 1 + (PREC(b)-1)/BITS_PER_MP_LIMB-q;            MPFR_ASSERTN(d <= MP_EXP_T_MAX);
   if (sh) mpn_rshift(PTR(bz), MANT(b)+e, q, BITS_PER_MP_LIMB-sh);            /* warning: (mp_exp_t) d rounds towards 0 */
   else MPN_COPY(PTR(bz), MANT(b)+e, q);            f = (mp_exp_t) d; /* f = floor(d) if d >= 0 and ceil(d) if d < 0 */
   bz->_mp_size=q;            if (f <= d)
               {
                 MPFR_ASSERTN(f < MP_EXP_T_MAX);
                 f++;
               }
           }
   
   /* computes the number of characters needed */        mpfr_init (a);
   q = neg + n + 2; /* n+1 may not be enough for 100000... */        mpfr_init (b);
   if (str==NULL) str0=str=(*_mp_allocate_func)(q);        mpz_init (bz);
   if (neg) *str++='-';  
   mpz_get_str(str, base, bz); /* n digits of mantissa */  
   if (strlen(str)==n+1) f++; /* possible due to rounding */  
   *expptr = f;  
   mpfr_clear(a); mpfr_clear(b); mpz_clear(bz);  
   return str0;  
 }  
   
         str0 = str;
   
         do
           {
             mp_prec_t prec, q;
   
             /* now the first n digits of the mantissa are obtained from
                rnd(op*base^(n-f)) */
             if (pow2)
               {
                 MPFR_ASSERTN(n <= MPFR_INTPREC_MAX / pow2);
                 prec = (mp_prec_t) n * pow2;
               }
             else
               {
                 double d;
   
                 d = (double) n / __mp_bases[base].chars_per_bit_exactly;
                 MPFR_ASSERTN(d <= MPFR_INTPREC_MAX - 1);
                 prec = (mp_prec_t) d + 1;
               }
   
             MPFR_ASSERTN(prec <= MPFR_INTPREC_MAX - ERR);
             /* one has to use at least q bits */
             q = ((prec + (ERR-1)) / BITS_PER_MP_LIMB + 1) * BITS_PER_MP_LIMB;
             mpfr_set_prec (a, q);
             mpfr_set_prec (b, q);
   
             while (1)
               {
                 mp_exp_unsigned_t p;
                 int div;
   
                 if (f < 0)
                   {
                     p = (mp_exp_unsigned_t) n - f;
                     MPFR_ASSERTN(p > n);
                     div = 0;
                   }
                 else if (n >= f)
                   {
                     p = n - f;
                     div = 0;
                   }
                 else
                   {
                     p = f - n;
                     div = 1;
                   }
   
                 if (pow2)
                   {
                     MPFR_ASSERTN(p <= ULONG_MAX / pow2);
                     MPFR_ASSERTN(p <= __mpfr_emax / pow2);
                     if (div)
                       mpfr_div_2ui (b, op, pow2*p, rnd_mode);
                     else
                       mpfr_mul_2ui (b, op, pow2*p, rnd_mode);
                   }
                 else
                   {
                     /* compute base^p with q bits */
                     mpfr_set_prec (b, q);
                     if (p == 0)
                       {
                         mpfr_set (b, op, rnd_mode);
                         mpfr_set_ui (a, 1, rnd_mode);
                       }
                     else
                       {
                         mp_rnd_t rnd1;
   
                         mpfr_set_prec (a, q);
                         if (div)
                           {
                             /* if div, we divide by base^p, so we have to invert
                                the rounding mode to compute base^p */
                             switch (rnd_mode)
                               {
                               case GMP_RNDN: rnd1 = GMP_RNDN; break;
                               case GMP_RNDZ: rnd1 = GMP_RNDU; break;
                               case GMP_RNDU: rnd1 = GMP_RNDZ; break;
                               case GMP_RNDD: rnd1 = GMP_RNDU; break;
                               default: MPFR_ASSERTN(0);
                               }
                           }
                         else
                           {
                             rnd1 = rnd_mode;
                           }
                         mpfr_ui_pow_ui (a, base, p, rnd1);
                         if (div)
                           {
                             mpfr_set_ui (b, 1, rnd_mode);
                             mpfr_div (a, b, a, rnd_mode);
                           }
                         /* now a is an approximation to 1/base^(f-n) */
                         mpfr_mul (b, op, a, rnd_mode);
                       }
                   }
   
                 if (neg)
                   MPFR_CHANGE_SIGN(b); /* put b positive */
   
                 if (prec <= (MPFR_INTPREC_MAX - BITS_PER_MP_LIMB) / 2 &&
                     q > 2 * prec + BITS_PER_MP_LIMB)
                   {
                     /* if the intermediate precision exceeds twice that of the
                        input, a worst-case for the division cannot occur */
                     rnd_mode = GMP_RNDN;
                     break;
                   }
                 else if (pow2 ||
                          mpfr_can_round (b, q-ERR, rnd_mode, rnd_mode, prec))
                   break;
   
                 MPFR_ASSERTN(q <= MPFR_INTPREC_MAX - BITS_PER_MP_LIMB);
                 q += BITS_PER_MP_LIMB;
               }
   
             {
               mp_rnd_t rnd = rnd_mode;
   
               if (neg)
                 switch (rnd_mode)
                   {
                   case GMP_RNDU: rnd = GMP_RNDZ; break;
                   case GMP_RNDD: rnd = GMP_RNDU; break;
                   }
   
               mpfr_round_prec (b, rnd, MPFR_EXP(b));
             }
   
             prec = MPFR_EXP(b); /* may have changed due to rounding */
   
             {
               mp_size_t n, e;
               int sh;
   
               /* now the mantissa is the integer part of b */
               n = 1 + (prec - 1) / BITS_PER_MP_LIMB;
               _mpz_realloc (bz, n);
               sh = prec % BITS_PER_MP_LIMB;
               e = 1 + (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
               MPFR_ASSERTN(e >= n);
               e -= n;
               if (sh != 0)
                 mpn_rshift (PTR(bz), MPFR_MANT(b) + e, n, BITS_PER_MP_LIMB - sh);
               else
                 MPN_COPY(PTR(bz), MPFR_MANT(b) + e, n);
               bz->_mp_size = n;
             }
   
             /* computes the number of characters needed */
             /* n+1 may not be enough for 100000... */
             if (str0 == NULL)
               str0 = (*__gmp_allocate_func) (neg + n + 2);
   
             str = str0; /* restore initial value in case we had to restart */
   
             if (neg)
               *str++ = '-';
   
             mpz_get_str (str, base, bz); /* n digits of mantissa */
             len = strlen(str);
           }
         while (len > n && (f++, 1));
   
         if (len == n - 1)
           {
             len++;
             f--;
             str[n-1]='0';
             str[n]='\0';
           }
   
         *expptr = f;
         mpfr_clear (a);
         mpfr_clear (b);
         mpz_clear (bz);
   
         /* if the given string was null, ensure we return a block
            which is exactly strlen(str0)+1 bytes long (useful for
            __gmp_free_func and the C++ wrapper) */
   
         /* NOTE: len != n + 1 is always satisfied; either this condition
            is useless or there is a bug somewhere */
         if (str_is_null && len != n + 1)
           str0 = (*__gmp_reallocate_func) (str0, neg + n + 2, neg + len + 1);
   
         return str0;
       }
   }
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.1.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>