Description
I don't know if the stdlib project would concern itself with standard intrinsic functions. Working with gfortran (10.3.1 on a Fedora Unix system) I found surprising behavior of the gamma intrinsic; in many cases it does not return the exact integer-valued result for small integer-valued argument, and in some cases even anint(gamma(z+1)) returns a wrong result although z! is exactly representable in the real kind used in my code. I provide a demonstration and testing code below. A much broader view of accuracy of implementations of intrinsic functions is provided in this working paper: [Vincenzo Innocente and Paul Zimmermann: Accuracy of Mathematical Functions in Single, Double, Double Extended, and Quadruple Precision. Web link: https://members.loria.fr/PZimmermann/papers/accuracy.pdf].
Test and demonstration code:
program TestFactorial
use iso_fortran_env, only : wp => real32
implicit none
! integer, parameter :: wp=10 ! Try this too
logical :: b0
integer :: m0, m1
write (*,'(a)') 'Testing the factorial function ...'
call Factorial_Limits (m0, m1)
! m0 and m1 are the largest values such that using real (kind=wp):
! For 0<=k<m0 Factorial(k) can be represented exactly.
! For 0<=k<m1 Factorial(k) is below the overflow limit.
write (*,'(a,2i5)') 'm0, m1:', m0, m1
b0 = .false.
! b0 has the value (we have found a discrepancy between two ways of
! computing the factorial).
block
integer :: i
real (kind=wp) :: t0, t1
do i = 0, m0-1
if (i==0) then
t0 = 1
else
t0 = i*t0
end if
t1 = gamma(i+1.0_wp)
! t0 = i! computed using the recursion.
! t1 = i! computed via the gamma intrinsic function.
if (t1/=t0) then
b0 = .true.
write (*,'(a,1x,i3,2(1x,1pg10.2))') 'error:', i, &
spacing(t0), (t1-t0)/spacing(t0)
end if
end do
end block
if (.not.b0) then
write (*,'(a)') '... PASS'
end if
stop
contains
subroutine Factorial_Limits (m0, m1)
integer, intent (out) :: m0, m1
! Determine the limits for computation of the Factorial function
! using real (kind=wp).
! m0: for 0<=k<m0, Factorial(k) can be represented exactly.
! m1: for 0<=k<m1, Factorial(k) is below the overflow limit.
integer :: m0loc, i, j
real (kind=wp) :: t0, t1
logical :: b0
i = 0 ; t0 = 1 ; t1 = 1 ; b0 = .true.
m0loc = huge(0)
! Invariant: 0 <= i.
! t0 = Factorial(i) (subject to roundoff for large i).
! t1 = Factorial(i)/2^k where 2^k is the largest power of two that
! leaves an (odd) integer result.
! b0 = ((i+1)*t0 will not overflow).
! m0loc is the smallest nonnegative integer value i for which we
! have deterimined that Factorial(i) cannot be represented exactly.
do while (b0)
i = i+1
t0 = i*t0
j = i
do while (modulo(j,2).eq.0)
j = j/2
end do
t1 = j*t1
if (1<spacing(t1)) then
! Factorial(i) cannot be represented exactly.
m0loc = min(i,m0loc)
end if
b0 = log(t0)+log(real(i+1,kind=wp)) < log(huge(t0))
end do
m0 = min(i+1,m0loc)
m1 = i+1
return
end subroutine Factorial_Limits
end program(The version posted here is for 32-bit reals; edit the definition of wp at the top to obtain another real format.)
The output for real kind real32 on my gfortran 10.3.1 system:
bash-5.0$ gfortran TestFactorial.f90
bash-5.0$ ./a.out
Testing the factorial function ...
m0, m1: 14 35
error: 2 2.38E-07 1.0
error: 3 4.77E-07 1.0
error: 5 7.63E-06 1.0
error: 9 3.12E-02 1.0
error: 11 4.0 1.0
error: 12 32. 1.0
error: 13 5.12E+02 2.0
In the output, m0 and m1 are the largest values such that using real (kind=wp), for 0<=k<m0 Factorial(k) can be represented exactly and for 0<=k<m1 Factorial(k) is below the overflow limit. The following lines show values of k in the range 0<=k<m0 for which the gamma function intrinsic computation of k! is not exact; the second numerical value in the line shows the local number spacing at Factorial(k) and the third numerical value shows the error in units of that spacing.
In particular I draw attention to the results for 11!, 12! and 13!, all of which are wrong by an integer offset although the values are exactly representable in the real32 kind.
The output for real kinds real64 and real128 shows the same issue, as does the output for "kind=10" that gives 80-bit reals on my system.