Hidden error in c program ( numerical computation )

Question

Here are 2 functions. I see no important difference between them, but the results are different. Why ?

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

long double GiveCx(long double x)
{  
double a = 5.3226927610784935E-01;
double b = 6.5410208763684241E-01;
double c = -1.4312869957125389E+00;
double d = 8.4710834303177074E-01;
return (c*atanl(expl(x-a)/b) + d);
}

double BurkardtCollectionBased_sech_cdf_Offset_model(double x_in)
{
 double temp;
 temp = 0.0;
 // coefficients
 double a = 5.3226927610784935E-01;
 double b = 6.5410208763684241E-01;
 double c = -1.4312869957125389E+00;
 double Offset = 8.4710834303177074E-01;
 temp = c * atan(exp((x_in-a)/b));
 temp += Offset;
 return temp;
}

int main()
{
 int ix;
 for (ix=0; ix<5; ix++)
  printf(" ix = %d ; c = %.20f  ;  %.20Lf \n", ix,      BurkardtCollectionBased_sech_cdf_Offset_model( (long double)ix), GiveCx((long double )ix));
return 0;
}

Results are :

ix = 0 ; c = 0.25000064588764425721 ; -0.20004050665796930359

ix = 1 ; c = -0.75000921765452766010 ; -0.84455584419096496618

ix = 2 ; c = -1.24993551212417064455 ; -1.18701705113792041978

ix = 3 ; c = -1.36825264215735509232 ; -1.32186430467910977205

ix = 4 ; c = -1.39401846938445195256 ; -1.37195787094497580628

Answer 1

You have expl(x-a)/b in one and exp((x_in-a)/b) in the other. Check if b needs to divide the exponent like in the latter case, or the result like in the former case.

You probably need expl((x-a)/b) in GiveCx().

Answer 2

The primary issue is as ja72 states, but note that, depending on your compiler, you may still end up with slightly differing results even if the functions are uniformalized.

http://ideone.com/f29ugQ

ix = 0 ; c = -0.20004050665796929698 ; -0.20004050665796930359
ix = 1 ; c = -0.84455584419096485504 ; -0.84455584419096496618
ix = 2 ; c = -1.18701705113792033153 ; -1.18701705113792041978
ix = 3 ; c = -1.32186430467910986941 ; -1.32186430467910977205
ix = 4 ; c = -1.37195787094497578806 ; -1.37195787094497580628

The cause is the differing precision of double and long double in compilers that don't alias long double to double and the inherent imprecision of floating-point computations. Even though the supplied inputs (aside from x versus x_in) have the same precision, you are using functions that perform calculations with differing precision (atanl and expl versus atan and exp). atanl and expl will both widen double inputs to long double and may also perform their calculations using that level of precision, depending on the compiler, while atan and exp should only perform calculations with double precision (though some compilers may have the option to make the computations more precise than the types of the inputs, as that could still result in more accurate results depending on the calculations performed).

This is why, if you need perfect decimal precision, you have to use fixed-point libraries that support any number of decimal places (though irrational values are still beyond such capabilities).