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