🚀 go-pugleaf

in the same vein as my earlier bug report, i ran into the following problem.
when line 222 in the following code is uncommented i get the following
compiler crash

================================================================================
fft.c: In function `main':
fft.c:226: internal error--insn does not satisfy its constraints:
(insn 4175 741 716 (set (reg:SI 68)
       (mem:SI (plus:SI (reg:SI 30)
               (const_int -380)))) 34 {movsi+3} (nil)
   (nil))
gcc: Internal compiler error: program cc1 got fatal signal 6
================================================================================

the offending (apparently) line is the call to exit(0);

================================================================================
#include <stdio.h>
#include <math.h>

#define L	32
#define LM	(L-1)
#define L3	(L*L*L)
#define FTOL	2.0e-5

main()
{
  extern void fourn();

  int a,b,c,d,e,f,r12;
  int nxm,nxp,nym,nyp,nzm,nzp,nxyz;
  int np, npl, nmin, it, itave;		/* number of poly balls */
  int density;				/* counter-ion density */
  int nimp;				/* number of counter ions */
  int locpol[3][20]; 			/* location of poly balls */
  int pbp[L],pbm[L];			/* indices for PB cond's */
  int nn[3], dim, sign;
  double amf,action,feqn,nit;
  double sumexp, sumexm, gammap, gammam;
  double cp, cm, ep, ef, phiold;
  double t1, t2;
  double phi[L3],vexc[L3];	 	/* field, poly-repulsions */
  double beta, z;
  double pil,il3;
  float dphi[L3][2],psq[L3];
  double eps,psqcut;

  pil = M_PI/L;
  il3 = 1.0/L3;
  fscanf(stdin,"%d\n",&np);
  for (a=0;a<np;a++)
    fscanf(stdin,"%d %d %d\n",&locpol[0][a],&locpol[1][a],&locpol[2][a]);
  np = 3;

  for (a=0;a<L;a++) {
    pbp[a] = (a+1)%L;
    pbm[a] = (a+LM)%L;
  }
  fscanf(stdin,"%le %le\n",&beta,&z);
  fscanf(stdin,"%d\n",&density);
  nimp=density*L3;
  /* nimp = 4000; */

  fscanf(stdin,"%le %le\n",&eps,&psqcut);
/* initialize needed stuff for fft acceleration */
  nn[0] = L;
  nn[1] = L;
  nn[2] = L;
  dim = 3;
  for (a=0;a<L;a++)
    for (b=0;b<L;b++)
      for (c=0;c<L;c++)
	if ( a+b+c == 0 ) psq[0] = 1.0;
	else psq[a+b*L+c*L*L]= 3.0/
	  (pow(sin(a*pil),2.0)+pow(sin(b*pil),2.0)+pow(sin(c*pil),2.0)+psqcut);

  nit = 0.0; itave = 0;
  for (a=0;a<L3;a++) vexc[a] = 0.0;
  locpol[0][0] = locpol[1][0] = locpol[2][0] = 0;
  locpol[0][1] = locpol[1][1] = locpol[2][1] = 0;
  locpol[0][2] = locpol[1][2] = locpol[2][2] = 0;

/* the first two colloids are position along the x axis, and the third is
   placed on all the other lattice points in the 1/4th plane nearest the
   orign.  the configurations with both "extra" colloids along the x axis
   are not double counted due to the test that locpol[0][2] > locpol[0][1].
   only one "plane" is needed since in principal the sysem has a c_inf
   axis along any pair of colloids (3 points determine a plane...)
   These configurations represent all the unique 3-body configurations. */
  for ( locpol[0][1] = 3; locpol[0][1] < L/2; locpol[0][1]++ )
  for ( locpol[1][2] = 0; locpol[1][2] < L/2; locpol[1][2]++ )
  for ( locpol[0][2] = 0; locpol[0][2] < L/2; locpol[0][2]++ )
  if ( locpol[1][2] || locpol[2][2] || (locpol[0][2] > locpol[0][1]) ) {
    for (a=0;a<np;a++)
      vexc[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]] = 250.0;
    for (a=0;a<L3;a++)
      phi[a]=0.0;
/* set (npl,nmin)=(total # + simple ions,total # - simple ions) */
    npl=np*z + nimp;
    nmin=nimp;
/* intialize gamma, based on DH level phi's: */
    sumexp=0.0;
    sumexm=0.0;
    for (a=0;a<L3;a++) {
      sumexm += exp(-phi[a]-vexc[a]);
      sumexp += exp(phi[a]-vexc[a]);
    }
    gammap = npl/sumexp;
    gammam = nmin/sumexm;
    printf("initial gammap= %e %e %d\n",gammap,sumexp,L3);
    printf("initial gammam= %e %e\n",gammam,sumexm);

    action=0.0;
    for(a=0;a<np;a++) 
      action -= phi[locpol[0][a]*L*L + locpol[1][a]*L + locpol[2][a]];
    action *= z;

    for (a=0,nxyz=0;a<L3;a+=L*L) {
      nxp = pbp[a/(L*L)]*L*L;
      nxm = pbm[a/(L*L)]*L*L;
      for (b=0;b<L*L;b+=L) {
        nyp = pbp[b/L]*L;
        nym = pbm[b/L]*L;
        for (c=0;c<L;c++,nxyz++) {
          nzp = pbp[c];
          nzm = pbm[c];
	  action += -0.5*beta*phi[nxyz]*
		   (phi[nxp+b+c] + phi[a+nyp+c] + phi[a+b+nzp] +
		    phi[nxm+b+c] + phi[a+nym+c] + phi[a+b+nzm] - 6*phi[nxyz]) +
		    gammap*exp(phi[nxyz]-vexc[nxyz]) +
		    gammam*exp(-phi[nxyz]-vexc[nxyz]);
	}
      }
    }
    printf("initial action = %le\n",action);

    feqn = 2*FTOL;
    nit += 1.0;
    it = 0;
    while ( fabs(feqn) > FTOL ){
      it++; itave++;

      for (a=f=0,nxyz=0;a<L3;a+=L*L,f++) {
        nxp = pbp[a/(L*L)]*L*L;
        nxm = pbm[a/(L*L)]*L*L;
        for (b=0;b<L*L;b+=L) {
          nyp = pbp[b/L]*L;
          nym = pbm[b/L]*L;
          for (c=e=0;c<L;c++,nxyz++,e+=L*L) {
            nzp = pbp[c];
            nzm = pbm[c];
            t1 = -beta*(phi[nxp+b+c]+phi[a+nyp+c]+phi[a+b+nzp]
                 + phi[nxm+b+c]+phi[a+nym+c]+phi[a+b+nzm] - 6*phi[nxyz]);
            for (d=0;d<np;d++)
              if ( nxyz == locpol[0][d]*L*L+locpol[1][d]*L+locpol[2][d] )
                t1 -= z;
            t1 += gammap*exp(phi[nxyz]-vexc[nxyz]) -
		       gammam*exp(-phi[nxyz]-vexc[nxyz]);
            dphi[f+b+e][0] = eps*t1;
            dphi[f+b+e][1] = 0;
          }
	}
      }
      sign = 1;
      fourn(dphi,nn,&dim,&sign);
/* fourier accelerate */
      for (a=0;a<L3;a++) {
        dphi[a][0] *= psq[a];
        dphi[a][1] *= psq[a];
      }
      sign = -1;
      fourn(dphi,nn,&dim,&sign);

      for (nxyz=f=0;f<L;f++)
        for (b=0;b<L*L;b+=L)
          for (d=0;d<L3;d+=L*L)
            phi[nxyz]=phi[nxyz]-il3*dphi[f+b+d][0];
      /* same for these four lines
      for (f=0;f<L;f++)
        for (b=0;b<L;b++)
          for (d=0;d<L;d++)
            phi[f*L*L+b*L+d]=phi[f*L*L+b*L+d]-il3*dphi[f+b*L+d*L*L][0];
      */

      action=0.0;
      sumexp=0.0;
      sumexm=0.0;
      feqn=0.0;

      for (a=0,nxyz=0;a<L3;a+=L*L) {
        nxp = pbp[a/(L*L)]*L*L;
        nxm = pbm[a/(L*L)]*L*L;
        for (b=0;b<L*L;b+=L) {
          nyp = pbp[b/L]*L;
          nym = pbm[b/L]*L;
          for (c=0;c<L;c++,nxyz++) {
            nzp = pbp[c];
            nzm = pbm[c];
	    t1 = exp(phi[nxyz]-vexc[nxyz]);
	    t2 = exp(-phi[nxyz]-vexc[nxyz]);
            ep = beta*(phi[nxp+b+c]+phi[a+nyp+c]+phi[a+b+nzp]
                 + phi[nxm+b+c]+phi[a+nym+c]+phi[a+b+nzm]-6*phi[nxyz]);
            action -= 0.5*phi[nxyz]*ep;
            ef = -ep + gammap*t1 - gammam*t2;
            for (d=0;d<np;d++)
              if( nxyz == locpol[0][d]*L*L+locpol[1][d]*L+locpol[2][d] )
                ef -= z;
            feqn += ef*ef;
            sumexp += t1;
            sumexm += t2;
          }
	}
      }
      feqn=feqn/L3;
      feqn=sqrt(feqn);
      action += gammap*sumexp + gammam*sumexm;
      for (a=0;a<np;a++)
        action -= z*phi[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]];
      printf("action=%le accuracy=%le\n",action,feqn);
      gammap=npl/sumexp;
      gammam=nmin/sumexm;
      printf("gammap=%le gammam=%le\n",gammap,gammam);
    }

/*
C  compute Helmholtz free energy (meanfield)
*/
    amf=0.0;
    for (a=0;a<L*L*L;a++)
      amf += gammap*(0.5*phi[a]-1.0)*exp(phi[a]-vexc[a])
		-gammam*(0.5*phi[a]+1.0)*exp(-phi[a]-vexc[a]);
    for (a=0;a<np;a++)
      amf += 0.5*z*phi[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]];
    amf += npl*log(gammap)+nmin*log(gammam);
    printf("%d %d %d %d %20.16le\n",locpol[0][1],locpol[0][2],
			       locpol[1][2],locpol[2][2],amf);
    for (a=0;a<np;a++)
      vexc[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]] = 0.0;
    /* commenting out the next line apparently solves the compile problem... */
    exit(0);
  }
      fprintf(stderr,"%le\n",itave/nit);
}