version 1.1.1.1, 2000/09/09 14:12:19 |
version 1.1.1.2, 2003/08/25 16:06:07 |
|
|
/* 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; |
|
} |
|
} |