Code crashes -- floating exception

[Rahul Siddharthan]

The following code crashes when it tries to write to a
file, and I can't figure out why. It says
Floating exception (core dumped)
at the fprintf statement, line 90. 
It works fine on other compilers. If I modify the
code so that it writes either suscfinal, or the other
expression, but not both, it works. If I write both with separate
fprintf expressions (even into separate files) it crashes.
If I comment out the fprintf expression it works.

Command line: cc -lm -O2 filename.c

If I don't optimize, it works.
I have been unable to simplify the code further and reproduce the
problem. For instance, if I just try assigning the calculated
values to the concerned variables and leave out all the
function calls which are supposed to calculate them, it works fine. 

VERSION
Red Hat Linux 5.0 on an alpha (Digital Personal Workstation 433a)
Linux kernel 2.0.30
gcc version egcs-2.91.57 19980901 (egcs-1.1 release)

I had reported a problem with slowness of some similar (but more
complicated) code some time back, and received a patch from Richard
Henderson which solved that problem. However, this seems to be
unrelated: it happens with the unpatched as well as patched compilers.

Rahul Siddharthan.



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


#define energy1(a,b) jsign*(1-(((a>>1) + (b>>3))&1)*2)  /* 3,1 */
#define energy2(a,b) jsign*(1-(((a>>1) + (b>>2))&1)*2)  /* 3,2 */
#define energy3(a,b) jsign*(1-(((a>>1) + b)&1)*2)       /* 3,4 */
#define energy4(a,b) jsign*(1-(((a>>2) + (b>>3))&1)*2)  /* 2,1 */
#define energy5(a,b) jsign*(1-((a + (b>>3))&1)*2)       /* 4,1 */
#define energy6(a,b) jsign*(1-((a + (b>>2))&1)*2)       /* 4,2 */


#define jsign 1
#define nn 6.45 
#define lmax 10
#define hoconst 6.022e23*1.38e-23/2048

#define randomint(a) floor(random()/(RAND_MAX+1.0)*a)

#define lx 8
#define ly 8
#define lz 8

main()
{
int a,b,c;
unsigned int arrindex;
unsigned long int ctr,n;
double E,S,T,Eavg,E2avg,Etot,E2tot,Savg,S2avg,Stot,S2tot,spfinal,suscfinal;
double energycell[16],spincell[16];
int lattice[lmax][lmax][lmax];

int randstate[256];
unsigned int randseed;
int randn;

double energylatticen();
double spinlattice();
double newstate();

FILE *ptvar;

randseed=134;
randn=256;
initstate(randseed,randstate,randn);


makespincell(spincell);
makeenergycell(energycell,spincell);


for (a=0; a<10; a++) {
  for (b=0; b<10; b++) {
   for (c=0; c<10; c++) {
     lattice[a][b][c]=9;
     if (((b+c)/2.0)!=floor((b+c)/2.0+0.1))
     lattice[a][b][c]=6;
   }
  }
}




E=energylatticen(lattice,energycell);

S=spinlattice(spincell,lattice);

T=300.0;
Etot=0; E2tot=0;
  for (n=1; n<=100; n++) {
  E=newstate(lattice,energycell,spincell,
                 E,&S,T,200,&Eavg,&E2avg,&Savg,&S2avg);
  if (n>40) { Stot=Stot+Savg; S2tot=S2tot+S2avg; 
               Etot=Etot+Eavg; E2tot=E2tot+E2avg;} 
  }
  Etot=Etot/60;
  E2tot=E2tot/60;
  Stot=Stot/60;
  S2tot=S2tot/60;
  suscfinal=(S2tot-Stot*Stot)/T/0.3/1.38e-16*6.022e23*2.25*(9.2741e-21)*(9.2741e-21)*36/2048; 
  spfinal=hoconst/T/T*(E2tot-Etot*Etot);
  ptvar=fopen("test","a");
  fprintf(ptvar,"%f %f %f \n",T*0.3,spfinal,suscfinal);
  fclose(ptvar);



}




