[Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90

[tkoenig at gcc dot gnu dot org]

------- Comment #18 from tkoenig at gcc dot gnu dot org  2007-11-04 19:51 -------

> Looking at Pugh's paper, he shows coefficients for n = 10 and
> r = 10.900511.  This is the same as you are using for the double
> case.  However, you seem to be missing coefficient d10.

Good catch, thanks!

The main problem with the Lanczos approximation are alternating-sign
terms with nearly cancel each other, which leads to massive precision
loss.

For z=16.5 and r=10.900511, the terms in the sum are

1 6425.81075191694890236249
2 -19919.53511610527857556008
3 24595.63902224190678680316
4 -15425.21437829293790855445
5 5196.45802113903846475296
6 -916.60901318718765651283
7   76.62541745991659070114
8   -2.45417886377221794447
9    0.01907340042601936639
10   -0.00001080118945830201

and the sum is

 33.24621718250205049117

so I'd expect about log2(24000/33) ~ 9.5 bits of precision loss.
Not good.

Some rearrangement can help here, possibly.  OTOH, maybe using
straight Netlib code would be better.

Ouch.

Maybe it's better to fall back on the Netlib code.


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698