needing help with: Floating point error

[jmonetti at uncu.edu.ar]

Hello List,

I'm experiencing a little problem with gmp, and hope
you can help me.
I tried to simplify my code and show you (below).  :)

The code declare some gmp variables, initializes,
and then start an iteration calculating the following
equations:


     d    = a + pow(0.05 , k );
     prod =(a*e)-(d*b);
     X    =(c*e-f*b)/prod;
     Y    =(a*f-d*c)/prod;

In each iteration k is incremented, so the pow
function returns values even smaller.
Finally it must operate with both big values,
and small values.

The problem I found that the program get a
Floating Point in iteration 29 (with my computer/
arquitectura/compiler/...).

Can you give me some help with it ?
Attached a simplified code and ouput.

Thanks in advanced.

Julio

CODE
----


#include <stdio.h>
#include <math.h>
#include <gmp.h>
#include <string.h>


  //declare
	mpf_t Y, X, prod;
	mpf_t a;
	mpf_t b;
	mpf_t e;
	mpf_t c;
	mpf_t d;
	mpf_t f;
	mpz_t aux0;
	mpf_t aux1;
	mpf_t aux2;
	mpf_t aux3;

	unsigned long int k;

 //initialize
	 mpf_init(Y);
	 mpf_init(X);
	 mpf_init(prod);
	 mpf_init( a );
	 mpf_init( b );
	 mpf_init( e );
	 mpf_init( c );
	 mpf_init( d );
	 mpf_init( f );
	 mpz_init( aux0 );
	 mpf_init( aux1 );
	 mpf_init( aux2 );
	 mpf_init( aux3 );


        //iterate
	for(k=0;k< atoi(argv[2]);k++) {

	 mpf_set_d(a, 0.123456789012345);
	 mpf_set_d(b, 0.987654321002340);
	 mpf_set_d(e, 0.987654321002340);
	 mpf_set_d(c, 0.1);
	 mpf_set_d(d, 0.0);
	 mpf_set_d(f, 0.93);
	 mpz_set_d(aux0, 0);
	 mpf_set_d(aux1, 0);
	 mpf_set_d(aux2, 0);
	 mpf_set_d(aux3, 0);

        //d= a + pow(0.05 , k );          //original equation
        mpf_set_d(aux2, 0.05);
        mpf_pow_ui(aux1 , aux2, k );     //this is critical!
        mpf_add(d, a, aux1) ;

       //prod=(a*e)-(d*b);
        mpf_mul(aux2,a,e);
        mpf_mul(aux3,d,b);
        mpf_sub( prod, aux2,aux3);

        //X=(c*e-f*b)/prod;
        mpf_mul(aux2,c,e);
        mpf_mul(aux3,f,b);
        mpf_sub(aux2,aux2,aux3);
        mpf_div( X, aux2, prod);

        //Y=(a*f-d*c)/prod;
        mpf_mul(aux2,a,f);
        mpf_mul(aux3,d,c);
        mpf_sub(aux2, aux2,aux3);
        mpf_div(Y,aux2,prod);

        printf("\niter %d ",k);
        gmp_printf(" %.Ff  \t %.Ff \t %.Ff",X,Y, prod );

        }  //for

	 mpf_clear( X );
	 mpf_clear( Y );
	 mpf_clear( prod );
	 mpf_clear( a );
	 mpf_clear( b );
	 mpf_clear( e );
	 mpf_clear( c );
	 mpf_clear( d );
	 mpf_clear( f );
	 mpz_clear( aux0 );
	 mpf_clear( aux1 );
	 mpf_clear( aux2 );
	 mpf_clear( aux3 );

  }
  return 0;
}
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: output.txt
URL: <http://gmplib.org/list-archives/gmp-discuss/attachments/20090626/b141a87c/attachment.txt>