makeenergycell(energycell,spincell)
double energycell[],spincell[];
{
int a;
for (a=0; a<=15; a++) {
if ((a==3)||(a==5)||(a==6)||(a==9)||(a==10)||(a==12))
 energycell[a]=-2*jsign;
 else if ((a==0) || (a==15))
 energycell[a]=6*jsign;
 else
 energycell[a]=0*jsign;
 energycell[a]=energycell[a];
}
}




makespincell(spincell)
double spincell[];
{
int a;
float m1[3],m2[3],m3[3],m4[3],cos1,cos2,cos3,cos4,hdir[3];

hdir[1]=1.0; 
hdir[2]=0.0;
hdir[3]=0.0;

m1[1]=sqrt(2.0/3.0);
m1[2]=sqrt(2.0/9.0);
m1[3]= 1.0/3.0;
m2[1]=-sqrt(2.0/3.0);
m2[2]=sqrt(2.0/9.0);
m2[3]=1.0/3.0;
m3[1]=0;
m3[2]=-sqrt(8.0/9.0);
m3[3]=1.0/3.0;
m4[1]=0.0;
m4[2]=0.0;
m4[3]=-1.0;
cos1=0;
cos2=0;
cos3=0;
cos4=0;
for (a=1; a<=3; a++){
 cos1=cos1+m1[a]*hdir[a];
 cos2=cos2+m2[a]*hdir[a];
 cos3=cos3+m3[a]*hdir[a];
 cos4=cos4+m4[a]*hdir[a];
}
for(a=0; a<=15;a++){
spincell[a]= (1-2*((a>>3)&1))*cos1 + (1-2*((a>>2)&1))*cos2
             +(1-2*((a>>1)&1))*cos3 + (1-2*(a&1))*cos4;
}
}







double energylatticen(lattice,energycell)
int lattice[][lmax][lmax];
double energycell[];
{
double elat;
int l1,l2,l3,l2m,l3m,a,a1,a2,a3,a4,a5,a6;
elat=0;
for (l1=0; l1<lx; l1++) {
 for (l2=0; l2<ly; l2++) {
  for (l3=0; l3<lz; l3++) {
   l2m=(l2-1+lx)&(lx-1);
   l3m=(l3-1+lx)&(lx-1);
   a=lattice[l1][l2][l3];
   a1=lattice[(l1+1)&(lx-1)][l2][l3];
   a2=lattice[l1][(l2+1)&(ly-1)][l3];
   a3=lattice[l1][l2][(l3+1)&(lz-1)];
   a4=lattice[(l1+1)&(lx-1)][l2m][l3];
   a5=lattice[(l1+1)&(lx-1)][l2][l3m];
   a6=lattice[l1][(l2+1)&(ly-1)][l3m];
   elat=elat+energycell[a]+energy1(a,a1)+energy2(a,a2)+energy3(a,a3)
      +energy4(a,a4)+energy5(a,a5)+energy6(a,a6);
  }
 }
}
return elat;
}  




double spinlattice(spincell,lattice)
int lattice[][lmax][lmax];
double spincell[];
{
 int l1,l2,l3;
 double stot;
 stot=0;
 for (l1=0; l1<lx; l1++){
  for (l2=0; l2<ly; l2++){
   for (l3=0; l3<lz; l3++){
   stot=stot+spincell[lattice[l1][l2][l3]];
   }
  }
 }
 return stot;
}





double newstate(lattice,energycell,spincell,E,S,T,Niters,Eavg,E2avg,
                Savg,S2avg)
