* inexact->exact bignum rounding
@ 2003-10-02 0:40 Kevin Ryde
2003-10-07 16:05 ` Marius Vollmer
0 siblings, 1 reply; 3+ messages in thread
From: Kevin Ryde @ 2003-10-02 0:40 UTC (permalink / raw)
I'd like to propose a new style scm_i_big2dbl, below (not really
tested yet).
This is of course used by exact->inexact, and also a few other places
where it makes sense to do the same as exact->inexact (in particular
make-rectangular for instance).
/* scm_i_big2dbl() rounds to the closest representable double, in accordance
with R5RS exact->inexact.
In GMP prior to 4.2, mpz_get_d rounding was unspecified. It followed the
hardware rounding mode, but on the absolute value actually. We get a 4.2
style truncate to zero by applying mpz_get_d to just the high
DBL_MANT_DIG bits (no rounding). This is not fast, but it's only for old
GMP.
In GMP 4.2 and up, mpz_get_d rounds towards zero. It basically takes the
high DBL_MANT_DIG bits, and we adjust to our desired round-to-nearest by
looking at the next bit below those.
Note that bignums exactly half way between representable doubles are
rounded to the next higher absolute value. This seems like an adequate
interpretation of "numerically closest", and it's easier to do than a
full "nearest-even" style. */
double
scm_i_big2dbl (SCM b)
{
double result;
size_t bits;
#if __GNU_MP_VERSION < 4 \
|| (__GNU_MP_VERSION == 4 && __GNU_MP_VERSION_MINOR < 2)
/* GMP prior to 4.2 */
mpz_t tmp;
mpz_init (tmp);
mpz_tdiv_q_2exp (tmp, SCM_I_BIG_MPZ (b), bits-DBL_MANT_DIG);
result = ldexp (mpz_get_d (tmp), bits-DBL_MANT_DIG);
mpz_clear (tmp);
#else
/* GMP 4.2 and up */
result = mpz_get_d (SCM_I_BIG_MPZ (b));
#endif
bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2);
if (bits > DBL_MANT_DIG)
if (mpz_tstbit (SCM_I_BIG_MPZ (b), bits-DBL_MANT_DIG-1))
result += ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b)),
bits-DBL_MANT_DIG-1);
scm_remember_upto_here_1 (b);
return result;
}
_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel
^ permalink raw reply [flat|nested] 3+ messages in thread
end of thread, other threads:[~2003-10-19 0:36 UTC | newest]
Thread overview: 3+ messages (download: mbox.gz follow: Atom feed
-- links below jump to the message on this page --
2003-10-02 0:40 inexact->exact bignum rounding Kevin Ryde
2003-10-07 16:05 ` Marius Vollmer
2003-10-19 0:36 ` Kevin Ryde
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).