fortran-lang/stdlib

Include factorials in `stdlib_constants` module

Opened this issue · 0 comments

Since factorials often arise in calculations, I suggest adding to stdlib_constants a named array constant factorial(1:20) of type int64 containing the factorials of 1 through 20, since that is what int64 can hold. Similar named constants such as double_factorial could also be added. An alternative to defining a single named constant factorial would be to define named array constants factorial_int8, factorial_int16, factorial_int32, factorial_int64. factorial_real32, factorial_real64, factorial_real128 holding the factorials that fit into the corresponding data types. The advantages of making factorial a named constant instead of a function is speed and avoiding overflow. Named constants for the absolute moments of the normal distribution (which involve double factorials) are useful to me, but maybe this is too specialized for stdlib.

The code

program xfactorial
use iso_fortran_env, only: int64
implicit none
integer, parameter :: dp = kind(1.0d0), ikind=int64
integer(kind=int64) :: i, fac
real(kind=dp) :: x
fac = 1_ikind
x = 1.0_dp
print*,huge(i)
do i=2,21
   fac = fac*i
   x = x*i
   print*,i,fac,x,fac/x
end do
end program xfactorial

gives

  9223372036854775807
                    2                    2   2.0000000000000000        1.0000000000000000     
                    3                    6   6.0000000000000000        1.0000000000000000     
                    4                   24   24.000000000000000        1.0000000000000000     
                    5                  120   120.00000000000000        1.0000000000000000     
                    6                  720   720.00000000000000        1.0000000000000000     
                    7                 5040   5040.0000000000000        1.0000000000000000     
                    8                40320   40320.000000000000        1.0000000000000000     
                    9               362880   362880.00000000000        1.0000000000000000     
                   10              3628800   3628800.0000000000        1.0000000000000000     
                   11             39916800   39916800.000000000        1.0000000000000000     
                   12            479001600   479001600.00000000        1.0000000000000000     
                   13           6227020800   6227020800.0000000        1.0000000000000000     
                   14          87178291200   87178291200.000000        1.0000000000000000     
                   15        1307674368000   1307674368000.0000        1.0000000000000000     
                   16       20922789888000   20922789888000.000        1.0000000000000000     
                   17      355687428096000   355687428096000.00        1.0000000000000000     
                   18     6402373705728000   6402373705728000.0        1.0000000000000000     
                   19   121645100408832000   1.2164510040883200E+017   1.0000000000000000     
                   20  2432902008176640000   2.4329020081766400E+018   1.0000000000000000     
                   21 -4249290049419214848   5.1090942171709440E+019  -8.3171103698537238E-002