mpz division to yield a double

[Sam Rawlins]

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