Description
Description
The gamma function implemented in cgamma_csp and cgamma_cdp in module stdlib_specialfunctions_gamma produce the wrong output when the argument is a pure imaginary. They return the complex conjugate of the output instead.
Expected Behaviour
An example where the stdlib implementation of gamma returns the complex conjugate of the expected result. The expected result here is calculated (with worse accuracy) by the user-defined function gammaWeierstrass, which implements
https://en.wikipedia.org/wiki/Gamma_function#19th_century:_Gauss,_Weierstrass_and_Legendre ,
and by exp(log_gamma(z)). The latter two agree up to some low accuracy, but are both different than gamma(z).
This is not expected, as exp(log_gamma(z)) should be the same as gamma(z).
test.f90:
program test
use iso_fortran_env, only: stdout => output_unit
use stdlib_specialfunctions_gamma, only: gamma, log_gamma
implicit none
integer, parameter :: sp = selected_real_kind(6)
integer, parameter :: dp = selected_real_kind(15)
complex(sp), parameter :: a = cmplx(0, 1, kind = sp)
complex(dp), parameter :: b = cmplx(0, 1, kind = dp)
write(stdout,'(9X," -- Stdlib Γ(z) -- ")',advance='no')
write(stdout,'(10X,"|"," -- Weierstrass Γ(z) -- ")', advance='no')
write(stdout,'(8X,"|"," -- exp(ln(Γ(z))) -- ")')
write(stdout,'(A,X,2E15.7," | ",2E15.7," | ",2E15.7)') "Γ(i): ", gamma(a), gammaWeierstrass(b), exp(log_gamma(b))
write(stdout,'(A,X,2E15.7," | ",2E15.7," | ",2E15.7)') "Γ(-i):", gamma(-a), gammaWeierstrass(-b), exp(log_gamma(-b))
write(stdout,'(A,X,2E15.7," | ",2E15.7," | ",2E15.7)') "Γ(i): ", gamma(b), gammaWeierstrass(b), exp(log_gamma(b))
write(stdout,'(A,X,2E15.7," | ",2E15.7," | ",2E15.7)') "Γ(-i):", gamma(-b), gammaWeierstrass(-b), exp(log_gamma(-b))
write(stdout,*)
contains
impure elemental function gammaWeierstrass(z) result(res)
implicit none
complex(dp), intent(in) :: z
complex(dp) :: res
real(dp), parameter :: EulerMascheroni = 0.57721566490153286060651209008240243104215933593992_dp
complex(dp), parameter :: one = cmplx(1, kind = dp)
integer :: k
complex(dp) :: kc
res = exp(-EulerMascheroni * z) / z
do k = 1, 100000
kc = cmplx(k, kind = dp)
res = res * exp( z / kc ) / ( 1 + z / kc )
end do
end function gammaWeierstrass
end program test
output:
-- Stdlib Γ(z) -- | -- Weierstrass Γ(z) -- | -- exp(ln(Γ(z))) --
Γ(i): -0.1549498E+00 0.4980157E+00 | -0.1549506E+00 -0.4980182E+00 | -0.1549498E+00 -0.4980157E+00
Γ(-i): -0.1549498E+00 -0.4980157E+00 | -0.1549506E+00 0.4980182E+00 | -0.1549498E+00 0.4980157E+00
Γ(i): -0.1549498E+00 0.4980157E+00 | -0.1549506E+00 -0.4980182E+00 | -0.1549498E+00 -0.4980157E+00
Γ(-i): -0.1549498E+00 -0.4980157E+00 | -0.1549506E+00 0.4980182E+00 | -0.1549498E+00 0.4980157E+00
Compiler: gfortran 13.2.1
Version of stdlib
Platform and Architecture
Artix Linux x86_64
Additional Information
It seems like this comes from the imaginary axis not being included in the right complex half-plane ( z % re > 0 ).
This causes the argument to be complex conjugated, which then conjugates (*) the output because Γ(z*) = Γ(z)*.
Suggested fix: include the imaginary axis in the right complex half-plane in gamma_cdp and gamma_csp.
diff --git a/src/stdlib_specialfunctions_gamma.f90 b/src/stdlib_specialfunctions_gamma.f90
index d808309..68f2ec5 100644
--- a/src/stdlib_specialfunctions_gamma.f90
+++ b/src/stdlib_specialfunctions_gamma.f90
@@ -341,14 +341,14 @@ contains
end if
- if(z % re > zero_k1) then
+ if(z % re < zero_k1) then
- y = z - one
+ x = cmplx(abs(z % re), - z % im, kind = sp)
+ y = x - one
else
- x = cmplx(abs(z % re), - z % im, kind = sp)
- y = x - one
+ y = z - one
end if
@@ -425,14 +425,14 @@ contains
end if
- if(z % re > zero_k1) then
+ if(z % re < zero_k1) then
- y = z - one
+ x = cmplx(abs(z % re), - z % im, kind = dp)
+ y = x - one
else
- x = cmplx(abs(z % re), - z % im, kind = dp)
- y = x - one
+ y = z - one
end if
This format of defining the left half-plane with a strict inequality, i.e. if(z % re < zero_k1), mirrors the behavior in l_gamma_cdp and l_gamma_csp, which is why log_gamma seems to not have this issue.