[gcc-help] variables set to NaN. Can't understand why.

[gcc-help] variables set to NaN. Can't understand why.

[Bernard Godard]

Hello,

I'm having problem porting a C program from SPARC/Solaris to x86/Linux.

The program compiles but sometimes (for some values of the inputs) crashes due 
to some variables being set to NaN (Not a Number). This problem occured on 
Slackware9.1/gcc3.2.3/Pentium4 and Mandrake9.2/gcc3.3.1/Pentium4(the same 
machine) but does not seem to happen with Redhat7.3/gcc2.96/Pentium3.

I tried to debug it with gdb, with no optimizations.

I couldn't really isolate the error (See the end of this message).

The description of what went wrong is at the end of this message.

Here is the incriminated piece of code (sorry for the length):

***************************************************
/*
 * Name:     SGP4_init
 * Purpose:  initializes time-invariant SGP4 variables
 * Input:    SGP4 sgp4
 * Output:   none
 */
void SGP4_init(SGP4 sgp4) {
  double a1, del1, a0, del0, perigee, sStar, eta_2, eta_3, eta_4, pinvsq, 
    temp1, temp2, temp3, x3thetam1, x1m5th, xhdot1, theta_4, zeta, zeta_2, 
    zeta_3, zeta_5, beta0, beta0_3, beta0_4, C2;

  /* recover original mean motion and semimajor axis from the elements */
  a1=pow(KE/sgp4->n0, TWOBYTHREE);
  del1=1.5 * (CK2/(a1*a1))*((3*cos(sgp4->i0)*cos(sgp4->i0)-1.0)/
    pow(1-sgp4->e0*sgp4->e0, 1.5));
  a0=a1*(1-del1/3.0-del1*del1-(134.0/81.0)*del1*del1*del1);
  del0=(del1*a1*a1)/(a0*a0);

  /* find original mean motion, n0dp (n0'') */
  sgp4->n0dp=sgp4->n0/(1+del0); 

  /* find semimajor axis, a0dp (a0'') */
  sgp4->a0dp=a0/(1-del0);

  /* find the perigee (in km) */
  perigee=(sgp4->a0dp*(1.0-sgp4->e0)-AE)*KM_PER_EARTH;
  
  /* make decisions on what value of s to use based on perigee */
  if (perigee<=156.0 && perigee>=98.0) {
    sStar=sgp4->a0dp*(1.0-sgp4->e0)-S-AE;
    sgp4->q0ms4=pow(pow(Q0MS4, 0.25) + S - sStar, 4);
    }
  else if (perigee<98.0) {
    sStar=20/KM_PER_EARTH + AE;
    sgp4->q0ms4=pow(pow(Q0MS4, 0.25) + S - sStar, 4);
    }
  else {
    sStar=S;
    sgp4->q0ms4=Q0MS4;
    }

  sgp4->theta=cos(sgp4->i0); sgp4->theta_2=sgp4->theta*sgp4->theta; 
  theta_4=sgp4->theta_2*sgp4->theta_2;
  zeta=1/(sgp4->a0dp-sStar);

  

  zeta_2=zeta*zeta; zeta_3=zeta_2*zeta; 

 

  sgp4->zeta_4=zeta_3*zeta;
  zeta_5=sgp4->zeta_4*zeta;

  sgp4->beta0_2=1-sgp4->e0*sgp4->e0;
  beta0=sqrt(sgp4->beta0_2); beta0_3=sgp4->beta0_2*beta0; 
  beta0_4=beta0_3*beta0;


  sgp4->eta=sgp4->a0dp*sgp4->e0*zeta;
  eta_2=sgp4->eta*sgp4->eta; eta_3=eta_2*sgp4->eta; eta_4=eta_3*sgp4->eta;
  
  C2=sgp4->q0ms4*(sgp4->zeta_4)*sgp4->n0dp*pow(1-eta_2, -7.0/2.0)*
    (sgp4->a0dp*(1+ ((3.0/2.0)*eta_2) + 4*sgp4->e0*sgp4->eta + 
    sgp4->e0*eta_3) + (3.0/2.0)*(CK2*zeta/(1-eta_2))*(-0.5 + 
    1.5*sgp4->theta_2)*(8.0+24*eta_2 + 3*eta_4));
  sgp4->C1=sgp4->bstar*C2; sgp4->C1_2=sgp4->C1*sgp4->C1; 
  sgp4->C1_3=sgp4->C1_2*sgp4->C1; sgp4->C1_4=sgp4->C1_3*sgp4->C1;
  sgp4->C3=(sgp4->q0ms4*(zeta_5)*A30*sgp4->n0dp*AE*sin(sgp4->i0))/
    (CK2*sgp4->e0);
  sgp4->C4=2*sgp4->n0dp*sgp4->q0ms4*(sgp4->zeta_4)*sgp4->a0dp*(sgp4->beta0_2)*
    pow(1-eta_2, -7.0/2.0)*( (2*sgp4->eta*(1+sgp4->e0*sgp4->eta)+
    0.5*sgp4->e0+0.5*eta_3)-2*CK2*zeta/(sgp4->a0dp*(1-eta_2))*
   (3*(1-3*sgp4->theta_2)*(1 + 1.5*eta_2 -2*sgp4->e0*sgp4->eta - 
   0.5*sgp4->e0*eta_3) + 0.75*(1-sgp4->theta_2)*(2*eta_2 - 
   sgp4->e0*sgp4->eta - sgp4->e0*eta_3)*cos(2*sgp4->w0)));
  sgp4->C5=2*sgp4->q0ms4*(sgp4->zeta_4)*sgp4->a0dp*(sgp4->beta0_2)*
    pow(1-eta_2, -7.0/2.0)*(1 + (11.0/4.0)*sgp4->eta*(sgp4->eta+sgp4->e0)+
    sgp4->e0*eta_3);

  sgp4->D2=4*sgp4->a0dp*zeta*sgp4->C1_2;
  sgp4->D3=(4.0/3.0)*sgp4->a0dp*zeta_2*(17*sgp4->a0dp + 
    sStar)*sgp4->C1_3;
  sgp4->D4=TWOBYTHREE*sgp4->a0dp*zeta_3*(221*sgp4->a0dp+
    31*sStar)*sgp4->C1_4;

  /* constants for secular effects calculation */
  pinvsq=1.0/(sgp4->a0dp*sgp4->a0dp*beta0_4);
  temp1=3.0*CK2*pinvsq*sgp4->n0dp;
  temp2=temp1*CK2*pinvsq;
  temp3=1.25*CK4*pinvsq*pinvsq*sgp4->n0dp;
  x3thetam1=3.0*sgp4->theta_2 - 1.0;
  x1m5th=1.0 - 5.0*sgp4->theta_2;
  xhdot1=-temp1*sgp4->theta;

  
  sgp4->mdot=sgp4->n0dp+0.5*temp1*beta0*x3thetam1+
    0.0625*temp2*beta0*(13.0 - 78.0*sgp4->theta_2+137.0*theta_4);

  sgp4->wdot=-0.5*temp1*x1m5th+0.0625*temp2*(7.0-114.0*sgp4->theta_2+
    395.0*theta_4)+temp3*(3.0-36.0*sgp4->theta_2+49.0*theta_4);
  
  sgp4->raandot=xhdot1+(0.5*temp2*(4.0-19.0*sgp4->theta_2)+2.0*temp3*
    (3.0-7.0*sgp4->theta_2))*sgp4->theta;
  }