double T,E,*S,*Eavg,*E2avg,*Savg,*S2avg;
double energycell[],spincell[];
int lattice[][lmax][lmax];
unsigned long int Niters;
/* Generates a new configuration by randomly changing the configuration
   of a randomly chosen tetrahedron, and does this Niter times */
{
  int anew,a,aold,atmp,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,l1,l2,l3,
             n,randint,an1;
  unsigned long int m;
  double dhn,dh,en1,en2,entot,en2tot,tw,twrand,spin,stot,s2tot;

 en1=E; entot=E; en2tot=E*E; spin=*S; stot=*S; s2tot=stot*stot;
 for (m=0; m<Niters; m++) {  

  l1=randomint(lx);
  l2=randomint(ly);
  l3=randomint(lz);  /* Coordinates of randomly chosen tetrahedron */

   
  a= lattice[l1][l2][l3];
  aold=a;
  randint=randomint(4);
  atmp=1 << randint;   /* New configuration of this tetrahedron */
  anew = (a&~atmp) | (~a&atmp); 

  
  if (atmp==8){

  a7=lattice[(l1-1+lx)&(lx-1)][l2][l3];
  a10=lattice[(l1-1+lx)&(lx-1)][(l2+1)&(ly-1)][l3];
  a11=lattice[(l1+lx-1)&(lx-1)][l2][(l3+1)&(lz-1)];
  
  dhn=energycell[a] +energy1(a7,a) +energy4(a10,a)+energy5(a11,a);

  a=anew;
  dhn=-dhn+energycell[a] +energy1(a7,a) +energy4(a10,a)+energy5(a11,a);

  }
  else
  if (atmp==4){
  
  a4=lattice[(l1+1)&(lx-1)][(l2-1+ly)&(ly-1)][l3];
  a8=lattice[l1][(l2-1+ly)&(ly-1)][l3];
  a12=lattice[l1][(l2-1+ly)&(ly-1)][(l3+1)&(lz-1)];
  
  dhn=energycell[a] +energy4(a,a4) +energy2(a8,a) +energy6(a12,a);
  a=anew;
  dhn=-dhn+energycell[a] +energy4(a,a4) +energy2(a8,a) +energy6(a12,a);

  }
  else
  if (atmp==2){
  
  a1=lattice[(l1+1)&(lx-1)][l2][l3];
  a2=lattice[l1][(l2+1)&(ly-1)][l3];
  a3=lattice[l1][l2][(l3+1)&(lz-1)];
  
  dhn=energycell[a]+energy1(a,a1)+energy2(a,a2)+energy3(a,a3);

  a=anew;
  dhn=-dhn+ energycell[a]+energy1(a,a1)+energy2(a,a2)+energy3(a,a3);

  }
  else
  {
  a5=lattice[(l1+1)&(lx-1)][l2][(l3-1+lz)&(lz-1)];
  a6=lattice[l1][(l2+1)&(ly-1)][(l3-1+lz)&(lz-1)];
  a9=lattice[l1][l2][(l3-1+lz)&(lz-1)];
  
  dhn=energycell[a]+ energy5(a,a5)+energy6(a,a6)+energy3(a9,a);

  a=anew;
  dhn= -dhn+energycell[a]+ energy5(a,a5)+energy6(a,a6)+energy3(a9,a);

  }

  if ((dhn)>0) {
    tw=exp(-(dhn)/T);
    twrand=random()/(RAND_MAX+1.0);
    if (twrand<=tw) {
      en1=en1+dhn; 
      lattice[l1][l2][l3]=anew; 
      spin=spin-spincell[aold]+spincell[anew];
      }
    }
  else {
    en1=en1+dhn;
    lattice[l1][l2][l3]=anew;
    spin=spin-spincell[aold]+spincell[anew];
  }
  entot=entot+en1;
  en2tot=en2tot+en1*en1;
  stot=stot+spin;
  s2tot=s2tot+spin*spin;
  }
  

  *Eavg=entot/(Niters+1.0);
  *Savg=stot/(Niters+1.0);
  *E2avg=en2tot/(Niters+1.0);
  *S2avg=s2tot/(Niters+1.0);
  *S=spin; 
  return en1;
}