Blame view
fvn_fnlib/fac.f
1.64 KB
38581db0c
git-svn-id: https...
|
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
|