***************************************************

And here is the sgp4 definition

***************************************************
/* the data structure itself */
struct sgp4_st {
  /* NORAD tle components */
  double epochTime, n0, n0dt, n0dt2, bstar, i0, raan, e0, w0, M0;

  /* derived variables */
  double n0dp, a0dp, q0ms4, theta, theta_2, zeta_4, beta0_2, eta, C1, 
    C1_2, C1_3, C1_4, C3, C4, C5, D2, D3, D4, mdot, wdot, raandot; 
  };

typedef struct sgp4_st * SGP4;
***************************************************

and memory allocation:

***************************************************
/*
 * Name:     SGP4_create
 * Purpose:  creates an instance of SGP4
 * Input:    none
 * Output:   SGP4 
 */
SGP4 SGP4_create(void) {
  return((SGP4)calloc(1, sizeof(struct sgp4_st)));
  }

***************************************************

The sgp4 variable has its memory allocated by SGP4_create. Then some of its 
members are assigned values by SGP4_set (Not shown in this message). SGP4_init 
(the piece of code that causes troubles) computes the other members of the sgp4 
structure.

Now for the funny part:

With the debugger, I saw that the values of the variables were correctly set 
except for C2 which was set to NaN and all the variables depending on C2.

I cut the calculation of C2 into several steps. Then C2 was correctly set and so 
were the variables depending on C2. But then sgp4->mdot was set to NaN (it was 
previously correct).

I cut the calculation of sgp4->mdot into two steps. Then sgp4->mdot was set 
correctly but sgp4->wdot was set to NaN (It was correct before).

Anyone has an idea of what is wrong in the code? Could this be a compiler bug?

Thank you.

     Bernard GODARD