From 80740eeedfe48594b6f18f56f6b6d92e37d4aae8 Mon Sep 17 00:00:00 2001 From: zimmerma Date: Wed, 24 Jul 2002 16:28:21 +0000 Subject: completely new version, written by Alain Delplanque and Paul Zimmermann. It now directly uses mpn_get_str, with subquadratic complexity. About 3 times faster than previous version in most cases. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@1993 280ebfd0-de03-0410-8827-d642c229c3f4 --- get_str.c | 898 +++++++++++++++++++++++++++++++++++++++++++------------------- 1 file changed, 623 insertions(+), 275 deletions(-) (limited to 'get_str.c') diff --git a/get_str.c b/get_str.c index e5799a2fb..600f48af6 100644 --- a/get_str.c +++ b/get_str.c @@ -1,6 +1,7 @@ /* mpfr_get_str -- output a floating-point number to a string Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc. +This function was contributed by Alain Delplanque and Paul Zimmermann. This file is part of the MPFR Library. @@ -19,333 +20,680 @@ 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 +static double _mpfr_ceil _PROTO ((double)); +static long mpn_exp _PROTO ((mp_limb_t *, mp_exp_t *, int, mp_exp_t, size_t)); +static int mpfr_get_str_aux _PROTO ((char *, mp_exp_t *, mp_limb_t *, + mp_size_t, mp_exp_t, long, int, size_t, mp_rnd_t)); + +static char num_to_text[] = "0123456789abcdefghijklmnopqrstuvwxyz"; + +/* log_b2[b-2] = floor(log(2)/log(b)) for 2 <= b <= 36 */ +static const double log_b2[35] = + {0.10000000000000000000E1, + 0.63092975357145743709E0, + 0.50000000000000000000E0, + 0.43067655807339305067E0, + 0.38685280723454158687E0, + 0.35620718710802217651E0, + 0.33333333333333333333E0, + 0.31546487678572871854E0, + 0.30102999566398119521E0, + 0.28906482631788785926E0, + 0.27894294565112984319E0, + 0.27023815442731974129E0, + 0.26264953503719354797E0, + 0.25595802480981548938E0, + 0.25000000000000000000E0, + 0.24465054211822603038E0, + 0.23981246656813144473E0, + 0.23540891336663823644E0, + 0.23137821315975917426E0, + 0.22767024869695299798E0, + 0.22424382421757543947E0, + 0.22106472945750374614E0, + 0.21810429198553155922E0, + 0.21533827903669652533E0, + 0.21274605355336315360E0, + 0.21030991785715247903E0, + 0.20801459767650945759E0, + 0.20584683246043445730E0, + 0.20379504709050619003E0, + 0.20184908658209985071E0, + 0.19999999999999999999E0, + 0.19823986317056053160E0, + 0.19656163223282260771E0, + 0.19495902189378630772E0, + 0.19342640361727079343E0}; + +/* copy most important limbs of {op, n2} in {rp, n1} */ +/* if n1 > n2 put 0 in low limbs of {rp, n1} */ +#define MPN_COPY2(rp, n1, op, n2) \ + if ((n1) <= (n2)) \ + { \ + MPN_COPY ((rp), (op) + (n2) - (n1), (n1)); \ + } \ + else \ + { \ + MPN_COPY ((rp) + (n1) - (n2), (op), (n2)); \ + MPN_ZERO ((rp), (n1) - (n2)); \ + } -/* - Convert op to a string in base 'base' with 'n' digits and writes the - mantissa in 'str', the exponent in 'expptr'. - The result is rounded wrt 'rnd_mode'. +static double +_mpfr_ceil (double x) +{ + double y = (double) (long int) x; + return ((y < x) ? y + 1.0 : y); +} - For op = 3.1416 we get str = "31416" and expptr=1. - */ -char* -mpfr_get_str (char *str, mp_exp_t *expptr, int base, size_t n, - mpfr_srcptr op, mp_rnd_t rnd_mode) +/* this function computes an approximation of b^e in {a, n}, with exponent + stored in exp_r. The computed value is rounded towards zero (truncated). + It returns an integer f such that the final error is bounded by 2^f ulps, + that is: + a*2^exp_r <= b^e <= 2^exp_r (a + 2^f), + where a represents {a, n}, i.e. the integer + a[0] + a[1]*B + ... + a[n-1]*B^(n-1) where B=2^BITS_PER_MP_LIMB */ +static long +mpn_exp (mp_limb_t *a, mp_exp_t *exp_r, int b, mp_exp_t e, size_t n) { - int neg; + mp_limb_t *c, B; + mp_exp_t f, h; + int i; + unsigned long t; /* number of bits in e */ + unsigned long bits; + size_t n1; + int erreur; /* (number - 1) of loop a^2b inexact */ + /* erreur == t meens no error */ + int err_s_a2 = 0; + int err_s_ab = 0; /* number of error when shift A^2, AB */ + TMP_DECL(marker); - if (base < 2 || base > 36) - return NULL; + MPFR_ASSERTN(e > 0); + MPFR_ASSERTN((2 <= b) && (b <= 36)); - 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; - } + TMP_MARK(marker); - /* Do not use MPFR_PREC_MIN as this applies to base 2 only. Perhaps we - should allow n == 1 for directed rounding modes and odd bases. */ - if (n < 2) - return NULL; + /* initialization of a, b, f, h */ - 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; - } + /* normalize the base */ + B = (mp_limb_t) b; + count_leading_zeros (h, B); - neg = MPFR_SIGN(op) < 0; /* 0 if positive, 1 if negative */ + bits = BITS_PER_MP_LIMB - h; - if (MPFR_IS_INF(op)) - { - char *str0; + B = B << h; + h = - h; - 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 */ - } + /* allocate space for A and set it to B */ + c = (mp_limb_t*) TMP_ALLOC(2 * n * BYTES_PER_MP_LIMB); + a [n - 1] = B; + MPN_ZERO (a, n - 1); + /* initial exponent for A: invariant is A = {a, n} * 2^f */ + f = h - (n - 1) * BITS_PER_MP_LIMB; + + /* determine number of bits in e */ + count_leading_zeros (t, (mp_limb_t) e); + + t = BITS_PER_MP_LIMB - t; /* number of bits of exponent e */ - /* op is a floating-point number */ + erreur = t; - if (MPFR_IS_ZERO(op)) + MPN_ZERO (c, 2 * n); + + for (i = t - 2; i >= 0; i--) { - 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 */ + /* determine precision needed */ + bits = n * BITS_PER_MP_LIMB - mpn_scan1 (a, 0); + n1 = (n * BITS_PER_MP_LIMB - bits) / BITS_PER_MP_LIMB; + + /* square of A : {c+2n1, 2(n-n1)} = {a+n1, n-n1}^2 */ + mpn_sqr_n (c + 2 * n1, a + n1, n - n1); + + /* set {c+n, 2n1-n} to 0 : {c, n} = {a, n}^2*K^n */ + + f = 2 * f + n * BITS_PER_MP_LIMB; + if ((c[2*n - 1] & MPFR_LIMB_HIGHBIT) == 0) + { + /* shift A by one bit to the left */ + mpn_lshift (a, c + n, n, 1); + a[0] |= mpn_lshift (c + n - 1, c + n - 1, 1, 1); + f --; + if (erreur != t) + err_s_a2 ++; + } + else + MPN_COPY (a, c + n, n); + + if ((erreur == t) && (2 * n1 <= n) && + (mpn_scan1 (c + 2 * n1, 0) < (n - 2 * n1) * BITS_PER_MP_LIMB)) + erreur = i; + + if (e & (1 << i)) + { + /* multiply A by B */ + c[2 * n - 1] = mpn_mul_1 (c + n - 1, a, n, B); + f += h + BITS_PER_MP_LIMB; + if ((c[2 * n - 1] & MPFR_LIMB_HIGHBIT) == 0) + { /* shift A by one bit to the left */ + mpn_lshift (a, c + n, n, 1); + a[0] |= mpn_lshift (c + n - 1, c + n - 1, 1, 1); + f --; + } + else + { + MPN_COPY (a, c + n, n); + if (erreur != t) + err_s_ab ++; + } + if ((erreur == t) && (c[n - 1] != 0)) + erreur = i; + + } } + + TMP_FREE(marker); + + *exp_r = f; + + if (erreur == t) + erreur = -1; /* exact */ 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 there are p loops after the first inexact result, with + j shifts in a^2 and l shifts in a*b, then the final error is + at most 2^(p+ceil((j+1)/2)+l+1)*ulp(res). + This is bounded by 2^(5/2*t-1/2) where t is the number of bits of e. + */ + erreur = erreur + err_s_ab + err_s_a2 / 2 + 3; + if ((erreur - 1) >= ((n * BITS_PER_MP_LIMB - 1) / 2)) + erreur = n * BITS_PER_MP_LIMB; /* completely wrong: this is very + unlikely to happen since erreur is + at most 5/2*log_2(e), and + n * BITS_PER_MP_LIMB is at least + 3*log_2(e) */ + } - str_is_null = str == NULL; + return erreur; +} - 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 */ - } - else - pow2 = 0; - - /* first determines the exponent */ - e = MPFR_EXP(op); - - /* 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 + +/* Input: an approximation r*2^f of an real Y, with |r*2^f-Y| <= 2^(e+f). + Returns if possible in the string s the mantissa corresponding to + the integer nearest to Y, within the direction rnd, and returns the + the exponent in exp. + n is the number of limbs of r. + e represents the maximal error in the approximation of Y + (e < 0 iff the approximation is exact, i.e. r*2^f = Y). + b is the wanted base (2 <= b <= 36). + m is the number of wanted digits in the mantissa. + rnd is the rounding mode. + It is assumed that b^(m-1) <= Y < b^(m+1), thus the returned value + satisfies b^(m-1) <= rnd(Y) < b^(m+1). + + Return value: + - the direction of rounding (-1, 0, 1) if rounding is possible + - 2 otherwise (rounding is not possible) +*/ +static int +mpfr_get_str_aux (char *str, mp_exp_t *exp, mp_limb_t *r, mp_size_t n, + mp_exp_t f, long e, int b, size_t m, mp_rnd_t rnd) +{ + int dir; /* direction of the rounded result */ + mp_limb_t ret = 0; /* possible carry in addition */ + mp_size_t i0, j0; /* number of limbs and bits of Y */ + char *str1; /* string of m+2 characters */ + size_t size_s1; /* length of str1 */ + mp_rnd_t rnd1; + int i; + char c; + int exact = (e < 0); + TMP_DECL(marker); + + /* if f > 0, then the maximal error 2^(e+f) is larger than 2 so we can't + determine the integer Y */ + MPFR_ASSERTN(f <= 0); + /* if f is too small, then r*2^f is smaller than 1 */ + MPFR_ASSERTN(f > (-n * BITS_PER_MP_LIMB)); + + TMP_MARK(marker); + + /* R = 2^f sum r[i]K^(i) + r[i] = (r_(i,k-1)...r_(i,0))_2 + R = sum r(i,j)2^(j+ki+f) + the bits from R are referenced by pairs (i,j) */ + + /* check if is possible to round r with rnd mode + where |r*2^f-Y| <= 2^(e+f) + the exponent of R is: f + n*BITS_PER_MP_LIMB + we must have e + f == f + n*BITS_PER_MP_LIMB - err + err = n*BITS_PER_MP_LIMB - e + R contains exactly -f bits after the integer point: + to determine the nearest integer, we thus need a precision of + n * BITS_PER_MP_LIMB + f */ + + if (exact || mpfr_can_round_raw (r, n, (mp_size_t) 1, + n * BITS_PER_MP_LIMB - e, GMP_RNDN, rnd, n * BITS_PER_MP_LIMB + f)) + { + /* compute the nearest integer to R */ + + /* bit of weight 0 in R has position j0 in limb r[i0] */ + i0 = (-f) / BITS_PER_MP_LIMB; + j0 = (-f) % BITS_PER_MP_LIMB; + + ret = mpfr_round_raw_generic (r + i0, r, n * BITS_PER_MP_LIMB, 0, + n * BITS_PER_MP_LIMB + f, rnd, &dir, 0); + + /* warning: mpfr_round_raw_generic returns 2 or -2 in case of even + rounding */ + if (dir > 0) /* when dir = MPFR_EVEN_INEX */ + dir = 1; + else if (dir < 0) /* when dir = -MPFR_EVEN_INEX */ + dir = -1; + + if (ret) /* Y is a power of 2 */ + { + if (j0) + r[n - 1] = MPFR_LIMB_HIGHBIT >> (j0 - 1); + else /* j0=0, necessarily i0 >= 1 otherwise f=0 and r is exact */ + { + r[n - 1] = ret; + r[--i0] = 0; /* set to zero the new low limb */ + } + } + else /* shift r to the right by (-f) bits (i0 already done) */ + { + if (j0) + mpn_rshift (r + i0, r + i0, n - i0, j0); + } + + /* now the rounded value Y is in {r+i0, n-i0} */ + + /* convert r+i0 into base b */ + str1 = TMP_ALLOC (m + 3); /* need one extra character for mpn_get_str */ + size_s1 = mpn_get_str (str1, b, r + i0, n - i0); + + /* round str1 */ + *exp = size_s1 - m; /* number of superfluous characters */ + + /* if size_s1 = m + 2, necessarily we have b^(m+1) as result, + and the result will not change */ + + /* so we have to double-round only when size_s1 = m + 1 and + (i) the result is inexact + (ii) or the last digit is non-zero */ + if ((size_s1 == m + 1) && ((dir != 0) || (str1[size_s1 - 1] != 0))) { - double d; - - 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) + /* rounding mode */ + rnd1 = rnd; + + /* round to nearest case */ + if (rnd == GMP_RNDN) { - MPFR_ASSERTN(f < MP_EXP_T_MAX); - f++; + if (2 * str1[size_s1 - 1] == b) + { + if (dir == -1) + rnd1 = GMP_RNDU; + else if (dir == 0) /* exact: even rounding */ + { + if (exact) + rnd1 = ((str1[size_s1-2] & 1) == 0) + ? GMP_RNDD : GMP_RNDU; + else + goto cannot_round; + } + else + rnd1 = GMP_RNDD; + } + else if (2 * str1[size_s1 - 1] < b) + rnd1 = GMP_RNDD; + else + rnd1 = GMP_RNDU; } + + /* now rnd1 is either GMP_RNDD or GMP_RNDZ -> truncate + or GMP_RDNU -> round towards infinity */ + + /* round away from zero */ + if (rnd1 == GMP_RNDU) + { + if (str1[size_s1 - 1] != 0) + { + /* the carry cannot propagate to the whole string, since + Y = x*b^(m-g) < 2*b^m <= b^(m+1)-b + where x is the input float */ + for (c = 1, i = size_s1 - 2; c == 1; i--) + { + str1[i] ++; + if ((c = (str1[i] == b))) + str1[i] = 0; + } + } + dir = 1; + } + /* round toward zero (truncate) */ + else + dir = -1; } - mpfr_init (a); - mpfr_init (b); - mpz_init (bz); + /* copy str1 into str and convert to ASCII */ + for (i = 0; i < m; i++) + str[i] = num_to_text[(int) str1[i]]; + str[m] = 0; + } + /* mpfr_can_round_raw failed: rounding is not possible */ + else + { + cannot_round: + dir = 2; + } + + TMP_FREE(marker); - str0 = str; + return (dir); +} - 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; +/* prints the mantissa of x in the string s, and writes the corresponding + exponent in e. + x is rounded with direction rnd, m is the number of digits of the mantissa, + b is the given base (2 <= b <= 36). - d = (double) n / __mp_bases[base].chars_per_bit_exactly; - MPFR_ASSERTN(d <= MPFR_INTPREC_MAX - 1); - prec = (mp_prec_t) d + 1; - } + Return value: + if s=NULL, allocates a string to store the mantissa, with + m characters, plus a final '\0', plus a possible minus sign + (thus m+1 or m+2 characters). - 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); + Important: when you call this function with s=NULL, don't forget to free + the memory space allocated, with free(s, strlen(s)). +*/ +char* +mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd) +{ + int exact; /* exact result */ + mp_exp_t exp, g; + mp_exp_t prec, log_2prec; /* precision of the computation */ + long err; + mp_limb_t *a; + mp_exp_t exp_a; + mp_limb_t *result; + mp_limb_t *xp, *x1; + mp_limb_t *reste; + size_t nx, nx1; + size_t n, i; + char *s0; + int neg; - while (1) - { - mp_exp_unsigned_t p; - int div; + /* if exact = 1 then err is undefined */ + /* otherwise err is such that |x*b^(m-g)-a*2^exp_a| < 2^(err+exp_a) */ - 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; - } + TMP_DECL(marker); - 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); - } - } + /* is the base valid? */ + if (b < 2 || b > 36) + return NULL; - if (neg) - MPFR_CHANGE_SIGN(b); /* put b positive */ + if (m == 0) + { + m = (size_t) _mpfr_ceil (__mp_bases[b].chars_per_bit_exactly + * (double) MPFR_PREC(x)); + if (m < 2) + m = 2; + } - 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; + /* Do not use MPFR_PREC_MIN as this applies to base 2 only. Perhaps we + should allow n == 1 for directed rounding modes and odd bases. */ + MPFR_ASSERTN (m >= 2); - MPFR_ASSERTN(q <= MPFR_INTPREC_MAX - BITS_PER_MP_LIMB); - q += BITS_PER_MP_LIMB; - } + if (MPFR_IS_NAN(x)) + { + if (s == NULL) + s = (*__gmp_allocate_func) (4); + strcpy (s, "NaN"); + return s; + } - { - mp_rnd_t rnd = rnd_mode; + neg = MPFR_SIGN(x) < 0; /* 0 if positive, 1 if negative */ - if (neg) - switch (rnd_mode) - { - case GMP_RNDU: rnd = GMP_RNDZ; break; - case GMP_RNDD: rnd = GMP_RNDU; break; - } + if (MPFR_IS_INF(x)) + { + if (s == NULL) + s = (*__gmp_allocate_func) (neg + 4); + strcpy (s, (neg) ? "-Inf" : "Inf"); + return s; + } - 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)); + /* x is a floating-point number */ - if (len == n - 1) + if (MPFR_IS_ZERO(x)) + { + if (s == NULL) + s = (*__gmp_allocate_func) (neg + m + 1); + s0 = s; + if (neg) + *s++ = '-'; + memset (s, '0', m); + s[m] = '\0'; + *e = 0; /* a bit like frexp() in ISO C99 */ + return s0; /* strlen(s0) = neg + m */ + } + + /* si x < 0, on se ram`ene au cas x > 0 */ + if (neg) + { + switch (rnd) + { + case GMP_RNDU : rnd = GMP_RNDD; break; + case GMP_RNDD : rnd = GMP_RNDU; break; + } + } + + if (s == NULL) + s = (*__gmp_allocate_func) (neg + m + 1); + s0 = s; + if (neg) + *s++ = '-'; + + xp = MPFR_MANT(x); + + if (IS_POW2(b)) + { + int pow2; + mp_exp_t f, r; + mp_limb_t *x1; + mp_size_t nb; + int inexp; + + count_leading_zeros (pow2, (mp_limb_t) b); + pow2 = BITS_PER_MP_LIMB - pow2 - 1; /* base = 2^pow2 */ + + /* set MPFR_EXP(x) = f*pow2 + r, 1 <= r <= pow2 */ + f = (MPFR_EXP(x) - 1) / pow2; + r = MPFR_EXP(x) - f * pow2; + if (r <= 0) + { + f --; + r += pow2; + } + + /* the first digit will contain only r bits */ + prec = (m - 1) * pow2 + r; + n = (prec - 1) / BITS_PER_MP_LIMB + 1; + + x1 = TMP_ALLOC((n + 1) * sizeof (mp_limb_t)); + nb = n * BITS_PER_MP_LIMB - prec; + /* round xp to the precision prec, and put it into x1 + put the carry into x1[n] */ + if ((x1[n] = mpfr_round_raw_generic (x1, xp, MPFR_PREC(x), MPFR_ISNEG(x), + prec, rnd, &inexp, 0))) { - len++; - f--; - str[n-1]='0'; - str[n]='\0'; + /* overflow when rounding x: x1 = 2^prec */ + if (r == pow2) /* prec = m * pow2, + 2^prec will need (m+1) digits in base 2^pow2 */ + { + /* divide x1 by 2^pow2, and increase the exponent */ + mpn_rshift (x1, x1, n + 1, pow2); + f ++; + } + else /* 2^prec needs still m digits, but x1 may need n+1 limbs */ + n ++; } + + /* it remains to shift x1 by nb limbs to the right, since mpn_get_str + expects a right-normalized number */ + if (nb != 0) + { + mpn_rshift (x1, x1, n, nb); + /* the most significant word may be zero */ + if (x1[n - 1] == 0) + n --; + } + + mpn_get_str (s, b, x1, n); + for (i=0; i n) + exact = mpn_scan1 (xp, 0) >= (nx - n) * BITS_PER_MP_LIMB; + err = !exact; + MPN_COPY2 (a, n, xp, nx); + exp_a = MPFR_EXP(x) - n * BITS_PER_MP_LIMB; + } + else if ((mp_exp_t) m > g) /* we have to multiply x by b^exp */ + { + /* a2*2^exp_a = b^e */ + err = mpn_exp (a, &exp_a, b, exp, n); + /* here, the error on a is at most 2^err ulps */ + exact = (err == -1); + + /* x = x1*2^(n*BITS_PER_MP_LIMB) */ + x1 = (nx >= n) ? xp + nx - n : xp; + nx1 = (nx >= n) ? n : nx; /* nx1 = min(n, nx) */ + + /* test si exact */ + if (nx > n) + exact = (exact && + ((mpn_scan1 (xp, 0) >= (nx - n) * BITS_PER_MP_LIMB))); + + /* we loose one more bit in the multiplication, + except when err=0 where we loose two bits */ + err = (err <= 0) ? 2 : err + 1; + + /* result = a * x */ + result = TMP_ALLOC ((n + nx1) * sizeof (mp_limb_t)); + mpn_mul (result, a, n, x1, nx1); + exp_a += MPFR_EXP (x); + if (mpn_scan1 (result, 0) < (nx1 * BITS_PER_MP_LIMB)) + exact = 0; + + /* normalize a and truncate */ + if ((result[n + nx1 - 1] & MPFR_LIMB_HIGHBIT) == 0) + { + mpn_lshift (a, result + nx1, n , 1); + a[0] |= result[nx1 - 1] >> (BITS_PER_MP_LIMB - 1); + exp_a --; + } + else + MPN_COPY (a, result + nx1, n); + } + else + { + /* a2*2^exp_a = b^e */ + err = mpn_exp (a, &exp_a, b, exp, n); + exact = (err == -1); + + /* allocate memory for x1, result and reste */ + x1 = TMP_ALLOC (2 * n * sizeof (mp_limb_t)); + result = TMP_ALLOC ((n + 1) * sizeof (mp_limb_t)); + reste = TMP_ALLOC (n * sizeof (mp_limb_t)); + + /* initialize x1 = x */ + MPN_COPY2 (x1, 2 * n, xp, nx); + if ((exact) && (nx > 2 * n) && + (mpn_scan1 (xp, 0) < (nx - 2 * n) * BITS_PER_MP_LIMB)) + exact = 0; + + /* result = x / a */ + mpn_tdiv_qr (result, reste, 0, x1, 2 * n, a, n); + exp_a = MPFR_EXP (x) - exp_a - 2 * n * BITS_PER_MP_LIMB; + + /* test if division was exact */ + if (exact) + exact = mpn_popcount (reste, n) == 0; + + /* normalize the result and copy into a */ + if (result[n] == 1) + { + mpn_rshift (a, result, n, 1); + a[n - 1] |= MPFR_LIMB_HIGHBIT;; + exp_a ++; + } + else + MPN_COPY (a, result, n); + + err = (err == -1) ? 2 : err + 2; + } + + /* check if rounding is possible */ + if (exact) + err = -1; + if (mpfr_get_str_aux (s, e, a, n, exp_a, err, b, m, rnd) != 2) + break; + + /* increment the working precision */ + prec += log_2prec; + + TMP_FREE(marker); } + + *e += g; + + TMP_FREE(marker); + + return s0; + } -- cgit v1.2.1