=================================================================== RCS file: /home/cvs/OpenXM_contrib/gmp/mpfr/Attic/get_str.c,v retrieving revision 1.1.1.1 retrieving revision 1.1.1.2 diff -u -p -r1.1.1.1 -r1.1.1.2 --- OpenXM_contrib/gmp/mpfr/Attic/get_str.c 2000/09/09 14:12:19 1.1.1.1 +++ OpenXM_contrib/gmp/mpfr/Attic/get_str.c 2003/08/25 16:06:07 1.1.1.2 @@ -1,33 +1,34 @@ /* 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. 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 -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. The MPFR Library is distributed in the hope that it will be useful, but 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. -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 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include -#include #include #include #include "gmp.h" #include "gmp-impl.h" #include "longlong.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 mantissa in 'str', the exponent in 'expptr'. @@ -35,173 +36,316 @@ MA 02111-1307, USA. */ For op = 3.1416 we get str = "31416" and expptr=1. */ -#if __STDC__ -char *mpfr_get_str(char *str, mp_exp_t *expptr, int base, size_t n, - mpfr_srcptr op, unsigned char 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 +char* +mpfr_get_str (char *str, mp_exp_t *expptr, int base, size_t n, + mpfr_srcptr op, mp_rnd_t rnd_mode) { - double d; long e, q, div, p, err, prec, sh; mpfr_t a, b; mpz_t bz; - char *str0; unsigned char rnd1; int f, pow2, ok=0, neg; + int neg; - if (base<2 || 36 36) + return NULL; - 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)) { - if (str==NULL) str0=str=(*_mp_allocate_func)(neg + n + 2); - if (SIGN(op)<0) *str++ = '-'; - for (f=0;f 0, then base = 2^pow2 */ + if (MPFR_IS_NAN(op)) + { + if (str == NULL) + str = (*__gmp_allocate_func)(4); + str[0] = 'N'; + str[1] = 'a'; + str[2] = 'N'; + str[3] = '\0'; + return str; + } - /* first determines the exponent */ - 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); + neg = MPFR_SIGN(op) < 0; /* 0 if positive, 1 if negative */ - do { - p = n-f; if ((div=(p<0))) p=-p; - rnd1 = rnd_mode; - if (div) { - /* if div we divide by base^p so we have to invert the rounding mode */ - switch (rnd1) { - 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_RNDZ; break; - } + if (MPFR_IS_INF(op)) + { + char *str0; + + if (str == NULL) + str = (*__gmp_allocate_func)(neg + 4); + str0 = str; + if (neg) + *str++ = '-'; + *str++ = 'I'; + *str++ = 'n'; + *str++ = 'f'; + *str = '\0'; + return str0; /* strlen(str0) = neg + 3 */ } - if (pow2) { - if (div) mpfr_div_2exp(b, op, pow2*p, rnd_mode); - else mpfr_mul_2exp(b, op, pow2*p, rnd_mode); - } - else { - /* compute base^p with q bits and rounding towards zero */ - mpfr_set_prec(b, q); - if (p==0) { mpfr_set(b, op, rnd_mode); mpfr_set_ui(a, 1, rnd_mode); } - else { - mpfr_set_prec(a, q); - 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 by default of 1/base^(f-n) */ - mpfr_mul(b, op, a, rnd_mode); - } + /* op is a floating-point number */ + + if (MPFR_IS_ZERO(op)) + { + char *str0; + + if (str == NULL) + str = (*__gmp_allocate_func)(neg + n + 1); + str0 = str; + if (neg) + *str++ = '-'; + memset(str, '0', n); + str[n] = '\0'; + *expptr = 0; /* a bit like frexp() in ISO C99 */ + return str0; /* strlen(str0) = neg + n */ } - 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) { - /* happens when just in the middle between two digits */ - n--; q-=BITS_PER_MP_LIMB; - if (n==0) { - fprintf(stderr, "cannot determine leading digit\n"); exit(1); + str_is_null = str == NULL; + + if (IS_POW2(base)) /* Is base a power of 2? */ + { + count_leading_zeros (pow2, (mp_limb_t) base); + pow2 = BITS_PER_MP_LIMB - pow2 - 1; /* base = 2^pow2 */ } - } - ok = pow2 || mpfr_can_round(b, q-err, rnd_mode, rnd_mode, prec); + else + pow2 = 0; - if (ok) { - if (pow2) { - 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; - } + /* first determines the exponent */ + e = MPFR_EXP(op); - 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 */ - mpz_init(bz); q=1+(prec-1)/BITS_PER_MP_LIMB; - _mpz_realloc(bz, q); - sh = prec%BITS_PER_MP_LIMB; - e = 1 + (PREC(b)-1)/BITS_PER_MP_LIMB-q; - if (sh) mpn_rshift(PTR(bz), MANT(b)+e, q, BITS_PER_MP_LIMB-sh); - else MPN_COPY(PTR(bz), MANT(b)+e, q); - bz->_mp_size=q; + d = mpfr_get_d3(op, 0, GMP_RNDN); + d = ((double) e + (double) _mpfr_floor_log2 (ABS(d))) + * __mp_bases[base].chars_per_bit_exactly; + MPFR_ASSERTN(d >= MP_EXP_T_MIN); + MPFR_ASSERTN(d <= MP_EXP_T_MAX); + /* warning: (mp_exp_t) d rounds towards 0 */ + f = (mp_exp_t) d; /* f = floor(d) if d >= 0 and ceil(d) if d < 0 */ + if (f <= d) + { + MPFR_ASSERTN(f < MP_EXP_T_MAX); + f++; + } + } - /* computes the number of characters needed */ - q = neg + n + 2; /* n+1 may not be enough for 100000... */ - if (str==NULL) str0=str=(*_mp_allocate_func)(q); - 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; -} + mpfr_init (a); + mpfr_init (b); + mpz_init (bz); + 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; + } +}