I must point out that,
conceptually, using a loop instead of a bunch of calls to `pow()`, inline, is correct, but you are dealing with floating point here and floating point aren't exact representations of fractional values. This simple test can show this:
Code:
#include <stdio.h>
#include <math.h>
double calc_abv2(double abw)
{
const double a = -0.000039705486746795932;
const double b = 1.2709666849144778;
const double c = -0.40926819348115739;
const double d = 2.0463351302912738;
const double f = -7.8964816507513707;
const double g = 15.009692673927390;
const double h = -15.765836469736477;
const double i = 8.8142267038252680;
const double j = -2.0695760421183493;
return a +
b * abw +
c * pow(abw, 2) +
d * pow(abw, 3) +
f * pow(abw, 4) +
g * pow(abw, 5) +
h * pow(abw, 6) +
i * pow(abw, 7) +
j * pow(abw, 8)
;
}
double calc_abv3(double abw)
{
static const double coefficients[] = {
-0.000039705486746795932,
1.2709666849144778,
-0.40926819348115739,
2.0463351302912738,
-7.8964816507513707,
15.009692673927390,
-15.765836469736477,
8.8142267038252680,
-2.0695760421183493
};
double abv = coefficients[0];
double abw_mult = abw;
for (int i = 1; i < sizeof(coefficients) / sizeof(coefficients[0]); ++i)
{
abv = abv + coefficients[i] * abw_mult;
abw_mult *= abw;
}
return abv;
}
int main( void )
{
double r1, r2;
// 10.2 as a test! Both functions should give you the same value, isn't it?
r1 = calc_abv2( 10.2 );
r2 = calc_abv3( 10.2 );
printf( "r1 = %.30f\nr2 = %.30f\nr1 %s r2\n",
r1, r2, r1 == r2 ? "==" : "!=" );
}
If you compile and run you'll see:
Code:
$ cc -o test test.c -lm
$ ./test
r1 = -157417092.855879038572311401367187500000
r2 = -157417092.855879098176956176757812500000
r1 != r2
Another known fact is that
double precision floating point structure gives you a garantee only about 16
decimal algarisms of precision, everything else is rounded or truncated.
The correct value is -157417092.855879180651745689527008
I've got this by using a multiprecision calculator (bc) with this script:
Code:
/* calc.bc */
scale=50
a=-0.000039705486746795932
b=1.2709666849144778
c=-0.40926819348115739
d=2.0463351302912738
f=-7.8964816507513707
g=15.009692673927390
h=-15.765836469736477
i=8.8142267038252680
j=-2.0695760421183493
x=10.2
y=a+b*x+c*x^2+d*x^3+f*x^4+g*x^5+h*x^6+i*x^7+j*x^8
y
quit
r1 is off by 1.420794342881598205*10⁻⁷
r2 is off by 8.24747895127691955*10⁻⁸
Tiny errors, huh? And the loop function, in this case, is more "precise" because it isn't using "pow()". But, the point is: floating point is always an approximation and calculations must be done with care.
Take a look at those 9 constant above. Let's view their TRUE values, shall we:
Code:
#include <stdio.h>
void showDouble( double v )
{
struct dblstruc {
unsigned long f:52;
unsigned long e:11;
unsigned long s:1;
};
struct dblstruc *p = (struct dblstruc *)&v;
printf( "%70.65f -> S=%d E=%-4d F=%lu\n", v, (int)p->s, (int)p->e - 1023, (unsigned long)p->f );
}
int main( void )
{
int i;
static const double values[] = {
-0.000039705486746795932,
1.2709666849144778,
-0.40926819348115739,
2.0463351302912738,
-7.8964816507513707,
15.009692673927390,
-15.765836469736477,
8.8142267038252680,
-2.0695760421183493
};
for ( i = 0; i < sizeof values / sizeof values[0]; i++ )
showDouble( values[i] );
}
Compiling and running:
Code:
$ cc -o test test.c
$ ./test
-0.00003970548674679593185651849118755762901855632662773132324218750 -> S=1 E=-15 F=1355895991351192
1.27096668491447783999603871052386239171028137207031250000000000000 -> S=0 E=0 F=1220325461210661
-0.40926819348115739405358226576936431229114532470703125000000000000 -> S=1 E=-2 F=2869120707254850
2.04633513029127378501925704767927527427673339843750000000000000000 -> S=0 E=1 F=104337437756972
-7.89648165075137065116450685309246182441711425781250000000000000000 -> S=1 E=2 F=4387048327594962
15.00969267392738970556820277124643325805664062500000000000000000000 -> S=0 E=3 F=3946106164285136
-15.76583646973647745426205801777541637420654296875000000000000000000 -> S=1 E=3 F=4371777278915676
8.81422670382526796117872436298057436943054199218750000000000000000 -> S=0 E=3 F=458368884992823
-2.06957604211834933494174038060009479522705078125000000000000000000 -> S=1 E=1 F=156671318679056
S, F and E are integer values used to satisfy the equation for normalized floating point (IEEE-754).
[]s
Fred