For future reference, Fortran has a gamma() function, which for positive whole numbers can be used to calculate factorials.
Factorials grow so rapidly that you will overflow a 32-bit INTEGER quite quickly.
Even typical DOUBLEPRECISION will overflow past 170! and Fortran is basically not required to indicate when overflow occurs (although compilers can commonly report it, depending on the compiler options used).
Of course factorials are often shown calculated via recursion, as it can be done with very little code, so you will often find recursive versions.
open calculate factorial in Fortran Playground
! GAMMA(X) computes Gamma of X. For positive whole number values of N the
! Gamma function can be used to calculate factorials, as (N-1)! ==
! GAMMA(REAL(N)). That is
!
! n! == gamma(real(n+1))
!
program demo_gamma
use, intrinsic :: iso_fortran_env, only: wp => real64, int64
implicit none
integer :: i, j
! gamma() is related to the factorial function
do i = 1, 171
! check value is not too big for default integer type
if (factorial(i) <= huge(0)) then
write (*, *) i, nint(factorial(i)), 'integer'
elseif (factorial(i) <= huge(0_int64)) then
write (*, *) i, nint(factorial(i),kind=int64), 'integer(kind=int64)'
else
write (*, *) i, factorial(i) , 'user factorial function'
write (*, *) i, product([(real(j, kind=wp), j=1, i)]), 'product'
write (*, *) i, gamma(real(i + 1, kind=wp)), 'gamma directly'
end if
end do
contains
function factorial(i) result(f)
integer, intent(in) :: i
real(kind=wp) :: f
if (i <= 0) then
write (*, '(*(g0))') '<ERROR> gamma(3) function value ', i, ' <= 0'
stop '<STOP> bad value in gamma function'
end if
f = gamma(real(i + 1,kind=wp))
end function factorial
end program demo_gamma