The following code demonstrates a problem I'm having where I occasionally get
NaN results when I compile without optimization, but which works correctly when
I compile with optimization.  Compiling with "gcc -o trial trial.c -lm" I get:

angle: 180
angle: nan
angle: nan
angle: nan

Compiling with "gcc -O -o trial trial.c -lm" I get:

angle: 180
angle: 180
angle: 180
angle: 180


This code behaves the same on a RedHat-9 box on an AthlonXP and on a RedHat 7.3
box on a Centaur VIA Samuel 2.  On the RedHat-9 box "gcc -v" gives:

Reading specs from /usr/lib/gcc-lib/i386-redhat-linux/3.2.2/specs
Configured with: ../configure --prefix=/usr --mandir=/usr/share/man
--infodir=/usr/share/info --enable-shared --enable-threads=posix
--disable-checking --with-system-zlib --enable-__cxa_atexit --host=i386-redhat-linux
Thread model: posix
gcc version 3.2.2 20030222 (Red Hat Linux 3.2.2-5)

On the Redhat-7.3 box "gcc -v" gives:

Reading specs from /usr/lib/gcc-lib/i386-redhat-linux/2.96/specs
gcc version 2.96 20000731 (Red Hat Linux 7.3 2.96-113)


Here's the code:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
 
 
double vector_angle1(const double *a, const double *b);
double vector_angle2(const double *a, const double *b);
 
 
int main(int argc, char *argv[])
{
        double r1[3], r2[3];
 
        r1[0] = -3928.0498049187122;
        r1[1] = -2499.319958930655;
        r1[2] = 0;
        r2[0] = 0.84369535291866071;
        r2[1] = 0.53682227176548536;
        r2[2] = 0;
 
        printf("angle: %g\n", 180.0*vector_angle1(r1, r2)/M_PI);
        printf("angle: %g\n", 180.0*vector_angle2(r1, r2)/M_PI);
 
        r1[0] = -3697.7209740416451;
        r1[1] = -2352.7674780004572;
        r1[2] = 0;
        r2[0] = 0.84369535291866071;
        r2[1] = 0.53682227176548536;
        r2[2] = 0;
 
        printf("angle: %g\n", 180.0*vector_angle1(r1, r2)/M_PI);
        printf("angle: %g\n", 180.0*vector_angle2(r1, r2)/M_PI);
 
        exit(EXIT_SUCCESS);
}
 
 
double vector_angle1(const double *a, const double *b)
/*
 * Calculates the angle between a & b and returns it (in rad).
 * Assumes |a| & |b| are not zero.
 */
{
        double l;
        l = sqrt((a[0] * a[0] + a[1] * a[1] + a[2] * a[2]) *
                 (b[0] * b[0] + b[1] * b[1] + b[2] * b[2]));
        return acos((a[0] * b[0] + a[1] * b[1] + a[2] * b[2]) / l);
}
 
 
double vector_angle2(const double *a, const double *b)
/*
 * Calculates the angle between a & b and returns it (in rad).
 * Assumes |a| & |b| are not zero.
 */
{
        return acos((a[0] * b[0] + a[1] * b[1] + a[2] * b[2]) /
                    (sqrt((a[0] * a[0] + a[1] * a[1] + a[2] * a[2]) *
                          (b[0] * b[0] + b[1] * b[1] + b[2] * b[2]))));
}