File: factnm.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (113 lines) | stat: -rw-r--r-- 3,460 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113

! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU Lesser General Public
! License. See the file COPYING for license details.

!BOP
! !ROUTINE: factnm
! !INTERFACE:
real(8) function factnm(n,m)
! !INPUT/OUTPUT PARAMETERS:
!   n : input (in,integer)
!   m : order of multifactorial (in,integer)
! !DESCRIPTION:
!   Returns the multifactorial
!   $$ n\underbrace{!!\,\cdots\,!}_{m\,{\rm times}}=
!    \prod_{\substack{i\ge 0\\ n-im>0}}(n-im) $$
!   for $n,\,m \ge 0$. $n$ should be less than 150.
!
! !REVISION HISTORY:
!   Created January 2003 (JKD)
!EOP
!BOC
implicit none
! arguments
integer, intent(in) :: n,m
! local variables
integer i,j
real(8) f1(24),f2(38)
data f1 / &
                       1.d0,                        2.d0,  &
                       6.d0,                       24.d0,  &
                     120.d0,                      720.d0,  &
                    5040.d0,                    40320.d0,  &
                  362880.d0,                  3628800.d0,  &
                39916800.d0,                479001600.d0,  &
              6227020800.d0,              87178291200.d0,  &
           1307674368000.d0,           20922789888000.d0,  &
         355687428096000.d0,         6402373705728000.d0,  &
      121645100408832000.d0,      2432902008176640000.d0,  &
    51090942171709440000.d0,   1124000727777607680000.d0,  &
 25852016738884976640000.d0, 620448401733239439360000.d0 /
data f2 / &
                       1.d0,                        2.d0,  &
                       3.d0,                        8.d0,  &
                      15.d0,                       48.d0,  &
                     105.d0,                      384.d0,  &
                     945.d0,                     3840.d0,  &
                   10395.d0,                    46080.d0,  &
                  135135.d0,                   645120.d0,  &
                 2027025.d0,                 10321920.d0,  &
                34459425.d0,                185794560.d0,  &
               654729075.d0,               3715891200.d0,  &
             13749310575.d0,              81749606400.d0,  &
            316234143225.d0,            1961990553600.d0,  &
           7905853580625.d0,           51011754393600.d0,  &
         213458046676875.d0,         1428329123020800.d0,  &
        6190283353629375.d0,        42849873690624000.d0,  &
      191898783962510625.d0,      1371195958099968000.d0,  &
     6332659870762850625.d0,     46620662575398912000.d0,  &
   221643095476699771875.d0,   1678343852714360832000.d0,  &
  8200794532637891559375.d0,  63777066403145711616000.d0 /
! fast return if possible
if (n.eq.0) then
  factnm=1.d0
  return
end if
if (m.eq.1) then
  if ((n.ge.1).and.(n.le.24)) then
    factnm=f1(n)
    return
  end if
end if
if (m.eq.2) then
  if ((n.ge.1).and.(n.le.38)) then
    factnm=f2(n)
    return
  end if
end if
if (n.lt.0) then
  write(*,*)
  write(*,'("Error(factnm): n < 0 : ",I8)') n
  write(*,*)
  stop
end if
if (m.le.0) then
  write(*,*)
  write(*,'("Error(factnm): m <= 0 : ",I8)') m
  write(*,*)
  stop
end if
if (n.gt.150) then
  write(*,*)
  write(*,'("Error(factnm): n out of range : ",I8)') n
  write(*,*)
  stop
end if
if (m.eq.1) then
  factnm=f1(24)
  do i=25,n
    factnm=factnm*dble(i)
  end do
else
  j=n/m
  if (mod(n,m).eq.0) j=j-1
  factnm=dble(n)
  do i=1,j
    factnm=factnm*dble(n-i*m)
  end do
end if
return
end function
!EOC