-
Notifications
You must be signed in to change notification settings - Fork 0
/
combi.f
93 lines (74 loc) · 2.57 KB
/
combi.f
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
module combi
implicit none
contains
c---------------------------------------------------------------------
integer function facult(n)
implicit none
integer, intent(in) :: n
integer i
facult=1
do i=2,n
facult=facult*i
end do
end function facult
c---------------------------------------------------------------------
function facult_large(n)
implicit none
integer, intent(in) :: n
integer(selected_int_kind(15)) :: facult_large
integer i
facult_large=1
do i=2,n
facult_large=facult_large*i
end do
end function facult_large
c---------------------------------------------------------------------
integer function dblfactorial(n)
implicit none
integer, intent(in) :: n
integer i
dblfactorial=1
do i=3,n,2
dblfactorial=dblfactorial*i
end do
end function
c---------------------------------------------------------------------
integer function n_over_k(n,k)
use defs, only : error_stop
implicit none
integer, intent(in) :: n,k
!
if (n.lt.0.or.k.lt.0) &
& call error_stop("Non-negative integer required")
!if (n.lt.k) call error_stop("n must be >= k")
if (n.lt.k) n_over_k=0
if (n.eq.k) n_over_k=1
if (n.gt.k) then
n_over_k=facult(n)/(facult(k)*facult(n-k))
end if
!
end function n_over_k
c---------------------------------------------------------------------
integer function n_over_k_large(n,k)
use defs, only : error_stop
implicit none
integer, intent(in) :: n,k
double precision facult_n,facult_k,facult_n_k
integer(selected_int_kind(15)) n_large,k_large
!
if (n.lt.0.or.k.lt.0) &
& call error_stop("Non-negative integer required")
!if (n.lt.k) call error_stop("n must be >= k")
if (n.lt.k) n_over_k_large=0
if (n==k) n_over_k_large=1
if (n.gt.k) then
!n_over_k=facult(n)/(facult(k)*facult(n-k))
facult_n=facult_large(n)
facult_k=facult_large(k)
facult_n_k=facult_large(n-k)
n_over_k_large=int(facult_n/(facult_k*facult_n_k))
end if
!
end function n_over_k_large
c---------------------------------------------------------------------
end module combi