Date: Tue, 16 Sep 1997 14:20:23 +0200 (MEST) From: Christoph Kukulies <kuku@gilberto.physik.RWTH-Aachen.DE> To: freebsd-hackers@freefall.FreeBSD.org Subject: FPE problem Message-ID: <199709161220.OAA05839@gil.physik.rwth-aachen.de>
next in thread | raw e-mail | index | archive | help
I'm having a weird problem with a f2c program which is causing an FPE
for a reason I cannot figure out.
It happens in a complex expression in a fortran expression
in a scientific (high energy physics) calculation program.
y=(x-1.D0)/x
z=-log(1.D0-y)
z2=z*z
rsp=z*(1.D0+a1*z*(1.D0+a2*z*(1.D0+a3*z2*(1.D0+a4*z2*
1 (1.D0+a5*z2*(1.D0+a6*z2*(1.D0+a7*z2*(1.D0+a8*z2*(1.D0+a9*z2*
2 (1.D0+a10*z2))))))))))
3 +zeta2-log(x)*log(1.D0-x)+0.5D0*log(x)**2
The corresponding C code (with debugging printf):
y = 1. - *x;
z__ = -log(1. - y);
z2 = z__ * z__;
printf("z__=%lf\n"
"spint_1.a1=%lf\n"
"spint_1.a2=%lf\n"
"spint_1.a3=%lf\n"
"spint_1.a4=%lf\n"
"spint_1.a5=%lf\n"
"spint_1.a6=%lf\n"
"spint_1.a7=%lf\n"
"spint_1.a8=%lf\n"
"spint_1.a9=%lf\n"
"spint_1.a10=%lf\n"
"spint_1.zeta2=%lf\n"
"z2=%lf\n"
"*x=%lf\n",
z__,
spint_1.a1,
spint_1.a2,
spint_1.a3,
spint_1.a4,
spint_1.a5,
spint_1.a6,
spint_1.a7,
spint_1.a8,
spint_1.a9,
spint_1.a10,
spint_1.zeta2,
z2,
*x);
fflush(stdout);
ret_val = -z__ * (spint_1.a1 * z__ * (spint_1.a2 * z__ * (spint_1.a3 *
z2 * (spint_1.a4 * z2 * (spint_1.a5 * z2 * (spint_1.a6 * z2 *
(spint_1.a7 * z2 * (spint_1.a8 * z2 * (spint_1.a9 * z2 * (
spint_1.a10 * z2 + 1.) + 1.) + 1.) + 1.) + 1.) + 1.) + 1.) +
1.) + 1.) + 1.) + spint_1.zeta2 + z__ * log(1. - *x);
return ret_val;
gdb output:
----------
Nucleon PDFs : Ngroup = 3, Nset = 35
, for MRS Set (H) (L230-MSb) Structure Functions
----------------------------------------------------------------------------------------------------------------
z__=0.029266
spint_1.a1=-0.250000
spint_1.a2=-0.111111
spint_1.a3=-0.010000
spint_1.a4=-0.017007
spint_1.a5=-0.019444
spint_1.a6=-0.020661
spint_1.a7=-0.021417
spint_1.a8=-0.021949
spint_1.a9=-0.022349
spint_1.a10=-0.022664
spint_1.zeta2=1.644934
z2=0.000857
*x=0.971158
Program received signal SIGFPE, Arithmetic exception.
0x1b420 in rsp_ (x=0xefbfd784) at disent.c:5165
5165 ret_val = -z__ * (spint_1.a1 * z__ * (spint_1.a2 * z__ * (spint_1.a3 *
A small C program fed with the above values evaluates the expression
without problems.
main()
{
double log(),
z__=0.029266,
a1=-0.250000,
a2=-0.111111,
a3=-0.010000,
a4=-0.017007,
a5=-0.019444,
a6=-0.020661,
a7=-0.021417,
a8=-0.021949,
a9=-0.022349,
a10=-0.022664,
zeta2=1.644934,
x=0.971158,
z2,ret_val ,
y;
y = 1. - x;
z__ = -log(1. - y);
z2 = z__ * z__;
ret_val = -z__ * (a1 * z__ * (a2 * z__ * (a3 *
z2 * (a4 * z2 * (a5 * z2 * (a6 * z2 *
(a7 * z2 * (a8 * z2 * (a9 * z2 * (
a10 * z2 + 1.) + 1.) + 1.) + 1.) + 1.) + 1.) + 1.) +
1.) + 1.) + 1.) + zeta2 + z__ * log(1. - x);
printf("%lf\n",ret_val);
}
Any ideas?
--
Chris Christoph P. U. Kukulies kuku@gil.physik.rwth-aachen.de
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?199709161220.OAA05839>