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); } _______________________________________________ gmp-devel mailing list gmp-devel@gmplib.org http://gmplib.org/mailman/listinfo/gmp-devel