sqrt Behaving Unexpectedly --- Short Program Example

[Keith]

The bits get twiddled every now and then.
Using gcc 3.2 and 3.4.

To build: cc -o test test.c -lm

---------------- 8< test.c ------------------------
#include <stdio.h>
#include <math.h>

double func(double);

int main() {

int ii ;
double x, r1, r2 ;

x = 2510771730459.374023 ;

for ( ii = 0 ; ii < 100 ; ii++ ) {

r1 = sqrt( 3*x*x ) ;
r2 = sqrt(func(3*x)*x);

if ( r1 != r2 ) {
fprintf(stderr, "\nIteration=%d\n", ii);
fprintf(stderr,"val1=%#llx\n",
*(unsigned long long*)(&r1)) ;
fprintf(stderr,"val2=%#llx\n",
*(unsigned long long*)(&r2)) ;
}
x = r1 ;
}

return 0 ;
}

double func( double temp ) {
return(temp);
}

 

[Keith]

Simpler program that shows the same thing:

-----------------8<-------------------------
#include <stdio.h>
#include <math.h>

int main() {

int ii ;
double x, y, r1, r2 ;

x = 2510771730459.374023 ;

for ( ii = 0 ; ii < 100 ; ii++ ) {

r1 = sqrt( 3*x*x ) ;
y = 3*x ;
r2 = sqrt(x*y);

if ( r1 != r2 ) {
fprintf(stderr, "\nIteration=%d\n", ii);
fprintf(stderr,"val1=%#llx\n",
*(unsigned long long*)(&r1)) ;
fprintf(stderr,"val2=%#llx\n",
*(unsigned long long*)(&r2)) ;
fprintf(stderr,"r=%lf\n", 3*r2);
}
x = r1 ;
}

return 0 ;
}
----------------------------------------------------

 

[Eric Sosman]

The bits get twiddled every now and then.
Using gcc 3.2 and 3.4.
[code snipped; see up-thread]

Part of what you're seeing is covered by Questions 14.1
and 14.4 of the comp.lang.c Frequently Asked Questions (FAQ)
list

http://www.eskimo.com/~scs/C-faq/top.html

In floating-point arithmetic, expressions that are
"algebraically identical" can sometimes yield computationally
different results. The whys and wherefores stray into the
particular characteristics of different C implementations,
which aren't really topical here. Suffice it to say that
on some implementations, in

double x, y, z;
x = ...;
y = 3 * x * x;
z = 3 * x;
z = z * x;

.... `y' and `z' will not have exactly the same value. On
other machines, they will: I ran your code on two different
machines and got equality on one, "twiddled bits" on the
other. Be comforted (and read Question 14.4 again): the
difference between `y' and `z', while not necessarily zero,
will be very small[*].

[*] "Very small" for these expressions, but not for all.
For example, try

double x = 1.0, y = 1.0e30;
printf ("%g\n%g\n", (x + y) - y, x + (y - y));

If you want to learn more about such effects, why they arise,
and how to avoid being bitten by them, read the first few
chapters of any textbook with the phrase "Numerical Analysis"
in the title.

 

[Keith]

Thanks Eric. I wrote this to the person that asked me the question
initially about a change in simulation response. I do not know what
group to pose gcc differences to.

I think that the key is where Eric said,
In floating-point arithmetic, expressions that are
"algebraically identical" can sometimes yield computationally
different results. The whys and wherefores stray into the
particular characteristics of different C implementations...


I suppose that is to say that gcc 3.2 and 3.4 (two different C
implementations) did something different which changed something
considered to be fundamental.

The change is "small[*]" as he put it, but neverthless makes
significant differences. It is especially significant when comparing
against a validated simulation response! Especially when a pointy head
is staring at you with an eyebrow raised.

I don't really know how to accept the changes that are imposed on us by
things we don't control. It is unsettling that a compiler upgrade
changes the angle of an end effector by a 10th of a degree or whatever.
What are you going to do?

