module QUAD_LNGAMMA_FUNCTION
use QUADRUPLE_PRECISION
implicit none

private
public :: Q_LNGM

contains

function Q_LNGM(X) result(B)

! Extended arithmetic calculation of the logarithm of the gamma
! function (N.B. gamma(x) = (x-1)! in factorial notation):
!  b = ln (gamma (x))
! where all quantities are in quadruple-precision.
! The result (b) may occupy the same location as the input value (x).
! x must not equal a negative integer; in this case a warning is given
! and no result is calculated.

! Algorithm:  The true Stirling's approximation (not the version derived by
! de Moivre which is usually quoted in textbooks) is used for x >= 16.
! For smaller arguments, a Lanczos-type approximation is used.

! Programmer: Alan Miller
!
! Latest Fortran 77 revision - 11 May 1988
! Fortran 90 version - 21 August 1997
! Version for the F compiler - 5 April 2000

type (QUAD), intent(in) :: X
type (QUAD)             :: B

! Local variables

logical                    :: LARGE
integer                    :: I, J
type (QUAD)                :: Z, TOTAL, TEMP, ZZ, TERM
real (kind=DP), parameter  :: ZERO = 0.0_DP, HALF = 0.5_DP, ONE = 1.0_DP
integer, parameter         :: G = 14

! Table of values of
!   a(2j) = B(2j).[2**(2j-1) - 1]/[2j.(2j-1).2**(2j-1)]
! where the B(2j)'s are the Bernouilli numbers.
! Stirling's approximation for log(gamma(x)) is then
!   log(gamma(x)) = 0.5*log(2.pi) + z.log(z) - z - sum[a(2j)/z**(2j-1)]
! where z = x - 0.5.

type (QUAD), parameter, dimension(16) :: ST_COEFF = (/  &
    QUAD(  0.4166666666666667e-01_DP, -0.4625929269271486e-17_DP),  &
    QUAD( -0.2430555555555556e-02_DP,  0.4722302795714642e-18_DP),  &
    QUAD(  0.7688492063492064e-03_DP, -0.8742455613057720e-19_DP),  &
    QUAD( -0.5905877976190476e-03_DP, -0.3820521941139396e-19_DP),  &
    QUAD(  0.8401067971380472e-03_DP, -0.2482348408384320e-19_DP),  &
    QUAD( -0.1916590625086719e-02_DP,  0.2283600073631585e-18_DP),  &
    QUAD(  0.6409473908253206e-02_DP, -0.7561615151693776e-18_DP),  &
    QUAD( -0.2954975178039152e-01_DP,  0.2751974392739155e-17_DP),  &
    QUAD(  0.1796430017910384e+00_DP,  0.5451863036096252e-16_DP),  &
    QUAD( -0.1392429561052216e+01_DP,  0.1922261747261675e-16_DP),  &
    QUAD(  0.1340285765318479e+02_DP, -0.1330122271048497e-14_DP),  &
    QUAD( -0.1568482659282295e+03_DP,  0.5486595470093716e-13_DP),  &
    QUAD(  0.2193103267973761e+04_DP, -0.1576457483073076e-12_DP),  &
    QUAD( -0.3610877098469368e+05_DP, -0.1128718158624283e-11_DP),  &
    QUAD(  0.6914722675633456e+06_DP, -0.9528010535925612e-11_DP),  &
    QUAD( -0.1523822153940742e+08_DP,  0.4711160925202246e-08_DP) /)

!     Lanczos-type coeffs. for g = 14

type (QUAD), parameter, dimension(16) :: A = (/  &
   QUAD(  0.1000000000000000E+01_DP,  0.4846757479531194E-21_DP),  &
   QUAD(  0.1069184011592089E+07_DP, -0.5763228334420908E-09_DP),  &
   QUAD( -0.4731034435392098E+07_DP,  0.1831263028328363E-09_DP),  &
   QUAD(  0.8831027007071320E+07_DP, -0.3562876555263793E-09_DP),  &
   QUAD( -0.9057605936125744E+07_DP, -0.8123034251454294E-09_DP),  &
   QUAD(  0.5575099541767454E+07_DP, -0.4791993277603296E-09_DP),  &
   QUAD( -0.2113664276594438E+07_DP,  0.2261448706472902E-09_DP),  &
   QUAD(  0.4882973818694724E+06_DP, -0.2584262485630417E-10_DP),  &
   QUAD( -0.6580271859490059E+05_DP,  0.6829382319982582E-11_DP),  &
   QUAD(  0.4754359892166024E+04_DP, -0.6321699947095538E-12_DP),  &
   QUAD( -0.1588514130362740E+03_DP,  0.1630647068074753E-13_DP),  &
   QUAD(  0.1878850329497677E+01_DP,  0.1651002503310329E-15_DP),  &
   QUAD( -0.4589872822043718E-02_DP,  0.3894567587921014E-18_DP),  &
   QUAD(  0.5931363621107231E-06_DP, -0.1388883823219028E-21_DP),  &
   QUAD(  0.1024789274886151E-11_DP,  0.2637600537816092E-27_DP),  &
   QUAD( -0.3355721053959592E-12_DP,  0.4593372363497046E-29_DP) /)

