------- 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