-O generates NaN

[Haruhiko Okumura <okumura at matsusaka-u dot ac dot jp>]

Hi,

Could someone tell why this code generates lots of NaN's when compiled
with any optimization flag (-O, -O2 etc) on?

$ g++ nan.cc         ... OK
$ g++ -O nan.cc      ... nan

$ g++ -v
Reading specs from /usr/local/lib/gcc-lib/i686-pc-linux-gnu/3.1/specs
Configured with: ../gcc/configure --enable-shared --enable-languages=c++,f77
gcc version 3.1 20010227 (experimental)

-- 
Haruhiko Okumura <okumura@matsusaka-u.ac.jp>
Matsusaka University, 1846 Kubo-cho, Matsusaka, 515-8511 Japan
Phone: +81-598-29-1122  Fax: +81-598-29-1014
http://www.matsusaka-u.ac.jp/~okumura/
// nan.cc
// g++ -O nan.cc generates NaN.
// g++ nan.cc is OK.
// Haruhiko Okumura <okumura@matsusaka-u.ac.jp>

#include <iostream>
#include <cstdlib>
#include <cmath>

using namespace std;

const double pi = 3.14159265358979323846;

// Specs of wire (strand) (in meters)

const double conductor_radius = 0.0001265;
const double wire_radius = 0.0001725;
const double twist_pitch1 = 0.028; // pitch of 1st-order wire
const double twist_pitch2 = 0.044; // pitch of 2nd-order wire (not used)

// Specs of Ultra-big Solenoid (lengths are in meters)

const double bobbin_radius = 0.183; // 366mm diameter
const double solenoid_length = 0.1; // 100mm length
const double turns = 74.5;      // inner 75, outer 74

// Parameters for computation

const int STEPRAD = 100;        // steps per radian

const double solenoid_pitch = solenoid_length / (turns + 1);
const double thetamax = 2 * turns * pi;
const int IMAX = int(STEPRAD * thetamax) + 1;

double xx[2][IMAX+1], yy[2][IMAX+1], zz[2][IMAX+1];
int maxiter = 20000;            // default maximum iterations in outerloop1()

double hypot(double x, double y, double z)
{
    return sqrt(x * x + y * y + z * z);
}

class Point {
public:
    double x, y, z;             // coordinates
    double dx, dy, dz;          // dx/dtheta, etc.
    int k;                      // 0 or 1
    void settheta(double theta) {
        int i = int(STEPRAD * theta);
        dx = STEPRAD * (xx[k][i+1] - xx[k][i]);
        dy = STEPRAD * (yy[k][i+1] - yy[k][i]);
        dz = STEPRAD * (zz[k][i+1] - zz[k][i]);
        double e = theta - double(i) / STEPRAD;
        x = xx[k][i] + dx * e;
        y = yy[k][i] + dy * e;
        z = zz[k][i] + dz * e;
    }
    void setk(int kk) {
        if (kk == 0 || kk == 1)
            k = kk;
        else
            cerr << "Illegal k" << endl;
    }
    Point() { k = 0; }
};

Point p1, p2;

double integrate(double left, double right, double integrand(double))
{
    double ans1 = 0;
    double ans2 = 0;
    double range = right - left;
    double sum = integrand(left) + integrand(right);
    double h = 0;
    for (int i = 2; i <= maxiter; i *= 3) {
        h = range / i;
        for (int j = 1; j < i; j++) {
            if (j % 3 == 0) {
                // do nothing
            } else if (j % 2 == 0) {
                sum += 2 * integrand(left + j * h);
            } else {
                sum += 4 * integrand(left + j * h);
            }
        }
        double ans = sum * h / 3;
        cout << " " << i << ": " << ans << endl;
        if (abs((ans - ans1) / ans) < 0.000001 &&
            abs((ans - ans2) / ans) < 0.000001) break;
        ans2 = ans1;
        ans1 = ans;
    }
    return sum * h / 3;
}

// total lengths of wire

double dlen(double t)
{
    p1.settheta(t);
    return hypot(p1.dx, p1.dy, p1.dz);
}

double totlen()
{
    return integrate(0, thetamax, dlen);
}

void wind0()
{
    double rmin, rmax, r, zmin, zmax, z;

    rmin = bobbin_radius + wire_radius;
    rmax = bobbin_radius + solenoid_pitch - wire_radius;
    zmin = wire_radius;
    zmax = solenoid_pitch - wire_radius;
    r = (rmin + rmax) / 2;
    z = (zmin + zmax) / 2;
    for (int i = 0; i <= IMAX; i++) {
        double theta = double(i) / STEPRAD;
        xx[0][i] = r * cos(theta);
        yy[0][i] = r * sin(theta);
        zz[0][i] = z + solenoid_pitch * theta / (2 * pi);
    }
    rmin = bobbin_radius + solenoid_pitch + wire_radius;
    rmax = bobbin_radius + 2 * solenoid_pitch - wire_radius;
    zmin = wire_radius;
    zmax = solenoid_pitch - wire_radius;
    r = (rmin + rmax) / 2;
    z = (zmin + zmax) / 2;
    for (int i = 0; i <= IMAX; i++) {
        double theta = double(i) / STEPRAD;
        xx[1][i] = r * cos(theta);
        yy[1][i] = r * sin(theta);
        zz[1][i] = solenoid_length - solenoid_pitch + z - solenoid_pitch * theta / (2 * pi);
    }
}

// #include <fpu_control.h>

int main()
{
    // fpu_control_t cw = 0x1372;  // interrupt on invalid/zero-divide/overflow
    // _FPU_SETCW(cw);

    wind0();
    p1.setk(0);
    double len1 = totlen();
    p1.setk(1);
    double len2 = totlen();
    cout << "len1 = " << len1 << ", len2 = " << len2 << endl;
}