------- Comment #19 from dave at hiauly1 dot hia dot nrc dot ca 2007-11-04
21:30 -------
Subject: Re: FAIL: gfortran.dg/gamma_5.f90
> 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.
I don't think this can be avoided without a huge performance hit
(i.e., use long double arithmetic for the sum). Possibly, the
float versions would benefit from doing the sum with doubles.
That shouldn't hurt.
> Some rearrangement can help here, possibly. OTOH, maybe using
> straight Netlib code would be better.
On systems with lgamma, it's possible to use exp (lgamma (x)) for
tgamma (x). I checked and gamma_5 doesn't fail on HP-UX with this
implementation. I retained the transformation for negative arguments.
Float variants could easily be implemented using lgamma. I think
using the system lgamma implementation might be best when it's available.
Found an old fortran implementation on one of my machines! I believe
that I wrote an implementation using a rational approximation (Pade?)
about 35 years ago, but I couldn't find it.
Dave
--
J. David Anglin [EMAIL PROTECTED]
National Research Council of Canada (613) 990-0752 (FAX: 952-6602)
double precision function dgamma(z)
c
c Source for routine: J.F. Hart et al, Computer Approximations,
c Wiley,1968.
c This is a double precision routine for the gamma function of
c real positive arguments. Values for negative arguments (non-
c integral or zero) may be obtained from this program and
c standard functional relations for the gamma function.
c argument range method of computation
c z>14 gamma(z)=dexp(ln(gamma(z))
c 3<= z <=14 argument reduction z*gam(z)=gam(z+1)
c 2<= z <3 rational approx. w/o argument reduction
c 0< z <2 argument expansion z*gam(z)=gam(z+1)
c
implicit double precision (a-h,o-z)
c
double precision nlsqp
c
dimension pnum(11),qden(11),phnum(4),qhden(4)
c
data pnum,qden,phnum,qhden/
1-.2983543278574342138830437659d+6,
2-.2384953970018198872468734423d+6,
3-.1170494760121780688403854445d+6,
4-.3949445048301571936421824091d+5,
5-.1046699423827521405330650531d+5,
6-.2188218110071816359394795998d+4,
7-.3805112208641734657584922631d+3,
8-.5283123755635845383718978382d+2,
9-.6128571763704498306889428212d+1,
9-.502801805441681246736419875d00,
9-.3343060322330595274515660112d-1,
1-.2983543278574342138830438524d+6,
2-.1123558608748644911342306408d+6,
3+.5332716689118142157485686311d+5,
4+.8571160498907043851961147763d+4,
5-.473486597702821170655681977d+4,
6+.1960497612885585838997039621d+3,
7+.1257733367869888645966647426d+3,
8-.2053126153100672764513929067d+2,
9+.1d1,
90.d0,
90.d0,
1+.12398282342474941538685913d0,
2+.670827838343321349614617d0,
3+.6450730291289920251389d0,
4+.666629070402007526d-1,
1+.148779388109699298468156d+1,
2+.809952718948975574728214d+1,
3+.7996691123663644194772d+1,
4+.1d1/
c
az=z
pi=4.d0*datan(1.d0)
accum=1.d0
if (az.gt.14.d0) go to 44
if (az.lt.2.d0) go to 333
if (az.lt.3.d0) go to 33
do 22 i=1,12
az=az-1.d0
accum=accum*az
if (az.lt.3.d0) go to 33
22 continue
333 do 222 i=1,12
accum=accum*(1.d0/az)
az=az+1.d0
if (az.ge.2.d0) go to 33
222 continue
33 az=az-2.d0
anum=0.d0
aden=0.d0
do 11 i=0,10
ii=i+1
if (i.eq.0) power=1.d0
if (i.eq.0) go to 55
power=az**i
55 anum=anum+pnum(ii)*power
aden=aden+qden(ii)*power
11 continue
dgamma=accum*(anum/aden)
go to 999
44 anum=0.d0
aden=0.d0
do 1 i=0,3
ii=i+1
if (i.eq.0) power=1.d0
if (i.eq.0) go to 5
power=(1.d0/(az*az))**i
5 anum=anum+phnum(ii)*power
aden=aden+qhden(ii)*power
1 continue
riemn=anum/(aden*az)
nlsqp=dlog(dsqrt(2.d0*pi))
daleg=(az-0.5d0)*dlog(az)-az+nlsqp+riemn
dgamma=dexp(daleg)
999 return
end
--
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698
