🚀 go-pugleaf

The following code causes gcc to get a fatal error if compiled with
the "-O" option.  This is on a Sun4/110 running SunOS 4.0.

Tom Quinn                 Canadian Institute for Theoretical Astrophysics
trq@moose.cita.toronto.edu
SOON TO BE trq@moose.cita.utoronto.ca
UUCP   - decvax!utgpu!moose!trq
BITNET - quinn@utorphys.bitnet
ARPA   - trq%moose.cita.toronto.edu@relay.cs.net

The compile:
gcc -g -v -O  -sun4 -c  nsat.c
gcc version 1.28
 /usr/local/lib/gcc-cpp -v -undef -D__GNU__ -D__GNUC__ -Dsparc -Dsun -Dunix -D__OPTIMIZE__ nsat.c /tmp/cca00727.cpp
GNU CPP version 1.28
 /usr/local/lib/gcc-cc1 /tmp/cca00727.cpp -quiet -dumpbase nsat.c -g -O -version -o /tmp/cca00727.s
GNU C version 1.28 (sparc) compiled by GNU C version 1.28.

nsat.c:80: warning: floating point number exceeds range of `double'
gcc: Program cc1 got fatal signal 6.

The code:
------------------------------------------------------------------------
extern double fabs(), floor();
extern double sqrt();
extern double  pow();
extern double sin(), cos();
extern double drem();	 
struct sat {
    double a;			 
    double e;			 
    double pomega;
    double mu;			 
    double omega;		 
    double phi;
    double p0;			 
    double yr, yi;		 
    int n;			 
    int index;			 
    int current;		 
};
static double dt;
main()
{
    struct sat *planet, test;	 
    int nplanet;
    int nsteps, i, print;
    double tmax, t;
    double p, deltay, gamma;
    double cphi, sphi;
    struct sat *sptr, *encptr, *nextpl();
    char *malloc();
    scanf("%d", &nplanet);
    if((planet = (struct sat *)malloc((unsigned) nplanet*sizeof(struct sat)))
       == 0 ) {
	exit(-1);
    }
    for(sptr = planet, i = 0; sptr < planet + nplanet; sptr++, i++) {
	scanf("%lf %lf %lf %lf", &(sptr->a), &(sptr->e), &(sptr->pomega),
	      &(sptr->mu));
	sptr->index = i;	 
	sptr->current = 0;	 
    }
    sptr = &test;
    scanf("%lf %lf %lf %lf", &(sptr->a), &(sptr->e), &(sptr->pomega),
	  &(sptr->mu));
    scanf("%d %d %lf %lf %lf", &nsteps, &print, &tmax, &t, &(sptr->phi));
    test.yr = test.e * cos(test.pomega);
    test.yi = test.e * sin(test.pomega);
    for(sptr = planet; sptr < planet + nplanet; sptr++) {
	sptr->p0 = (test.a - sptr->a)/sptr->a;
	sptr->n = 0;
	sptr->omega = pow(sptr->a, -1.5);
	sptr->phi = sptr->omega*t;
	sptr->phi -= (2.0*3.14 ) *(floor(sptr->phi/(2.0*3.14 ) ));
	sptr->yr = sptr->e * cos(sptr->pomega);
	sptr->yi = sptr->e * sin(sptr->pomega);
    }
    encptr = planet;
    for(i = 0; i < nsteps && t < tmax; i++) {
	test.omega = pow(test.a , -1.5);
	if(i % print == 0 && i != 0) {
	    double pphi = encptr->omega*3.14 /fabs(encptr->omega - test.omega);
	    pphi += test.phi;
	    cphi = cos(pphi);
	    sphi = sin(pphi);
	    printf("%g %g %g %.18g %g %d\n", test.a,
		   test.yr*cphi + test.yi*sphi,
		   -test.yr*sphi + test.yi*cphi, test.phi, t, encptr->index);
	}
	if((encptr = nextpl(planet, nplanet, test)) == 0 ) {  
	    break;
	}
	encptr->n++;
	test.phi += test.omega * dt;
	test.phi = drem(test.phi, (2.0*3.14 ) );
	if(test.phi < 0.0) test.phi = (2.0*3.14 )  + test.phi;
	test.phi -= (2.0*3.14 ) *(floor(test.phi/(2.0*3.14 ) ));
	for(sptr = planet; sptr < planet + nplanet; sptr++) {
	    sptr->phi += sptr->omega * dt;
	    sptr->phi = drem(sptr->phi, (2.0*3.14 ) );
	    if(sptr->phi < 0.0) sptr->phi = (2.0*3.14 )  + sptr->phi;
	    if(dt > 0.0) sptr->current = 0;
	}
	encptr->current++;
	t += dt;
	if(test.a > encptr->a)
	    deltay = encptr->mu * -2.23  /(encptr->p0 * encptr->p0);
	else
	    deltay = (-encptr->mu) * -2.23  /(encptr->p0 * encptr->p0);
	p = (test.a - encptr->a)/encptr->a;
	gamma = 0.75*((p)*(p))  - ((encptr->yr - test.yr)*(encptr->yr - test.yr)) 
	  - ((encptr->yi - test.yi)*(encptr->yi - test.yi)) ;
	test.yr += deltay*sin(test.phi);
	test.yi -= deltay*cos(test.phi);
	test.e = sqrt(((test.yr)*(test.yr))  + ((test.yi)*(test.yi)) );
	p = sqrt((4.0/3.0) * (gamma + ((test.yr - encptr->yr)*(test.yr - encptr->yr)) 
				    + ((test.yi - encptr->yi)*(test.yi - encptr->yi)) ));
	if(test.a > encptr->a)
	    test.a = encptr->a * (1.0 + p);
	else
	    test.a = encptr->a * (1.0 - p);
	if((test.a > encptr->a
	    && (1.0 - test.e)*test.a <= (1.0 + encptr->e)*encptr->a) 
	   || (test.a < encptr->a
	       && (1.0 + test.e)*test.a >= (1.0 - encptr->e)*encptr->a)) {
	    break;
	}
    }	
}