[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

inexact->exact bignum rounding

From: Kevin Ryde
Subject: inexact->exact bignum rounding
Date: Thu, 02 Oct 2003 10:40:53 +1000
User-agent: Gnus/5.1003 (Gnus v5.10.3) Emacs/21.3 (gnu/linux)

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

   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.  */

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);
  /* GMP 4.2 and up */
  result = mpz_get_d (SCM_I_BIG_MPZ (b));

  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)),

  scm_remember_upto_here_1 (b);
  return result;

reply via email to

[Prev in Thread] Current Thread [Next in Thread]