26 lines
1.3 KiB
Forth
26 lines
1.3 KiB
Forth
2 constant (gamma-shift) \ don't change this
|
|
\ an approximation of the d(x) function
|
|
: ~d(x) ( f1 -- f2)
|
|
fdup 10 s>f f< \ use first symmetrical sigmoidal
|
|
if \ for range 1-10
|
|
-2705443e-8 fswap 2280802e-6 f/ 1428045e-6 f** 1 s>f f+ f/ 3187831e-8 f+
|
|
else \ use second symmetrical sigmoidal
|
|
-29372563e-9 fswap 1841693e-6 f/ 1052779e-6 f** 1 s>f f+ f/ 3330828e-8 f+
|
|
then 333333333e-10 fover f< if fdrop 1 s>f 30 s>f f/ then
|
|
; \ perform some sane clipping to infinity
|
|
|
|
: (ramanujan) ( f1 -- f2)
|
|
fdup fdup f* 4 s>f f* ( n 4n2)
|
|
fover fover f* fdup f+ f+ fover f+ ( n 8n3+4n2+n)
|
|
fover ~d(x) f+ ( n 8n3+4n2+n+d[x])
|
|
1 s>f 6 s>f f/ f** ( n 8n3+4n2+n+d[x]^1/6)
|
|
fswap fdup 2.7182818284590452353e f/ ( 8n3+4n2+n+d[x]^1/6 n n/e)
|
|
fswap f** f* pi fsqrt f* ( f)
|
|
;
|
|
|
|
: gamma ( f1 -- f2)
|
|
fdup f0< >r fdup f0= r> or abort" Gamma less or equal to zero"
|
|
fdup (gamma-shift) 1- s>f f+ (ramanujan) fswap
|
|
1 s>f (gamma-shift) 0 do fover i s>f f+ f* loop fswap fdrop f/
|
|
;
|