[Gcl-devel] Bug in mpz_get_d

[Camm Maguire]

Greetings, and thanks for your reply!

Enrique Perez-Terron <enrio at online.no> writes:

> 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
                    ^
                    This is the problem.  Even rounding toward 0 at
                    all times this should be ....673 from the digit
                    string above.  This is presuming that rounding
                    toward zero in any base should not produce a number 
                    lower by more than one least significant digit in
                    any other base.  (true?)

>
> 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   
>

(in lisp (1- a) = a - 1, so the sign issue is correct here).

>
> 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.
>

This was my thought -- do you agree that the mpz_get_d algorithm should
produce at least 73?

I have a temporary ugly workaround which calls mpz_get_str, and sscanf
for big integers, which is working for now, but clearly less than
optimal. 


> Gcl already links with mpfr, no?
>

Alas, no.  We just use mpz.  For now at least.

Many thanks for the below!!!

Take care,

> 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));
>
> }
>
>
>
>
>
>

-- 
Camm Maguire			     		    camm at maguirefamily.org
==========================================================================
"The earth is but one country, and mankind its citizens."  --  Baha'u'llah