mpq_get_d_nearest

[Marc Nieper-Wißkirchen]

Hi,

There was a discussion on the mailing list roughly ten years ago about
adding a function mpq_get_d_nearest, which would be more useful than
mpq_get_d in many use cases:

https://gmplib.org/list-archives/gmp-devel/2013-April/003223.html

What happened to this suggestion?  Or did it come to nothing?

I attached an attempt to handle subnormal numbers in Paul's code.

Thanks,

Marc

**

#if FLT_RADIX != 2
#error "unsupported float radix"
#endif

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^DBL_MANT_DIG, no overflow nor underflow can
occur */
  if (sa <= DBL_MANT_DIG && sb <= DBL_MANT_DIG)
    return mpz_get_d (a) / mpz_get_d (b);

  /* same if a = m*2^e with m representable on DBL_MANT_DIG 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 <= DBL_MAX_EXP && na <= DBL_MANT_DIG && sb <= DBL_MAX_EXP && nb <=
DBL_MANT_DIG)
    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 ++;
    }

  mp_bitcnt_t bits = DBL_MANT_DIG + 1;

  /* handle subnormal numbers */
  mp_bitcnt_t e = sb + DBL_MIN_EXP + 1;
  if (sa < e)
    {
      e -= sa;
      if (e >= bits)
        bits = 1;
      else
        bits -= e;
    }

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

  mpz_tdiv_qr (aa, bb, aa, bb);

  /* the quotient aa should have exactly BITS 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);
}

**