046. Use GAMMA to compute factorials#

topic: Math intrinsic functions

Why does Fortran not have a factorial function? It does, as for positive whole values of \(n\), the gamma function simplifies to the factorial function for \(n-1\). That is,

n! == gamma(n+1)
factorial.f90 | | 0 | Godbolt Compiler Explorer logo | Fortran logo#
program factorial
   use, intrinsic :: iso_fortran_env, only: wp => real64

   integer, parameter :: n(*) = [0, 1, 5, 11, 170]
   integer            :: j

   do j = 1, size(n)
      print '(*(g0,1x))', 'factorial of', n(j), ' is ', &
         product([(real(i, kind=wp), i=1, n(j))]), &
         ' or ', &
         gamma(real(n(j) + 1, kind=wp))
   end do
   print *, huge(0._wp)

end program factorial
Output1#
factorial of 0  is  1.0000000000000000  or  1.0000000000000000
factorial of 1  is  1.0000000000000000  or  1.0000000000000000
factorial of 5  is  120.00000000000000  or  120.00000000000000
factorial of 11  is  39916800.000000000  or  39916800.000000000
factorial of 170  is  0.72574156153079940E+307  or  0.72574156153079990E+307
   1.7976931348623157E+308

Intrinsic function gamma was added in the F2008 standard.

We can see that \(170!\) is close to the real64 limit. \(171!\) would be too large.

Also used above:

  • implied do loop for computing the array passed to product

  • assumed-size array (n(*)) – the size is “assumed” from the size of the array literal on the right-hand side of the equals sign

  • g0 edit descriptor and * infinite repeat count – see Tip 041

Note

Thanks to @urbanjost for this tip and code.

1

Compiled using GNU Fortran (Ubuntu 11.3.0-1ubuntu1~22.04) 11.3.0 with no flags