# mpq_get_d_nearest

Zimmermann Paul Paul.Zimmermann at inria.fr
Fri Apr 12 11:13:49 CEST 2013

```       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)
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
}
else
{
if (mpz_sgn (aa) > 0)
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);
}
```