Incorrect result for gamma functions of pure imaginary argument

[banana-bred]

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

a4ff2f054b640857c76dfb0d625d038dd94c9baa

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.

[jvdp1]

@banana-bred Thank you for reporting this bug. Could it be a special case when z%re is equal to 0 (I never use complex values)? Because gamma(x) and exp(log_gamma( x)) report (almost) equal values for cmplx(0.25, 0.25), cmplx(-0.25, 0.25), cmplx(0.25, -0.25).

[banana-bred]

I think so. The cases you mentioned seem to be handled appropriately, but z%re = 0 currently gets treated the same as z%re < 0, which is probably (?) not intended and is what seems to cause this issue.

[jvdp1]

Ok. Thank you for the clarification. Would you like to submit a PR with the appropriate fix and test, please?

Le lun. 14 août 2023 à 22:16, banana-bred @.***> a écrit :

I think so. The cases you mentioned seem to be handled appropriately, but z%re = 0 currently gets treated the same as z%re < 0, which is probably (?) not intended and is what seems to cause this issue.

— Reply to this email directly, view it on GitHub https://github.com/fortran-lang/stdlib/issues/729#issuecomment-1677997997, or unsubscribe https://github.com/notifications/unsubscribe-auth/AD5RO7EA7ZPICKQX62JDWJTXVKBTHANCNFSM6AAAAAA3QCLKCU . You are receiving this because you commented.Message ID: @.***>

[banana-bred]

sure

[jvdp1]

Fixed by #730