talexander
/
freebsd_amp_hwpstate
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
freebsd_amp_hwpstate/contrib/libgmp/mpf/get_str.c

501 lines
13 KiB
C
Raw Blame History

/* 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
Reference in New Issue View Git Blame Copy Permalink