Re: A problem with tgamma function

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

Brian Dessent wrote:

> That is a bug report against glibc.  You cannot draw any inferences
> between that and newlib/Cygwin because they are totally separate and
> unrelated code bases.

> That would be expected given Lev's indication of a fix in newlib

What I wished to stress was that the 'tgamma' function of newlib gives the
same results of the glibc-tgamma, for which a bug report has been sent.

And this regardless the "Lev's indication of a fix in newlib".

So it looks that the newlib-tgamma has similar problems to that of
glibc-tgamma.

Building the following test case with a recent CVS-GFortran-4.3.0
(that now uses tgamma), the result is

  Segment fault   with Cygwin-1.5.24-2

    fail          with Cygwin-snapshot 20070920

    fail          on GNU/Linux Kubuntu 7.04


Linking it with Cernlib (that has its own gamma implemented), it is
completed successfully.



Cheers,

   Angelo.




==========================================================================
!
! gfortran -O0 -funroll-loops -fomit-frame-pointer -fno-second-underscore
!          -fno-automatic gfc_gamma_test_c302m.F90 -o gfc_gamma_test_c302m
!
! Uncommenting the lines
!
!    !real :: gamma
!    !real(DP) :: dgamma
!    !external gamma,dgamma
!
! you must build in this way:
!
! gfortran -O0 -funroll-loops -fomit-frame-pointer -fno-second-underscore
!          -fno-automatic gfc_gamma_test_c302m.F90 `cernlib -gfc mathlib`
!          -o gfc_gamma_test_c302m
!
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 :: gamma
    !real(DP) :: dgamma
    !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 = gamma(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 = dgamma(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=gamma(-sngl(HALF))
       if (jf == 2) dr = dgamma(-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
====================================================================

--
Unsubscribe info:      http://cygwin.com/ml/#unsubscribe-simple
Problem reports:       http://cygwin.com/problems.html
Documentation:         http://cygwin.com/docs.html
FAQ:                   http://cygwin.com/faq/