mpfr-3.1.5.tar.gz

/* mpfr_get_str -- output a floating-point number to a string Copyright 1999-2016 Free Software Foundation, Inc. Contributed by the AriC and Caramba projects, INRIA. This file is part of the GNU MPFR Library. The GNU MPFR Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. The GNU 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 Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ #define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" static int mpfr_get_str_aux (char *const, mpfr_exp_t *const, mp_limb_t *const, mp_size_t, mpfr_exp_t, long, int, size_t, mpfr_rnd_t); /* The implicit \0 is useless, but we do not write num_to_text[62] otherwise g++ complains. */ static const char num_to_text36[] = "0123456789abcdefghijklmnopqrstuvwxyz"; static const char num_to_text62[] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ" "abcdefghijklmnopqrstuvwxyz"; /* 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)); \ } #define MPFR_ROUND_FAILED 3 /* Input: an approximation r*2^f of a 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 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 <= 62). 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). Rounding may fail for two reasons: - the error is too large to determine the integer N nearest to Y - either the number of digits of N in base b is too large (m+1), N=2*N1+(b/2) and the rounding mode is to nearest. This can only happen when b is even. Return value: - the direction of rounding (-1, 0, 1) if rounding is possible - -MPFR_ROUND_FAILED if rounding not possible because m+1 digits - MPFR_ROUND_FAILED otherwise (too large error) */ static int mpfr_get_str_aux (char *const str, mpfr_exp_t *const exp, mp_limb_t *const r, mp_size_t n, mpfr_exp_t f, long e, int b, size_t m, mpfr_rnd_t rnd) { const char *num_to_text; 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 */ unsigned char *str1; /* string of m+2 characters */ size_t size_s1; /* length of str1 */ mpfr_rnd_t rnd1; size_t i; int exact = (e < 0); MPFR_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 * GMP_NUMB_BITS)); MPFR_TMP_MARK(marker); num_to_text = b < 37 ? num_to_text36 : num_to_text62; /* 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*GMP_NUMB_BITS we must have e + f == f + n*GMP_NUMB_BITS - err err = n*GMP_NUMB_BITS - e R contains exactly -f bits after the integer point: to determine the nearest integer, we thus need a precision of n * GMP_NUMB_BITS + f */ if (exact || mpfr_can_round_raw (r, n, (mp_size_t) 1, n * GMP_NUMB_BITS - e, MPFR_RNDN, rnd, n * GMP_NUMB_BITS + f)) { /* compute the nearest integer to R */ /* bit of weight 0 in R has position j0 in limb r[i0] */ i0 = (-f) / GMP_NUMB_BITS; j0 = (-f) % GMP_NUMB_BITS; ret = mpfr_round_raw (r + i0, r, n * GMP_NUMB_BITS, 0, n * GMP_NUMB_BITS + f, rnd, &dir); MPFR_ASSERTD(dir != MPFR_ROUND_FAILED); /* warning: mpfr_round_raw_generic returns MPFR_EVEN_INEX (2) or -MPFR_EVEN_INEX (-2) in case of even rounding */ 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 = (unsigned char*) MPFR_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 */ MPFR_ASSERTN(size_s1 >= m); *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))) { /* rounding mode */ rnd1 = rnd; /* round to nearest case */ if (rnd == MPFR_RNDN) { if (2 * str1[size_s1 - 1] == b) { if (dir == 0 && exact) /* exact: even rounding */ { rnd1 = ((str1[size_s1 - 2] & 1) == 0) ? MPFR_RNDD : MPFR_RNDU; } else { /* otherwise we cannot round correctly: for example if b=10, we might have a mantissa of xxxxxxx5.00000000 which can be rounded to nearest to 8 digits but not to 7 */ dir = -MPFR_ROUND_FAILED; MPFR_ASSERTD(dir != MPFR_EVEN_INEX); goto free_and_return; } } else if (2 * str1[size_s1 - 1] < b) rnd1 = MPFR_RNDD; else rnd1 = MPFR_RNDU; } /* now rnd1 is either MPFR_RNDD or MPFR_RNDZ -> truncate, or MPFR_RNDU or MPFR_RNDA -> round toward infinity */ /* round away from zero */ if (rnd1 == MPFR_RNDU || rnd1 == MPFR_RNDA) { 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 */ MPFR_ASSERTN(size_s1 >= 2); i = size_s1 - 2; while (str1[i] == b - 1) { MPFR_ASSERTD(i > 0);