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

[dave at hiauly1 dot hia dot nrc dot ca]

------- 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                                  dave.anglin@nrc-cnrc.gc.ca
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