Re: Problems with GAMMA functions?

[Angelo Graziosi <Angelo dot Graziosi at roma1 dot infn dot it>]

On Sat, 13 Oct 2007, Thomas Koenig wrote:

> I have submitted this inaccuracy as a bug to the glibc people
> as http://sources.redhat.com/bugzilla/show_bug.cgi?id=5159 .
> Gcc, g++ and gfortran on glibc-based systems will all
> benefit if a solution is found there.
>
> If this is important for you, I would suggest that you contribute
> on this bug report.  In the meantime, you are probably better
> off using the external gamma function from CERNLIB.
>

I wrote: http://gcc.gnu.org/ml/fortran/2007-10/msg00197.html


Obviously the following passes too!

   Angelo.

============================================================================
/*
   dgammf_fortran.c
*/

double dgammf_(double *x)
{
   static const double PI = 3.14159265358979324;

   static const double C[] =
   {
      3.65738772508338244E0,
      1.95754345666126827E0,
      0.33829711382616039E0,
      0.04208951276557549E0,
      0.00428765048212909E0,
      0.00036521216929462E0,
      0.00002740064222642E0,
      0.00000181240233365E0,
      0.00000010965775866E0,
      0.00000000598718405E0,
      0.00000000030769081E0,
      0.00000000001431793E0,
      0.00000000000065109E0,
      0.00000000000002596E0,
      0.00000000000000111E0,
      0.00000000000000004E0
   };

   static int i;
   static double u,f,h,alfa,b0,b1,b2;

   u = (*x);

   if (u <= 0)
   {
      /* For non positive integer there are poles, return 0 is a choice */
      if (u == (int)(*x))
      {

	 h = 0;

	 goto LABEL_9;
      }
      else
      {
	 u = 1-(*x);
      }
   }

   f = 1;

   if (u < 3)
   {
      i = (int)(4-u);
      while (i--)
      {
	 f /= u;
	 u += 1;
      }
   }
   else
   {
      i = (int)(u-3);
      while (i--)
      {
	 u -= 1;
	 f *= u;
      }
   }

   h = u+u-7;
   alfa = h+h;
   b1 = 0;
   b2 = 0;

   for (i = 15; i >= 0; i--)
   {
      b0 = C[i]+alfa*b1-b2;
      b2 = b1;
      b1 = b0;
   }

   h = f*(b0-h*b2);

   if ((*x) < 0) h = PI/(sin(PI*(*x))*h);

   LABEL_9:
      return h;
}

----------------------------------------------------------------------------
/*
   gammf_fortran.c
*/
#define SROUND(D) ((D)+((D)-((float)(D))))

float gammf_(float *x)
{
   double dgammf_(double*),y = (double)(*x);

   return (SROUND(dgammf_(&y)));
}

----------------------------------------------------------------------------
!
!    test_c302m.F90
!
!
! gcc -c dgammf_fortran.c
! gcc -c gammf_fortran.c
! gfortran -c test_c302m.F90
! gfortran -o test test_c302m.o dgammf_fortran.o gammf_fortran.o
!
! ./test
!
program gfc_gamma_test_c302m
  implicit none
  integer, parameter :: LOUT = 10
  integer :: ntest = 0
  call c302m()
contains
  subroutine c302m
    implicit none
    logical :: ltest=.true.,ltest1=.true.,ltest2=.true.
    call header('C302',0)
    call c302d(ltest1)
    ltest=ltest .and. ltest1
    call test_result('C302',ltest)
    call pagend('C302')
  end subroutine c302m
  !
  subroutine c302d(ltest1)
    implicit none
    ! implicit double precision (A-H,O-Z)
    logical, intent(inout) :: ltest1
    !
    integer, parameter :: DP = kind(1.D0)
    real(DP), parameter :: HALF = 5D-1, PI=3.14159265358979324D0

    ! Set maximum error allowed for test to be considered successful
    real(DP), parameter :: TOL(2) = (/1D-6,5D-14/)
    !
    integer :: jf,n,nf
    character(len=6) :: tfunc(2) = (/'GAMMA ','DGAMMA'/)
    real(DP) :: dr,er,etol,errmax,rerrmax,t,x,c(0:20)
    real :: gammf
    real(DP) :: dgammf
    !external gamma,dgamma

    LTEST1=.true.

    c(0)=1
    do  n = 1,20
       c(n)=(2*n-1)*c(n-1)/2
    enddo

    nf=2

    do jf = 1,nf
       errmax= 0.0D0
       rerrmax= 0.0E0
       write(LOUT,'(/10X,''Test of C302 '',A)') tfunc(jf)
       write(LOUT,'(/9X,''X  '',7X,''Exact'',25X,''Calculated'',&
            &14X,''Rel. Error'')')

       do n =  1,20
          x = n+HALF
          t = c(abs(n))*sqrt(PI)

          if (jf == 1) then
             dr = gammf(sngl(x))
             if (dr /= 0) er = abs(sngl(dr-t)/sngl(dr))
             write(LOUT,'(1X,F10.1,2E27.9,5X,E10.1)')&
                  sngl(x),sngl(t),sngl(dr),sngl(er)
          endif
          if (jf == 2) then
             dr = dgammf(x)
             if (dr /= 0) er = abs((dr-t)/dr)
             write(LOUT,'(1X,F10.1,2E27.18,5X,E10.1)') x,t,dr,er
          endif

          errmax= max(errmax,er)
       enddo

       etol = TOL(jf)

       write(LOUT,'(/''Largest relative error was'',E10.1)') errmax
       ltest1 = ltest1.and.(errmax <= etol)

       write(LOUT,'(/''Testing error messages:'')')

       if (jf == 1) dr=gammf(-sngl(HALF))
       if (jf == 2) dr = dgammf(-HALF)
    enddo
    return
  end subroutine c302d
  !
  subroutine header(code,mode)
    ! This routine prints a page header with a testing routine name
    ! message.
    implicit none
    character(len=*), intent(in) :: code
    integer, intent(in) :: mode
    ntest = ntest+1
    if(mode == 1) &
         write(*,'(" Test#",I3," ( ",A," ):     started")') ntest,code
    write(LOUT,'("1",25X,30("*")/26X,"**   Testing Routine ",A,3X,"**"/&
         &26X,30("*"))') code
  end subroutine header
  !
  subroutine pagend(code)
    ! This subroutine prints a page end message
    implicit none
    character(len=*), intent(in) :: code
    write(LOUT,'(/26X,30("*")/26X,"**   End of Test of  ",A,3X,"**"/&
         &26X,30("*"))') code
  end subroutine pagend
  !
  subroutine test_result(code,test)
    character(len=*),intent(in) :: code
    logical, intent(in) :: test
    character(len=*), parameter ::&
         fmt1 = '(" Test#",I3," ( ",A," ):     completed successfully")',&
         fmt2 = '(" Test#",I3," ( ",A," ): *** failed ***")'
    if (test) then
       write(*,fmt1) ntest,code
       if (LOUT /= 6) write(LOUT,fmt1) ntest,code
    else
       write(*,fmt2) ntest,code
       if (LOUT /= 6) write(LOUT,fmt2) ntest,code
    endif
  end subroutine test_result
end program

============================================================================