BUG in gmp-4.2.1: mpf_get_str on large arg

abbott abbott at dima.unige.it
Thu Jan 4 22:53:46 CET 2007


I'd like to report a bug in gmp-4.2.1; I imagine the bug exhibits itself
on any 32-bit platform.

There is a problem converting large mpf_t values into strings.  The bug can be
triggered by running the following program:
/*******************************************************/
#include <stdio.h>
#include "gmp.h"

/* Program should print out 0.1e4900000001 */
int main()
{
  mpf_t overflow;
  mpf_init2(overflow, 64);
  mpf_set_ui(overflow, 10);
  mpf_pow_ui(overflow, overflow, 70000);
  mpf_pow_ui(overflow, overflow, 70000);
  mpf_out_str(stdout, 10, 20, overflow); /* gives SEGV */
  return 0;
}
/*******************************************************/


The cause appears to be around line 223 in the file mpf/get_str.c:
      e = (unsigned long) n_less_limbs_needed * (GMP_NUMB_BITS * mp_bases[base].chars_per_bit_exactly);

The cast to unsigned long should only be attempted if the value actually
fits into an unsigned long; in this case, 70000*70000 > 2^32.  About 10
lines later the program tries to TMP_ALLOC_LIMBS almost 2Gbytes of space,
and the call to mpn_tdiv_qr five lines later crashes.

Further comments: one possible fix could be to be determine the upper bits
of the exponent (and then divide accordingly) then deal with the lower bits
of the exponent.  Things might get tricky on a 64-bit platform if the double
value computed has less than 64 bits in the mantissa (e.g. 53 bits); still,
I believe that dealing first with the upper part of the exponent and then
the lower part should allow everything to work correctly (even allowing for
the limited precision in chars_per_bit_exactly).


In contrast if I call mpf_get_d_2exp instead of mpf_out_str then the
exponent wraps silently.

John Abbott.
PS in case you're curious, I wanted to compute factorial(10^9) in an mpf_t.
PPS Happy New Year!


More information about the gmp-bugs mailing list