1 ! Copyright (C) 2008 Doug Coleman, Slava Pestov, Aaron Schaefer.
2 ! See http://factorcode.org/license.txt for BSD license.
3 USING: combinators.short-circuit kernel math math.constants math.functions
4 math.vectors sequences ;
9 ! http://www.rskey.org/gamma.htm "Lanczos Approximation"
10 ! n=6: error ~ 3 x 10^-11
12 CONSTANT: gamma-g6 5.15
16 2.50662827563479526904 225.525584619175212544 -268.295973841304927459
17 80.9030806934622512966 -5.00757863970517583837 0.0114684895434781459556
20 : gamma-z ( x n -- seq )
21 [ + recip ] with map 1.0 0 pick set-nth ;
23 : (gamma-lanczos6) ( x -- log[gamma[x+1]] )
25 [ 0.5 + dup gamma-g6 + [ log * ] keep - ]
26 [ 6 gamma-z gamma-p6 v. log ] bi + ;
28 : gamma-lanczos6 ( x -- gamma[x] )
29 #! gamma(x) = gamma(x+1) / x
30 [ (gamma-lanczos6) exp ] keep / ;
32 : gammaln-lanczos6 ( x -- gammaln[x] )
33 #! log(gamma(x)) = log(gamma(x+1)) - log(x)
34 [ (gamma-lanczos6) ] keep log - ;
36 : gamma-neg ( gamma[abs[x]] x -- gamma[x] )
37 dup pi * sin * * pi neg swap / ; inline
42 #! gamma(x) = integral 0..inf [ t^(x-1) exp(-t) ] dt
43 #! gamma(n+1) = n! for n > 0
44 dup { [ 0.0 <= ] [ 1.0 mod zero? ] } 1&& [
47 [ abs gamma-lanczos6 ] keep dup 0 > [ drop ] [ gamma-neg ] if
50 : gammaln ( x -- gamma[x] )
51 #! gammaln(x) is an alternative when gamma(x)'s range
56 [ abs gammaln-lanczos6 ] keep dup 0 > [ drop ] [ gamma-neg ] if
59 : nth-root ( n x -- y )
62 ! Forth Scientific Library Algorithm #1
64 ! Evaluates the Real Exponential Integral,
65 ! E1(x) = - Ei(-x) = int_x^\infty exp^{-u}/u du for x > 0
66 ! using a rational approximation
68 ! Collected Algorithms from ACM, Volume 1 Algorithms 1-220,
69 ! 1980; Association for Computing Machinery Inc., New York,
72 ! (c) Copyright 1994 Everett F. Carter. Permission is granted by the
73 ! author to use this software for any application provided the
74 ! copyright notice is preserved.
77 #! For real values of x only. Accurate to 7 decimals.
79 dup 0.00107857 * 0.00976004 -
114 ! James Stirling's approximation for N!:
115 ! http://www.csse.monash.edu.au/~lloyd/tildeAlgDS/Numerical/Stirling/
117 : stirling-fact ( n -- fact )
120 [ 12 * recip 1+ ] tri * * ;