mpf_sub bug (essentially same as reldiff bug)

[abbott at module.dima.unige.it]

Hi,

SUMMARY:
  mpf_sub wrongly gives zero as the difference of 1 and 1+eps
  where eps is 2^(-129).

PLATFORM:
Linux point 2.4.22 #1 Thu Sep 11 12:43:28 UTC 2003 i686 unknown
[actually Athlon XP2200 running slackware]


COMPILER:
Reading specs from /usr/lib/gcc-lib/i386-slackware-linux/3.2.2/specs
Configured with: ../gcc-3.2.2/configure --prefix=/usr --enable-shared
--enable-threads=posix --enable-__cxa_atexit --disable-checking
--with-gnu-ld --verbose --target=i386-slackware-linux
--host=i386-slackware-linux
Thread model: posix
gcc version 3.2.2


PROGRAM EXHIBITING BUG: (C++ source)
#include <gmp.h>
#include <iostream>

using namespace std;

int main()
{
  mpf_t a, b, rd;
  mpf_init2(a, 256);        // ask for 256 bits' precision
  mpf_init2(b, 256);        // again
  mpf_init2(rd, 8);         // I don't need much precision here

  mpf_set_ui(a, 1);         // a = 1
  mpf_set_ui(b, 1);
  mpf_div_2exp(b, b, 129);
  mpf_add(b, b, a);         // b = 1 + 2^(-129)
  {
    // print out mantissa just to prove it is not equal to 1
    char buffer[99];
    long exp;
    mpf_get_str(buffer, &exp, 10, 60, b);
    cout << "mantissa=" << buffer << endl;
  }

  mpf_sub(rd, a, b);   // should be about 2^(-129), but gives zero
  double mant;
  long exp;
  mant = mpf_get_d_2exp(&exp, rd);
  cout << "diff is approx " << mant << " times 2^(" << exp << ")" << endl;

  return 0;
}