fac.f
1.64 KB
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
function fac (n)
c june 1977 edition. w. fullerton, c3, los alamos scientific lab.
dimension facn(26)
external r9lgmc
c
data facn( 1) / 1.0e0 /
data facn( 2) / 1.0e0 /
data facn( 3) / 2.0e0 /
data facn( 4) / 6.0e0 /
data facn( 5) / 24.0e0 /
data facn( 6) / 120.0e0 /
data facn( 7) / 720.0e0 /
data facn( 8) / 5040.0e0 /
data facn( 9) / 40320.0e0 /
data facn(10) / 362880.0e0 /
data facn(11) / 3628800.0e0 /
data facn(12) / 39916800.0e0 /
data facn(13) / 479001600.0e0 /
data facn(14) / 6227020800.0e0 /
data facn(15) / 87178291200.0e0 /
data facn(16) / 1307674368000.0e0 /
data facn(17) / 20922789888000.0e0 /
data facn(18) / 355687428096000.0e0 /
data facn(19) / 6402373705728000.0e0 /
data facn(20) / .12164510040883200e18 /
data facn(21) / .24329020081766400e19 /
data facn(22) / .51090942171709440e20 /
data facn(23) / .11240007277776077e22 /
data facn(24) / .25852016738884977e23 /
data facn(25) / .62044840173323944e24 /
data facn(26) / .15511210043330986e26 /
c
data sq2pil / 0.9189385332 0467274e0/
data nmax / 0 /
c
if (nmax.ne.0) go to 10
call r9gaml (xmin, xmax)
nmax = xmax - 1.
c
10 if (n.lt.0) call seteru (
1 47hfac factorial of negative integer undefined, 47, 1, 2)
c
if (n.le.25) fac = facn(n+1)
if (n.le.25) return
c
if (n.gt.nmax) call seteru (
1 39hfac n so big factorial(n) overflows, 39, 2, 2)
c
x = n + 1
fac = exp ( (x-0.5)*alog(x) - x + sq2pil + r9lgmc(x) )
c
return
end