Re: Problem with Floating-point functions

[Ronald Bruck]

Armando Garcia B. wrote:

> I have builded a little program using GNU MP in a Linux RedHat 6.2
> system runing in a Pentium II 400 MHz machine (Dell Precision
> Worksatation 410).
>
> This simple program just asign a decimal number to a mpf_t variable and
> then the variable content is converted to a string. Finally the string
> is printed to the standard output.
>
> The program is:
>
> #include <stdio.h>
> #include <gmp.h>
>
> main()
> {
>   mpf_t x;
>   mp_exp_t exp;
>
>   mpf_init_set_d( x, 0.001 );   /* x  = 0.001 */
>
>   mpf_get_str( NULL, &exp, 10, 0, x );
>   printf( "x(mpf_get_str)=%s(%ld)\n", mpf_get_str( NULL, &exp, 10, 0, x
> ), exp );
>
>   exit( 0 );
> }


Let me try this one, as one novice to another :-)

This is a good lesson why one should NOT use the set_d functions when 
you need precision.  The string 0.001 is being converted to a double, 
which in IEEE 754 arithmetic has only 53 bits' significance (52 
"genuine" bits, and one implicit bit, the leading "1"), which is then 
being converted to a 64-bit mantissa.  Result:  You lose precision.

Use "mpf_init_set_str" instead.  Even better:  unless you're wedded 
to decimal for emotional reasons, use 1/1024 instead of 1/1000.  I do 
this a lot; when I need to set a tolerance of about 1000 * machine 
epsilon, I do something like

  mpf_init2(tol, precision);
  mpf_set_ui(tol, 1);
  mpf_div_2exp(tol, tol, precision-10);

That's 1024 * machine epsilon, but it's precise in binary.  And if I 
want 1024 * sqrt(machine_epsilon) I do

  mpf_div_2exp(tol, tol, precision/2-10)

instead.


About the earlier problem I reported:  yup, I was overwriting an 
array bound.  Ain't efence wonderful?

--Ron Bruck