! Check negative arguments.

if (X%HI <= ZERO) then
  write(unit=*, fmt=*) " *** Negative argument in q_lngm ***"
  return
end if

if (X%HI == ONE .or. X%HI == 2.0_DP) then
  if (X%LO == ZERO) then
    B%HI = ZERO
    B%LO = ZERO
    return
  end if
end if

if (X%HI < 16.0_DP) then

!  Fit a Lanczos-type approximation if x%hi < 16.
!  i.e. ln(gamma(x)) = lnsqrt(2.pi) + (x-.5)ln(x+g-.5) - (x+g-.5) + ln(A(x))
!  where A(x) = a0 + a1/x + a2/(x+1) + ... + ak/(x+g).
!  g = 14 here to achieve about 29 significant digits accuracy, except in the
!  vicinity of 1 and 2, where ln(gamma) = 0, where the absolute accuracy is
!  about 30 decimal digits.

  Z = A(1)
  do I = 2, G + 2
    Z = Z + A(I) / (X + (I-2))
  end do
  TEMP = X - HALF
  ZZ = TEMP + G
  B = LNSQRT2PI + TEMP*log(ZZ) - ZZ + log(Z)

  return
end if

! Stirling's approximation.

TOTAL = LNSQRT2PI
Z = X - HALF
TOTAL = TOTAL + Z * log(Z) - Z
ZZ = Z * Z
TERM = Z
LARGE = .true.
do J = 1, 16
  if (LARGE) then
    TEMP = ST_COEFF(J) / TERM
    TOTAL = TOTAL - TEMP
    TERM = TERM * ZZ
    if (abs(TEMP%HI) < abs(TOTAL%LO)) then
      LARGE = .false.
    end if
  else
    TEMP%LO = ST_COEFF(J)%HI / TERM%HI
    TOTAL%LO = TOTAL%LO - TEMP%LO
    if (abs(TEMP%LO) < 1.0e-30) then
      cycle
    end if
    TERM%HI = TERM%HI * ZZ%HI
  end if
end do

! The smallest bit of sum%hi & the largest bit of sum%lo may be out
! of alignment, so add zero to re-align.

B = TOTAL + ZERO

return
end function Q_LNGM

end module QUAD_LNGAMMA_FUNCTION



program TEST_Q_LNGM
! Test the quadruple-precision lngamma function using
!    lngamma(m+x) - lngamma(n+x) = log[ (m+x-1).(m+x-2) .. (n+x) ]

use QUADRUPLE_PRECISION
use QUAD_LNGAMMA_FUNCTION
implicit none

integer                :: M, N, I
type (QUAD)            :: X, LNGMM, LNGMN, DIFF, LOG_PROD, ERROR, TEMP
type (QUAD), parameter :: QONE = QUAD( 1.0_DP, 0.0_DP)

do
  write(unit=*, fmt="(a)", advance="NO") " Enter 2 different positive integers: "
  read(unit=*, fmt=*) M, N
  if (M == N) then
    write(unit=*, fmt=*) "The integers must be different - try again!"
    cycle
  end if
  if (N > M) then
    I = M
    M = N
    N = I
  end if
  call random_number(X%HI)
  X%LO = 0.0_DP
  LNGMM = Q_LNGM(X + real(M,kind=DP))
  write(unit=*, fmt="(' lngamma(', f20.15, ') = ', f22.16)")  &
                X%HI+real(M,kind=DP), LNGMM%HI
  LNGMN = Q_LNGM(X + real(N,kind=DP))
  write(unit=*, fmt="(' lngamma(', f20.15, ') = ', f22.16)")  &
                X%HI+real(N,kind=DP), LNGMN%HI
  DIFF = LNGMM - LNGMN
  TEMP = X + real(M-1,kind=DP)
  LOG_PROD = log(TEMP)
  do I = 1, M-N-1
    TEMP = TEMP - QONE
    LOG_PROD = LOG_PROD + log(TEMP)
  end do
  ERROR = DIFF - LOG_PROD
  write(unit=*, fmt="(' Diff =', f20.15, '  log_prod =', f20.15, '  Error =', g11.3)")  &
                DIFF%HI, LOG_PROD%HI, ERROR%HI
end do

stop
end program TEST_Q_LNGM