Commit 63526985 authored by aarograh's avatar aarograh Committed by Henderson, Shane
Browse files

Add single precision version of GAMMAF 

Squash branch 'hotfix_gammf' into 'master'

* Add single precision version of GAMMAF

<!-- Replace this a detailed description of changes -->

<!-- Include a link to VERA development issues if appropriate, or delete this line -->

**Developer Checklist:**
- [x] Have you done a self-review after creating the merge request?
- [x] Have you filled in the Merge Request information (title, description) thoroughly?
- [x] Have you updated the relevant tickets (if this MR is linked to any VERA-dev tickets)?
- [x] Have you addressed all suggested feedback and commented on it to let the reviewer know? (Do not resolve discussions that the reviewer started)

**Reviewer Checklist:**
- [x] Have you confirmed all discussions were adequately addressed and resolved them all?
- [x] Does it conform to formatting guidelines?
- [x] Are there adequate and clear comments?
- [x] Is the design clean and sensible?
- [x] Are the changes optimal/efficient?
- [x] Were sufficient DBC checks added?
- [x] Are there unit tests? (if necessary)
- [x] Is the MR description clear, including a link to the VERA-Dev issue if appropriate?

**PSM Checklist**
- [x] Have you confirmed that all discussions were addressed, or that follow-on issues have been created for them?
- [x] Have you confirmed sufficient testing was conducted?
- [x] Does this impact other repositories?
- [x] Does the MR have an adequate description?
- [ ] If the MR has multiple commits, did you set the MR to squash merge?

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

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