Compiling with warning maybe-uninitialized show a potential bug in gamma functions

[cyrilgandon]

Description

In

stdlib/src/stdlib_specialfunctions_gamma.fypp

Line 580 in 68524b3

${t2}$ :: res, p_lim, a, b, g, c, d, y, ss

the value ss is declared:

and then use here:

stdlib/src/stdlib_specialfunctions_gamma.fypp

Line 673 in 68524b3

m = nint(ss)

The ss variable is never initalized, and generates 2 warnings with

fpm build --verbose --flag '-Wmaybe-uninitialized'

..\src\temp\stdlib_specialfunctions_gamma.f90:1424:24:

1424 | m = nint(ss)
| ^
Warning: 'ss' may be used uninitialized [-Wmaybe-uninitialized]
..\src\temp\stdlib_specialfunctions_gamma.f90:1331:52:

1331 | real(qp) :: res, p_lim, a, b, g, c, d, y, ss
| ^
note: 'ss' was declared here
..\src\temp\stdlib_specialfunctions_gamma.f90:1290:24:

1290 | m = nint(ss)
| ^
Warning: 'ss' may be used uninitialized [-Wmaybe-uninitialized]
..\src\temp\stdlib_specialfunctions_gamma.f90:1197:52:

1197 | real(dp) :: res, p_lim, a, b, g, c, d, y, ss
| ^
note: 'ss' was declared here
[ 12%] stdlib_specialfunctions_gamma. done.

Expected Behaviour

Either the value to be initialized, of the chunk of code deleted.

Version of stdlib

68524b3

Platform and Architecture

Windows

Additional Information

No response

[perazz]

Thanks @cyrilgandon, from a quick look, it really looks like a bug.
The mentioned reference can be found at this link ("Algorithm 2").

From a quick look, we have variable names:

paper | stdlib 
    t | a
    c | c
    d | d
    s | b
    l | n

Not only ss, but also m should be unused and it seems like both should be replaced by p

[perazz]

From what I see, that chunk of code works for real x<0 and real p, there is no test program for that. Instead, the real x<0, integer p implementation seems correct:

stdlib/src/stdlib_specialfunctions_gamma.fypp

Lines 810 to 832 in 68524b3

	a = -x
	c = one / a
	d = p - 1
	b = c * (a - d)
	n = 1
	do
	c = d * (d - one) / (a * a)
	d = d - 2
	y = c * ( a - d)
	b = b + y
	n = n + 1
	if(int(n, ${k1}$) > (p - two) / two .or. y < b * tol_${k2}$) exit
	end do
	if(y >= b * tol_${k2}$ .and. mod(p, two) /= zero_k1) &
	b = b + d * c / a
	g = ((-1) ** p * exp(-a + l_gamma(p, one) - (p - 1) * log(a)) &
	+ b ) / a