mpq_get_d_nearest

[Zimmermann Paul]

       Hi,

the mpq_get_d function rounds towards zero (i.e., truncates).
In several applications, people usually prefer rounding to nearest.
Is it planned to provide a function (say mpq_get_d_nearest) for that?
I could contribute it, based on the code below (which does not deal with
subnormal numbers yet).

We could also have rounding towards -infinity and +infinity, which would be
useful for people doing interval arithmetic.

Paul

double
mpq_get_d_nearest (mpq_t q)
{
  mpz_ptr a = mpq_numref (q);
  mpz_ptr b = mpq_denref (q);
  size_t sa = mpz_sizeinbase (a, 2);
  size_t sb = mpz_sizeinbase (b, 2);
  size_t na, nb;
  mpz_t aa, bb;
  double d;
  
  /* easy case: |a|, |b| < 2^53, no overflow nor underflow can occur */
  if (sa <= 53 && sb <= 53)
    return mpz_get_d (a) / mpz_get_d (b);

  /* same if a = m*2^e with m representable on 53 bits, idem for b, but beware
     that both a and b do not give an overflow */
  na = sa - mpz_scan1 (a, 0);
  nb = sb - mpz_scan1 (b, 0);
  if (sa <= 1024 && na <= 53 && sb <= 1024 && nb <= 53)
    return mpz_get_d (a) / mpz_get_d (b);

  /* hard case */
  mpz_init (aa);
  mpz_init (bb);
  
  if (sa >= sb)
    {
      mpz_set (aa, a);
      mpz_mul_2exp (bb, b, sa - sb);
    }
  else
    {
      mpz_mul_2exp (aa, a, sb - sa);
      mpz_set (bb, b);
    }

  /* q = aa/bb*2^(sa-sb) */

  if (mpz_cmpabs (aa, bb) >= 0)
    {
      mpz_mul_2exp (bb, bb, 1);
      sa ++;
    }

  mpz_mul_2exp (aa, aa, 54);
  sb += 54;

  mpz_tdiv_qr (aa, bb, aa, bb);

  /* the quotient aa should have exactly 54 bits */

  if (mpz_tstbit (aa, 0) == 0)
    {
    }
  else if (mpz_cmp_ui (bb, 0) != 0)
    {
      if (mpz_sgn (aa) > 0)
        mpz_add_ui (aa, aa, 1);
      else
        mpz_sub_ui (aa, aa, 1);
    }
  else /* mid case: round to even */
    {
      if (mpz_tstbit (aa, 1) == 0)
        {
          if (mpz_sgn (aa) > 0)
            mpz_sub_ui (aa, aa, 1);
          else
            mpz_add_ui (aa, aa, 1);
        }
      else
        {
          if (mpz_sgn (aa) > 0)
            mpz_add_ui (aa, aa, 1);
          else
            mpz_sub_ui (aa, aa, 1);
        }
    }
  
  d = mpz_get_d (aa); /* exact */
  mpz_clear (aa);
  mpz_clear (bb);
  return ldexp (d, (long) sa - (long) sb);
}