An interesting object lesson (Ada vs C)

[dewar@merv.cs.nyu.edu (Robert Dewar)]

Recently during one of my classes (on microprocessor architecture, where
I asked them to do tests on some particular aspect of MP performance), I
received the message below from a student who really knows what he is
doing, and got completely stumped, to the point, where, as you will see
from the message, he was sure he had found a gcc or Pentium bug. These
kind of dangers that come from using an untyped language waste a huge
amount of time. Needless to say if the student had used Ada, he would
have saved himself a lot of time.

I will leave it as an excercise for the reader to see the error in
the student's C code:

Classmates & Professor Dewar,

I am comparing the performance of floating point for single, double and
extended double precision numbers.  I am using the code below, with only the
types changed for the three different versions.  This version, the extended
double (long double) version, is returning strange results.  In fact, I'm
convinced I've found a bug in either gcc or the Pentium.

Here is the code:

#include <stdio.h>
#include <time.h>

long double mysqrt( long double number );

long double a, b, y, x = 0.5;
int c;
clock_t ticks1, ticks2;

int main() {
    ticks1 = clock();
    for ( c = 0; c < 100000; c++ ) {
	a = 1.0 / mysqrt(1.0 + x*x);
	b = 1.0;
	while ( (a - b < 0.0 ? b - a : a - b) > 10e-9 ) {
	    a = 0.5 * (a + b);
	    b = mysqrt(a * b);
	}
	y = x / ( mysqrt(1.0 + x*x) * a );
    }
    ticks2 = clock();
    printf( "arctan %g = %g\n", x, y );
    printf( "Clock ticks:  %d\n",
	ticks2 - ticks1 );

    return 0;
}

long double mysqrt( long double number ) {
    long double r = number / 2.0;
    long double old_r = 0.0;
    while ( ( r > old_r ? r - old_r : old_r - r ) / r > 10e-9 ) {
	old_r = r;
	r = 0.5 * ( r + number / r );
    }

    return r;             
}

Notice that x is only assigned to once, when it is initialized.

Here are the results:

arctan 0 = 3.25435e-229
Clock ticks:  2160

How did x get set to 0!?!?!?!?