mpz division to yield a double
Sam Rawlins
sam.rawlins at gmail.com
Mon Jun 18 20:27:56 CEST 2012
On Sun, Jun 17, 2012 at 5:21 AM, Zimmermann Paul
<Paul.Zimmermann at loria.fr>wrote:
> Sam,
>
> > If, however, one of the numbers "isinf," (example: [2]) then a second
> > method is called, which gets into the internals of the bignum
> > implementation to perform the division. The obvious solution would be to
> > initialize mpf_t's out of my divisor and dividend, but with (the current
> > version of) mini-gmp, I only have mpz_*, and mpn_* methods.
> >
> > [2] An example that falls into this "isinf" category would be something
> > like: 3*10^100 / 2*10^100, and I want the result as a double: 1.5.
>
> in that case, assume the numerator has n bits and the denominator has d
> bits,
> multiply the numerator by 2^max(0,d-n+54), divide the scaled numerator by
> the denominator (as mpz_t), then convert the integer quotient to double,
> and scale back dividing by 2^max(0,d-n+54) using the ldexp function for
> example.
>
> In your example both the numerator and denominator have 334 bits, you
> scale the
> numerator by 2^54, dividing by d gives 27021597764222976, dividing by 2^54
> gives 1.5.
>
Thanks Paul, and others who responded directly to me. I figured there was
something elementary I wasn't figuring out :/ . So I've implemented such
function now, and I'm posting it here, for completeness. It is exactly
Paul's method:
/* In mini-gmp, only some mpz_* and mpn_* functions are available. In
* order to carry out a division between two mpz_t's, either of which might
not
* fit into a double, and return a double, I implement the following
function:
*
* mpz_div_d(n, d):
* q = n * 2^(d_bits - n_bits + 54)
* q = n / d, using mpz_tdiv_q, so result is mpz_t
* q_d = mpz_get_d(q), converting q into a float
* return inf if isinf(q_d)
* q_d = q_d / 2^(d_bits - n_bits + 54)
* return q_d
*/
double mpz_div_d(mpz_t n, mpz_t d) {
mpz_t q;
size_t n_bits, d_bits;
double q_d;
int exponent;
mpz_init(q);
n_bits = mpz_sizeinbase(n, 2);
d_bits = mpz_sizeinbase(d, 2);
/* scale n out by 54 bits, plus the difference between their sizes */
mpz_mul_2exp(q, n, max(0, d_bits-n_bits + 54));
mpz_tdiv_q(q, q, d);
q_d = mpz_get_d(q);
/* mpz_get_d() may result in an infinity. From GMP manual:
*
* "If the exponent from the conversion is too big, the result is
system
* dependent. An infinity is returned where available. A hardware
overflow
* trap may or may not occur."
*
* So I'm just going to stop doing math here and return infinity if that's
* where things stand.
*/
if (isinf(q_d)) {
return q_d;
}
/* scale the quotient back */
exponent = max(0, -(d_bits-n_bits + 54));
q_d = ldexp(q_d, exponent); /* note: ldexp() can also return +/-HUGE_VAL
*/
return q_d;
}
I've tested it, not extensively.
> Paul Zimmermann
>
--
Sam Rawlins
More information about the gmp-discuss
mailing list