Hmmm. I suppose you have to ask what is truth. Seems like things aren't
so black and white. How was the first response "vailidated"? How did it
become "accepted"?

Seems like the weather to me. Butterfly wings causing weather changes.

If I was the ruler of the world I might suggest that when considering
simulation responses and issues of bit twiddling changing responses,
that modelers should keep the response within an epsilon. After all, it
is simulation.

If the validated response IS reality, and NO deviation can occur then
we are all in a heap of trouble. RK2 or RK4 or the best stinking
look-ahead-multi-parallel-heuristic-stochastic integrator in the
universe will fail. We'll have to scramble for f(t) which describes the
station and shuttle in time. And we'll have no need for Trick. We just
need f(t). And I won't have a job.

It is chaos. The nature of things.
Yep. Too much green tea this morning.

Keith

 

[Eric Sosman]

Thanks Eric. I wrote this to the person that asked me the question
initially about a change in simulation response. I do not know what
group to pose gcc differences to.

I think that the key is where Eric said,
In floating-point arithmetic, expressions that are
"algebraically identical" can sometimes yield computationally
different results. The whys and wherefores stray into the
particular characteristics of different C implementations...

I suppose that is to say that gcc 3.2 and 3.4 (two different C
implementations) did something different which changed something
considered to be fundamental.

