Commit 8e57d5c1 authored by Henderson, Shane's avatar Henderson, Shane
Browse files

Merge branch 'hotfix_gammf' into 'master' 

Add single precision version of GAMMAF

See merge request https://code.ornl.gov/futility/Futility/-/merge_requests/411
parents 4111896c 63526985

src/ExtendedMath.f90

+45 −3
Original line number Diff line number Diff line
@@ -27,6 +27,11 @@ PUBLIC :: RotateQtrClockwise 
PUBLIC :: ThirdOrderBickley

PUBLIC :: RegularizedIncompleteBeta



INTERFACE GAMMAF

  MODULE PROCEDURE GAMMAF_single

  MODULE PROCEDURE GAMMAF_double

ENDINTERFACE GAMMAF



INTERFACE LeastCommonMultiple

  MODULE PROCEDURE LeastCommonMultiple_scalar

  MODULE PROCEDURE LeastCommonMultiple_A1
@@ -35,7 +40,44 @@ ENDINTERFACE LeastCommonMultiple 
CONTAINS

!

!-------------------------------------------------------------------------------

!> @brief Computes the GAMMA function

!> @brief Computes the single precision GAMMA function

!> @param z a single precision real to compute the value of gamma for

!> @returns g the value of the GAMMA function at z

!>

!> This routine uses the Lanczos approximation to compute the gamma function

!> The coefficients @c p are the coefficients used by the GNU Scientific

!> Library. It is a single precision function and should be accurate to 15

!> digits of precision. It agrees with the intrinsic function supplied by the

!> Intel compiler.

!>

  PURE RECURSIVE FUNCTION GAMMAF_single(z) RESULT(g)

    INTEGER(SIK),PARAMETER :: cg=7

    REAL(SSK),DIMENSION(0:8),PARAMETER :: p=(/0.99999999999980993_SSK, &

        676.5203681218851_SSK,-1259.1392167224028_SSK,771.32342877765313_SSK, &

        -176.61502916214059_SSK,12.507343278686905_SSK, &

        -0.13857109526572012_SSK,9.9843695780195716e-6_SSK, &

        1.5056327351493116e-7_SSK/)

    REAL(SSK),INTENT(IN) :: z

    REAL(SSK) :: g

    INTEGER(SIK) :: i

    REAL(SSK) :: t,w,x



    x=z

    IF(x < 0.5_SSK) THEN

      g=PI/(SIN(PI*x)*GAMMAF(1.0_SSK-x))

    ELSE

      x=x-1.0_SSK

      t=p(0)

      DO i=1,cg+1

        t=t+p(i)/(x+REAL(i,SSK))

      ENDDO

      w=x+REAL(cg,SSK)+0.5_SSK

      g=SQRT(2.0_SSK*PI)*(w**(x+0.5_SSK))*EXP(-w)*t

    ENDIF

  ENDFUNCTION GAMMAF_single

!

!-------------------------------------------------------------------------------

!> @brief Computes the double precision GAMMA function

!> @param z a double precision real to compute the value of gamma for

!> @returns g the value of the GAMMA function at z

!>
@@ -45,7 +87,7 @@ CONTAINS 
!> digits of precision. It agrees with the intrinsic function supplied by the

!> Intel compiler.

!>

  PURE RECURSIVE FUNCTION GAMMAF(z) RESULT(g)

  PURE RECURSIVE FUNCTION GAMMAF_double(z) RESULT(g)

    INTEGER(SIK),PARAMETER :: cg=7

    REAL(SDK),DIMENSION(0:8),PARAMETER :: p=(/0.99999999999980993_SDK, &

        676.5203681218851_SDK,-1259.1392167224028_SDK,771.32342877765313_SDK, &
@@ -69,7 +111,7 @@ CONTAINS 
      w=x+REAL(cg,SDK)+0.5_SDK

      g=SQRT(2.0_SDK*PI)*(w**(x+0.5_SDK))*EXP(-w)*t

    ENDIF

  ENDFUNCTION GAMMAF

  ENDFUNCTION GAMMAF_double

!

!-------------------------------------------------------------------------------

!> @brief Computes a Rational Fraction less than a tolerance