needing help with: Floating point error

[jmonetti at uncu.edu.ar]

Hello Torbjörn,

Thanks for responding.

>Unfortunately, your code is simplified to the point where it cannot be
>compiled.  There is no main, nor any function header for the function
>you're having problem with.  I don't see the line where you setup the
>precision, which is critical.

Maybe you didn't see my original post with the full code.  I
attach it again below...

I appreciate your help.

Julio.



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

int main(int argc, char **argv) {


	double Y[40], X[40], prod[40];
	double a=0.123456789012345;
	double b=0.98765432402340;
	double e=b;
	double c=0.1;
	double d=0.0;
	double f=0.93;
	int k;
	for(k=0;k<30;k++) {
	        d= a + pow(0.05 , k );
	        prod[k+1]=(a*e)-(d*b);
	        X[k+1]=(c*e-f*b)/prod[k+1];
	        Y[k+1]=(a*f-d*c)/prod[k+1];
	        printf("\n %lf \t %lf \t %lf", X[k+1],Y[k+1],prod[k+1]);
        }
}
else
{

  //Con precisión arbitraria
	mpf_t Y[40], X[40], prod[40];
	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;

 //inicialización
	 for(k=0;k< atoi(argv[2]);k++) {
	 mpf_init(Y[k]);
	 mpf_init(X[k]);
	 mpf_init(prod[k]);
	 }


	for(k=0;k< atoi(argv[2]);k++) {
	 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 );
	 mpf_set_d(a, 0.123456789012345);
	 mpf_set_d(b, 0.98765432402340);
	 mpf_set_d(e, 0.98765432402340);
	 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 );
        mpf_set_d(aux2, 0.05);
        mpf_pow_ui(aux1 , aux2, k );
        mpf_add(d, a, aux1) ;

       //prod[k+1]=(a*e)-(d*b);
        mpf_mul(aux2,a,e);
        mpf_mul(aux3,d,b);
        mpf_sub( prod[k+1], aux2,aux3);


        //X[k+1]=(c*e-f*b)/prod[k+1];
        mpf_mul(aux2,c,e);
        mpf_mul(aux3,f,b);
        mpf_sub(aux2,aux2,aux3);
        mpf_div( X[k+1], aux2, prod[k+1]);

        //Y[k+1]=(a*f-d*c)/prod[k+1];
        mpf_mul(aux2,a,f);
        mpf_mul(aux3,d,c);
        mpf_sub(aux2, aux2,aux3);
        mpf_div(Y[k+1],aux2,prod[k+1]);

        //printf("\n %lf \t %lf \t %lf", X[k+1],Y[k+1],prod[k+1]);
        printf("\niter %d ",k);
        gmp_printf(" %.Ff  \t %.Ff \t %.Ff",X[k+1],Y[k+1], prod[k+1] );

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