What's *probably* going on is that the two compilers
generate different code sequences for the floating-point
calculations. Some machines use more precision for F-P
values in their internal registers than when those values
are stored back into memory; x86 designs, for example, use
80-bit F-P while calculating but round the result to 64
bits when storing it back to a `double' variable. Hence
the compiler's decisions about which intermediate results
to hold in registers and which to round and store back to
memory can affect the outcome of the computation. Change
the compiler version, change the optimization level, write
expressions in different although "algebraically equivalent"
ways, and you'll get different decisions from the compiler.

The change is "small[*]" as he put it, but neverthless makes
significant differences. It is especially significant when comparing
against a validated simulation response! Especially when a pointy head
is staring at you with an eyebrow raised.

The "small" effect I observed (on an x86 machine, as it
happens) amounted to a change in the least significant bit
of the result. If a relative error of fifty-five billionths
of a billionth (55.5e-18) raises eyebrows, the heads those
eyebrows decorate must come to very sharp points indeed.

Hmmm. I suppose you have to ask what is truth. Seems like things aren't
so black and white. How was the first response "vailidated"? How did it
become "accepted"?

"That is the question." -- Hamlet, Prince of Denmark

 

[Lawrence Kirby]

On Wed, 02 Feb 2005 08:05:30 -0800, Keith wrote:

....

I suppose that is to say that gcc 3.2 and 3.4 (two different C
implementations) did something different which changed something
considered to be fundamental.

The change is "small[*]" as he put it, but neverthless makes
significant differences.

If the differences are significant enough to be a problem then the code is
probably faulty to start with. Perhaps data is not being held to a high
enough precision for the purpose, or algorithms and formulae are being
used that have bad behaviour in this respect.

It is especially significant when comparing
against a validated simulation response! Especially when a pointy head
is staring at you with an eyebrow raised.

Where floating point is concerned a validated result should be specified
within error bounds. If your new result is within those bounds there is no
problem, if it isn't then the program is probably buggy.

I don't really know how to accept the changes that are imposed on us by
things we don't control. It is unsettling that a compiler upgrade
changes the angle of an end effector by a 10th of a degree or whatever.
What are you going to do?

Hmmm. I suppose you have to ask what is truth. Seems like things aren't
so black and white. How was the first response "vailidated"? How did it
become "accepted"?

Right. Is there any reason to believe that the original result is "better"
than the new one?

Seems like the weather to me. Butterfly wings causing weather changes.

Don't believe this myself, a butterfly wingbeat will cause the position
of individual molecules in the air to be increasingly different as time
goes on, an effect which will propagate around the world, but the
difference in macroscopic effects such as those responsible for weather
are likely to decrease. OTOH errors almost invariably increase as you
progress through a calculation.

If I was the ruler of the world I might suggest that when considering
simulation responses and issues of bit twiddling changing responses,
that modelers should keep the response within an epsilon. After all, it
is simulation.

Error analysis need to be performed, yes. Without that there's no
guarantee that the results have any meaning at all.

....

Lawrence

 

[Keith]

If a relative error of fifty-five billionths
of a billionth (55.5e-18) raises eyebrows, the heads those
eyebrows decorate must come to very sharp points indeed.

To defend the pointy heads, which I am not likely to do very often, the
billionths wouldn't raise the eyebrow of the most pointiest of pointy
heads. The problem is that over time the algorithms slowly begin to
change the outcome to a semi-significant level. At this
semi-significant level, the pointy head's eyebrow rises. Defending the
pointy head again, the response is important and may deserve attention.
The answer to the pointy heads is more philosophic when it comes to
what you consider to be a "truth model" when using imprecise algorithms
and imprecise machinery. I have already posted my thoughts over the
nature of what is "right" so no need to say that the discrete nature of
how computers represent real numbers make programs imprecise.

The beginning issue was whether the precision should change from C
compiler to C compiler with such a seemingly innocuous issue as the
3*x*x deal. You've answered that question when you said that execution
is compiler implementation dependent. So that helps a lot. The
details you gave about what is *probably* going on between the
compilers help too. Thanks.

 

[Keith]

Error analysis need to be performed, yes. Without that there's no
guarantee that the results have any meaning at all.

Yep. That sums it up.

changes.

Don't believe this myself, a butterfly wingbeat
will cause the position
of individual molecules in the air to be increasingly
different as time goes on, an effect which will propagate
around the world, but thedifference in macroscopic effects
such as those responsible for weather are likely to decrease.
OTOH errors almost invariably increase as you
progress through a calculation.

That chaos example of a butterfly seems far fetched, huh?
I have seen the error increase with "higher fidelity" algorithms just
because of calculation induced error.

The fact that the x86 uses 80 bit registers *could* explain why there
are differences between a function call and a macro. The function call
would lose precision since the double passed would be reduced to fit on
the function stack.

 

[Gregory Toomey]

The bits get twiddled every now and then.
Using gcc 3.2 and 3.4.

To build: cc -o test test.c -lm

---------------- 8< test.c ------------------------
#include <stdio.h>
#include <math.h>

double func(double);

int main() {

int ii ;
double x, r1, r2 ;

x = 2510771730459.374023 ;

for ( ii = 0 ; ii < 100 ; ii++ ) {

r1 = sqrt( 3*x*x ) ;
r2 = sqrt(func(3*x)*x);

if ( r1 != r2 ) {
fprintf(stderr, "\nIteration=%d\n", ii);
fprintf(stderr,"val1=%#llx\n",
*(unsigned long long*)(&r1)) ;
fprintf(stderr,"val2=%#llx\n",
*(unsigned long long*)(&r2)) ;
}
x = r1 ;
}

return 0 ;
}

double func( double temp ) {
return(temp);
}

Double has approx 15 significant digits in the IEEE floating point standard,
used by all modern computers. You will get problems like this as the
calculations differ in one or two bits.

Instead compare the absolute value of the difference.

gtoomey

 

[Lawrence Kirby]

Yep. That sums it up.


That chaos example of a butterfly seems far fetched, huh?

Just a misapplication of the principles of Chaos Theory. A little
knowledge is often far more dangerous than no knowledge.

I have seen the error increase with "higher fidelity" algorithms just
because of calculation induced error.

That's why error analysis is so important. A result is meaningless unless
you know what the error bounds are.

The fact that the x86 uses 80 bit registers *could* explain why there
are differences between a function call and a macro. The function call
would lose precision since the double passed would be reduced to fit on
the function stack.

In the error analysis you would assume the lowest possible precision at
all steps. If the precision happens to be higher it's not a problem,
except in the case of some specific algorithms that require a consistent
precision to work.

Lawrence