diff options
Diffstat (limited to 'contrib/libgmp/mpf/get_str.c')
| -rw-r--r-- | contrib/libgmp/mpf/get_str.c | 500 |
1 files changed, 500 insertions, 0 deletions
diff --git a/contrib/libgmp/mpf/get_str.c b/contrib/libgmp/mpf/get_str.c new file mode 100644 index 0000000..bfee18d --- /dev/null +++ b/contrib/libgmp/mpf/get_str.c @@ -0,0 +1,500 @@ +/* mpf_get_str (digit_ptr, exp, base, n_digits, a) -- Convert the floating + point number A to a base BASE number and store N_DIGITS raw digits at + DIGIT_PTR, and the base BASE exponent in the word pointed to by EXP. For + example, the number 3.1416 would be returned as "31416" in DIGIT_PTR and + 1 in EXP. + +Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP 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 +option) any later version. + +The GNU MP 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 +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP 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 "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +/* + New algorithm for converting fractions (951019): + 0. Call the fraction to convert F. + 1. Compute [exp * log(2^BITS_PER_MP_LIMB)/log(B)], i.e., + [exp * BITS_PER_MP_LIMB * __mp_bases[B].chars_per_bit_exactly]. Exp is + the number of limbs between the limb point and the most significant + non-zero limb. Call this result n. + 2. Compute B^n. + 3. F*B^n will now be just below 1, which can be converted easily. (Just + multiply by B repeatedly, and see the digits fall out as integers.) + We should interrupt the conversion process of F*B^n as soon as the number + of digits requested have been generated. + + New algorithm for converting integers (951019): + 0. Call the integer to convert I. + 1. Compute [exp * log(2^BITS_PER_MP_LIMB)/log(B)], i.e., + [exp BITS_PER_MP_LIMB * __mp_bases[B].chars_per_bit_exactly]. Exp is + the number of limbs between the limb point and the least significant + non-zero limb. Call this result n. + 2. Compute B^n. + 3. I/B^n can be converted easily. (Just divide by B repeatedly. In GMP, + this is best done by calling mpn_get_str.) + Note that converting I/B^n could yield more digits than requested. For + efficiency, the variable n above should be set larger in such cases, to + kill all undesired digits in the division in step 3. +*/ + +char * +#if __STDC__ +mpf_get_str (char *digit_ptr, mp_exp_t *exp, int base, size_t n_digits, mpf_srcptr u) +#else +mpf_get_str (digit_ptr, exp, base, n_digits, u) + char *digit_ptr; + mp_exp_t *exp; + int base; + size_t n_digits; + mpf_srcptr u; +#endif +{ + mp_size_t usize; + mp_exp_t uexp; + unsigned char *str; + size_t str_size; + char *num_to_text; + long i; /* should be size_t */ + mp_ptr rp; + mp_limb_t big_base; + size_t digits_computed_so_far; + int dig_per_u; + mp_srcptr up; + unsigned char *tstr; + mp_exp_t exp_in_base; + TMP_DECL (marker); + + TMP_MARK (marker); + usize = u->_mp_size; + uexp = u->_mp_exp; + + if (base >= 0) + { + if (base == 0) + base = 10; + num_to_text = "0123456789abcdefghijklmnopqrstuvwxyz"; + } + else + { + base = -base; + num_to_text = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; + } + + /* Don't compute more digits than U can accurately represent. + Also, if 0 digits were requested, give *exactly* as many digits + as can be accurately represented. */ + { + size_t max_digits = (((u->_mp_prec - 1) * BITS_PER_MP_LIMB) + * __mp_bases[base].chars_per_bit_exactly); + if (n_digits == 0 || n_digits > max_digits) + n_digits = max_digits; + } + + if (digit_ptr == 0) + { + /* We didn't get a string from the user. Allocate one (and return + a pointer to it) with space for `-' and terminating null. */ + digit_ptr = (char *) (*_mp_allocate_func) (n_digits + 2); + } + + if (usize == 0) + { + *exp = 0; + *digit_ptr = 0; + return digit_ptr; + } + + str = (unsigned char *) digit_ptr; + + /* Allocate temporary digit space. We can't put digits directly in the user + area, since we almost always generate more digits than requested. */ + tstr = (unsigned char *) TMP_ALLOC (n_digits + 3 * BITS_PER_MP_LIMB); + + if (usize < 0) + { + *digit_ptr = '-'; + str++; + usize = -usize; + } + + digits_computed_so_far = 0; + + if (uexp > usize) + { + /* The number has just an integral part. */ + mp_size_t rsize; + mp_size_t exp_in_limbs; + mp_size_t msize; + mp_ptr tp, xp, mp; + int cnt; + mp_limb_t cy; + mp_size_t start_str; + mp_size_t n_limbs; + + n_limbs = 2 + ((mp_size_t) (n_digits / __mp_bases[base].chars_per_bit_exactly) + / BITS_PER_MP_LIMB); + + /* Compute n such that [u/B^n] contains (somewhat) more than n_digits + digits. (We compute less than that only if that is an exact number, + i.e., exp is small enough.) */ + + exp_in_limbs = uexp; + + if (n_limbs >= exp_in_limbs) + { + /* The number is so small that we convert the entire number. */ + exp_in_base = 0; + rp = (mp_ptr) TMP_ALLOC (exp_in_limbs * BYTES_PER_MP_LIMB); + MPN_ZERO (rp, exp_in_limbs - usize); + MPN_COPY (rp + (exp_in_limbs - usize), u->_mp_d, usize); + rsize = exp_in_limbs; + } + else + { + exp_in_limbs -= n_limbs; + exp_in_base = (((exp_in_limbs * BITS_PER_MP_LIMB - 1)) + * __mp_bases[base].chars_per_bit_exactly); + + rsize = exp_in_limbs + 1; + rp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); + tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); + + rp[0] = base; + rsize = 1; + + count_leading_zeros (cnt, exp_in_base); + for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--) + { + mpn_mul_n (tp, rp, rp, rsize); + rsize = 2 * rsize; + rsize -= tp[rsize - 1] == 0; + xp = tp; tp = rp; rp = xp; + + if (((exp_in_base >> i) & 1) != 0) + { + cy = mpn_mul_1 (rp, rp, rsize, (mp_limb_t) base); + rp[rsize] = cy; + rsize += cy != 0; + } + } + + mp = u->_mp_d; + msize = usize; + + { + mp_ptr qp; + mp_limb_t qflag; + mp_size_t xtra; + if (msize < rsize) + { + mp_ptr tmp = (mp_ptr) TMP_ALLOC ((rsize+1)* BYTES_PER_MP_LIMB); + MPN_ZERO (tmp, rsize - msize); + MPN_COPY (tmp + rsize - msize, mp, msize); + mp = tmp; + msize = rsize; + } + else + { + mp_ptr tmp = (mp_ptr) TMP_ALLOC ((msize+1)* BYTES_PER_MP_LIMB); + MPN_COPY (tmp, mp, msize); + mp = tmp; + } + count_leading_zeros (cnt, rp[rsize - 1]); + cy = 0; + if (cnt != 0) + { + mpn_lshift (rp, rp, rsize, cnt); + cy = mpn_lshift (mp, mp, msize, cnt); + if (cy) + mp[msize++] = cy; + } + + { + mp_size_t qsize = n_limbs + (cy != 0); + qp = (mp_ptr) TMP_ALLOC ((qsize + 1) * BYTES_PER_MP_LIMB); + xtra = qsize - (msize - rsize); + qflag = mpn_divrem (qp, xtra, mp, msize, rp, rsize); + qp[qsize] = qflag; + rsize = qsize + qflag; + rp = qp; + } + } + } + + str_size = mpn_get_str (tstr, base, rp, rsize); + + if (str_size > n_digits + 3 * BITS_PER_MP_LIMB) + abort (); + + start_str = 0; + while (tstr[start_str] == 0) + start_str++; + + for (i = start_str; i < str_size; i++) + { + tstr[digits_computed_so_far++] = tstr[i]; + if (digits_computed_so_far > n_digits) + break; + } + exp_in_base = exp_in_base + str_size - start_str; + goto finish_up; + } + + exp_in_base = 0; + + if (uexp > 0) + { + /* The number has an integral part, convert that first. + If there is a fractional part too, it will be handled later. */ + mp_size_t start_str; + + rp = (mp_ptr) TMP_ALLOC (uexp * BYTES_PER_MP_LIMB); + up = u->_mp_d + usize - uexp; + MPN_COPY (rp, up, uexp); + + str_size = mpn_get_str (tstr, base, rp, uexp); + + start_str = 0; + while (tstr[start_str] == 0) + start_str++; + + for (i = start_str; i < str_size; i++) + { + tstr[digits_computed_so_far++] = tstr[i]; + if (digits_computed_so_far > n_digits) + { + exp_in_base = str_size - start_str; + goto finish_up; + } + } + + exp_in_base = str_size - start_str; + /* Modify somewhat and fall out to convert fraction... */ + usize -= uexp; + uexp = 0; + } + + if (usize <= 0) + goto finish_up; + + /* Convert the fraction. */ + { + mp_size_t rsize, msize; + mp_ptr rp, tp, xp, mp; + int cnt; + mp_limb_t cy; + mp_exp_t nexp; + + big_base = __mp_bases[base].big_base; + dig_per_u = __mp_bases[base].chars_per_limb; + + /* Hack for correctly (although not efficiently) converting to bases that + are powers of 2. If we deem it important, we could handle powers of 2 + by shifting and masking (just like mpn_get_str). */ + if (big_base < 10) /* logarithm of base when power of two */ + { + int logbase = big_base; + if (dig_per_u * logbase == BITS_PER_MP_LIMB) + dig_per_u--; + big_base = (mp_limb_t) 1 << (dig_per_u * logbase); + /* fall out to general code... */ + } + +#if 0 + if (0 && uexp == 0) + { + rp = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB); + up = u->_mp_d; + MPN_COPY (rp, up, usize); + rsize = usize; + nexp = 0; + } + else + {} +#endif + uexp = -uexp; + if (u->_mp_d[usize - 1] == 0) + cnt = 0; + else + count_leading_zeros (cnt, u->_mp_d[usize - 1]); + + nexp = ((uexp * BITS_PER_MP_LIMB) + cnt) + * __mp_bases[base].chars_per_bit_exactly; + + if (nexp == 0) + { + rp = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB); + up = u->_mp_d; + MPN_COPY (rp, up, usize); + rsize = usize; + } + else + { + rsize = uexp + 2; + rp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); + tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); + + rp[0] = base; + rsize = 1; + + count_leading_zeros (cnt, nexp); + for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--) + { + mpn_mul_n (tp, rp, rp, rsize); + rsize = 2 * rsize; + rsize -= tp[rsize - 1] == 0; + xp = tp; tp = rp; rp = xp; + + if (((nexp >> i) & 1) != 0) + { + cy = mpn_mul_1 (rp, rp, rsize, (mp_limb_t) base); + rp[rsize] = cy; + rsize += cy != 0; + } + } + + /* Did our multiplier (base^nexp) cancel with uexp? */ +#if 0 + if (uexp != rsize) + { + do + { + cy = mpn_mul_1 (rp, rp, rsize, big_base); + nexp += dig_per_u; + } + while (cy == 0); + rp[rsize++] = cy; + } +#endif + mp = u->_mp_d; + msize = usize; + + tp = (mp_ptr) TMP_ALLOC ((rsize + msize) * BYTES_PER_MP_LIMB); + if (rsize > msize) + cy = mpn_mul (tp, rp, rsize, mp, msize); + else + cy = mpn_mul (tp, mp, msize, rp, rsize); + rsize += msize; + rsize -= cy == 0; + rp = tp; + + /* If we already output digits (for an integral part) pad + leading zeros. */ + if (digits_computed_so_far != 0) + for (i = 0; i < nexp; i++) + tstr[digits_computed_so_far++] = 0; + } + + while (digits_computed_so_far <= n_digits) + { + /* For speed: skip trailing zeroes. */ + if (rp[0] == 0) + { + rp++; + rsize--; + if (rsize == 0) + { + n_digits = digits_computed_so_far; + break; + } + } + + cy = mpn_mul_1 (rp, rp, rsize, big_base); + if (digits_computed_so_far == 0 && cy == 0) + { + abort (); + nexp += dig_per_u; + continue; + } + /* Convert N1 from BIG_BASE to a string of digits in BASE + using single precision operations. */ + { + unsigned char *s = tstr + digits_computed_so_far + dig_per_u; + for (i = dig_per_u - 1; i >= 0; i--) + { + *--s = cy % base; + cy /= base; + } + } + digits_computed_so_far += dig_per_u; + } + if (exp_in_base == 0) + exp_in_base = -nexp; + } + + finish_up: + + /* We can have at most one leading 0. Remove it. */ + if (tstr[0] == 0) + { + tstr++; + digits_computed_so_far--; + exp_in_base--; + } + + /* We should normally have computed too many digits. Round the result + at the point indicated by n_digits. */ + if (digits_computed_so_far > n_digits) + { + /* Round the result. */ + if (tstr[n_digits] * 2 >= base) + { + digits_computed_so_far = n_digits; + for (i = n_digits - 1; i >= 0; i--) + { + unsigned int x; + x = ++(tstr[i]); + if (x < base) + goto rounded_ok; + digits_computed_so_far--; + } + tstr[0] = 1; + digits_computed_so_far = 1; + exp_in_base++; + rounded_ok:; + } + } + + /* We might have fewer digits than requested as a result of rounding above, + (i.e. 0.999999 => 1.0) or because we have a number that simply doesn't + need many digits in this base (i.e., 0.125 in base 10). */ + if (n_digits > digits_computed_so_far) + n_digits = digits_computed_so_far; + + /* Remove trailing 0. There can be many zeros. */ + while (n_digits != 0 && tstr[n_digits - 1] == 0) + n_digits--; + + /* Translate to ascii and null-terminate. */ + for (i = 0; i < n_digits; i++) + *str++ = num_to_text[tstr[i]]; + *str = 0; + *exp = exp_in_base; + TMP_FREE (marker); + return digit_ptr; +} + +#if COPY_THIS_TO_OTHER_PLACES + /* Use this expression in lots of places in the library instead of the + count_leading_zeros+expression that is used currently. This expression + is much more accurate and will save odles of memory. */ + rsize = ((mp_size_t) (exp_in_base / __mp_bases[base].chars_per_bit_exactly) + + BITS_PER_MP_LIMB) / BITS_PER_MP_LIMB; +#endif |
