BUG in gmp-4.2.1: mpf_get_str on large arg

[Torbjorn Granlund]

abbott <abbott at dima.unige.it> writes:

  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;
  }
  /*******************************************************/
  
GMP does not handle overflows well, not here and not elsewhere.

You can also silently get exponent overflow during arithmetic.  I am
not sure what the right approach would be for fixing this.  Overflow
is irrelevant for most applications, and checking for it in GMP when
doing exponent arithmetic would be costly.  And how should GMP handle
it?  Abort?  I suppose a sticky inf would be cleanest, but would also
be costliest.

Does MPFR handle this better?  It wouldn't surprise me if it did.

Wrt the overflow condition you ran into.  It is systematic overflow in
the get_str function.  I played with making the exponent calculation
like this,

      n_less_limbs_needed = ue - n_limbs_needed;
      ed = n_less_limbs_needed * (GMP_NUMB_BITS * mp_bases[base].chars_per_bit_exactly);
      eh = ed / (1.0 + (~CNST_LIMB(0)));
      el = ed - eh * (1.0 + (~CNST_LIMB(0)));

i.e., using two exponent words, eh and el.  (The 1.0 + (~CNST_LIMB(0))
expression computes the base of limbs as an floating point number.)
Then I fixed mpn_pow_1_highpart to take two exponent words.

We then run into that we cannot represent the exponent in the output
base ("exp_in_base") and have to make that into two words too...
Unfortunately, mpf_get_str is a documented interface, so changing it
cannot be done.  An mpf_get_str_abbott will be needed, with a two word
exponent.

I suppose my approach will also not work for huge bases on 64-bit
machines, which you no doubt would report next if I went with this
solution.

I leave this for now, to get it right and complete will take more time
than what I have at this point.  I don't consider the bug critical,
given the general lack of overflow checking in mpf.  If this is a
critical bug, then all of mpf is a critical bug.  :-)
  
  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).
  
Perhaps that will work, I haven't thought it through.  An alternative
might be to store the logarithms (mp_bases[base].chars_per_bit_exactly)
as fractions of limbs, f/B, where B is the limbs base.

I'll wait for your nice fix to this problem.  :-)

Please check that the other conversion functions work as expected too,
e.g., mpf_set_str.  From looking at its source, I feel pretty confident
it needs fixing too.  It reads the exponent with a dumb old strtol call.
Of course, we need our own two word strtol2...  mpf_set_str's
mpn_pow_1_highpart will need a two word exponent.

Please tell me if you want my half-done hack.

-- 
Torbjörn