1 ! Copyright (C) 2008 Doug Coleman, Slava Pestov.
2 ! See http://factorcode.org/license.txt for BSD license.
3 USING: kernel math math.constants math.functions math.intervals
4 math.vectors namespaces sequences combinators.short-circuit ;
9 ! http://www.rskey.org/gamma.htm "Lanczos Approximation"
10 ! n=6: error ~ 3 x 10^-11
12 : gamma-g6 5.15 ; inline
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 + dup [ log * ] dip - ]
26 [ 6 gamma-z gamma-p6 v. log ] bi + ;
28 : gamma-lanczos6 ( x -- gamma[x] )
29 #! gamma(x) = gamma(x+1) / x
30 dup (gamma-lanczos6) exp swap / ;
32 : gammaln-lanczos6 ( x -- gammaln[x] )
33 #! log(gamma(x)) = log(gamma(x+1)) - log(x)
34 dup (gamma-lanczos6) swap 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 dup abs gamma-lanczos6 swap dup 0 > [ drop ] [ gamma-neg ] if
50 : gammaln ( x -- gamma[x] )
51 #! gammaln(x) is an alternative when gamma(x)'s range
56 dup abs gammaln-lanczos6 swap dup 0 > [ drop ] [ gamma-neg ] if
59 : nth-root ( n x -- y )
60 [ recip ] dip swap ^ ;
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 )