[Gcl-devel] Bug in mpz_get_d

[Enrique Perez-Terron]

On Fri, 2009-12-11 at 15:12 -0500, Camm Maguire wrote:
> Greetings!  Workaround?
> 
> >(setq a (numerator **))
> 
> 10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069375
> 
> >(float a)
> 
> 1.0715086071862672E301
> 
> >(float (1+ a))
> 
> 1.0715086071862673E301
> 
> >(float (1- a))
> 
> 1.0715086071862672E301

Testing GMP and MPFR directly
$ ./test-gmp 
Using GMP and mpf_get_d:
================
num = 10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069375
1.071508607186267202E+301    

1 + num = 10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376
1.0715086071862673209E+301   

1 - num = -10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069374
-1.071508607186267202E+301   


Using MPFR:
==================

Num = 1.0715086071862673209E+301   

1 + num = 1.0715086071862673209E+301   

1 - num = -1.0715086071862673209E+301  

I don't know why float(1-a) turned positive. 

It is easy to miss that there is one correct and two wrong digit
strings, not the other way around.
Numbers ending in 7202 are the wrong ones, correct ending is 73209.

Gcl already links with mpfr, no?

Below is the C code. Compile with -lmpfr -lgmp.
-Enrique
------------------------8<-----------------------------
#include <stdio.h>
#include <gmp.h>
#include <mpfr.h>

char *mynum = "10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069375";

main ()
{
  printf("Using GMP and mpf_get_d:\n================\n");
  mpz_t num;
  mpz_init_set_str(num, mynum, 10);
  printf("num = ");
  mpz_out_str(stdout, 10, num);
  putchar('\n'); 

  mpf_t myfloat;
  mpf_init(myfloat);
  mpf_set_z(myfloat, num);
  printf("%-29.20lG\n\n", mpf_get_d(myfloat));

  mpz_t other;
  mpz_init_set_str(other, "1", 10);
  mpz_add(other, other, num);
  printf("1 + num = ");
  mpz_out_str(stdout, 10, other);
  putchar('\n'); 

  mpf_set_z(myfloat, other);
  printf("%-29.20lG\n\n", mpf_get_d(myfloat));

  mpz_set_ui(other, 1);
  mpz_sub(other, other, num);
  printf("1 - num = ");
  mpz_out_str(stdout, 10, other);
  putchar('\n'); 

  mpf_set_z(myfloat, other);
  printf("%-29.20lG\n", mpf_get_d(myfloat));

  printf("\n\nUsing MPFR:\n==================\n\n");
  mpfr_t mpfr;
  mpfr_init(mpfr);
  printf("Num = ");
  mpfr_set_z(mpfr, num, GMP_RNDN);
  printf("%-29.20lG\n\n", mpfr_get_d(mpfr, GMP_RNDN));

  mpz_set_ui(other, 1);
  mpz_add(other, other, num);
  printf("1 + num = ");

  mpfr_set_z(mpfr, other, GMP_RNDN);
  printf("%-29.20lG\n\n", mpfr_get_d(mpfr, GMP_RNDN));

  mpz_set_ui(other, 1);
  mpz_sub(other, other, num);
  printf("1 - num = ");
  mpfr_set_z(mpfr, other, GMP_RNDN);
  printf("%-29.20lG\n", mpfr_get_d(mpfr, GMP_RNDN));

}