1 USING: kernel math math.constants math.functions math.intervals
2 math.vectors namespaces sequences ;
7 ! http://www.rskey.org/gamma.htm "Lanczos Approximation"
8 ! n=6: error ~ 3 x 10^-11
10 : gamma-g6 5.15 ; inline
14 2.50662827563479526904 225.525584619175212544 -268.295973841304927459
15 80.9030806934622512966 -5.00757863970517583837 0.0114684895434781459556
18 : gamma-z ( x n -- seq )
19 [ + recip ] with map 1.0 0 pick set-nth ;
21 : (gamma-lanczos6) ( x -- log[gamma[x+1]] )
23 dup 0.5 + dup gamma-g6 + dup >r log * r> -
24 swap 6 gamma-z gamma-p6 v. log + ;
26 : gamma-lanczos6 ( x -- gamma[x] )
27 #! gamma(x) = gamma(x+1) / x
28 dup (gamma-lanczos6) exp swap / ;
30 : gammaln-lanczos6 ( x -- gammaln[x] )
31 #! log(gamma(x)) = log(gamma(x+1)) - log(x)
32 dup (gamma-lanczos6) swap log - ;
34 : gamma-neg ( gamma[abs[x]] x -- gamma[x] )
35 dup pi * sin * * pi neg swap / ; inline
40 #! gamma(x) = integral 0..inf [ t^(x-1) exp(-t) ] dt
41 #! gamma(n+1) = n! for n > 0
42 dup 0.0 <= over 1.0 mod zero? and [
45 dup abs gamma-lanczos6 swap dup 0 > [ drop ] [ gamma-neg ] if
48 : gammaln ( x -- gamma[x] )
49 #! gammaln(x) is an alternative when gamma(x)'s range
54 dup abs gammaln-lanczos6 swap dup 0 > [ drop ] [ gamma-neg ] if
57 : nth-root ( n x -- y )
58 over 0 = [ "0th root is undefined" throw ] when >r recip r> swap ^ ;
60 ! Forth Scientific Library Algorithm #1
62 ! Evaluates the Real Exponential Integral,
63 ! E1(x) = - Ei(-x) = int_x^\infty exp^{-u}/u du for x > 0
64 ! using a rational approximation
66 ! Collected Algorithms from ACM, Volume 1 Algorithms 1-220,
67 ! 1980; Association for Computing Machinery Inc., New York,
70 ! (c) Copyright 1994 Everett F. Carter. Permission is granted by the
71 ! author to use this software for any application provided the
72 ! copyright notice is preserved.
75 #! For real values of x only. Accurate to 7 decimals.
77 dup 0.00107857 * 0.00976004 -