2015-04-01 17:03:53 -04:00
|
|
|
! Copyright (C) 2011 John Benediktsson
|
|
|
|
|
! See http://factorcode.org/license.txt for BSD license
|
|
|
|
|
|
|
|
|
|
USING: combinators kernel locals math math.constants
|
|
|
|
|
math.functions sequences ;
|
|
|
|
|
|
|
|
|
|
IN: picomath
|
|
|
|
|
|
|
|
|
|
<PRIVATE
|
|
|
|
|
|
|
|
|
|
CONSTANT: a1 0.254829592
|
|
|
|
|
CONSTANT: a2 -0.284496736
|
|
|
|
|
CONSTANT: a3 1.421413741
|
|
|
|
|
CONSTANT: a4 -1.453152027
|
|
|
|
|
CONSTANT: a5 1.061405429
|
|
|
|
|
CONSTANT: p 0.3275911
|
|
|
|
|
|
|
|
|
|
PRIVATE>
|
|
|
|
|
|
|
|
|
|
! Standalone error function erf(x)
|
|
|
|
|
! http://www.johndcook.com/blog/2009/01/19/stand-alone-error-function-erf/
|
|
|
|
|
:: erf ( x -- value )
|
|
|
|
|
x 0 >= 1 -1 ? :> sign
|
|
|
|
|
x abs :> x!
|
2015-05-13 17:47:15 -04:00
|
|
|
p x * 1 + recip :> t
|
2015-04-01 17:03:53 -04:00
|
|
|
a5 t * a4 + t * a3 + t * a2 + t * a1 + t *
|
|
|
|
|
x x neg * e^ * 1 swap - :> y
|
|
|
|
|
sign y * ;
|
|
|
|
|
|
|
|
|
|
:: expm1 ( x -- value )
|
|
|
|
|
x abs 1e-5 < [ x sq 0.5 * x + ] [ x e^ 1.0 - ] if ;
|
|
|
|
|
|
|
|
|
|
! Standalone implementation of phi(x)
|
|
|
|
|
:: phi ( x -- value )
|
|
|
|
|
x 0 >= 1 -1 ? :> sign
|
|
|
|
|
x abs 2 sqrt / :> x!
|
2015-05-13 17:47:15 -04:00
|
|
|
p x * 1 + recip :> t
|
2015-04-01 17:03:53 -04:00
|
|
|
a5 t * a4 + t * a3 + t * a2 + t * a1 + t *
|
|
|
|
|
x x neg * e^ * 1 swap - :> y
|
|
|
|
|
sign y * 1 + 2 / ;
|
|
|
|
|
|
|
|
|
|
<PRIVATE
|
|
|
|
|
|
|
|
|
|
CONSTANT: lf {
|
|
|
|
|
0.000000000000000
|
|
|
|
|
0.000000000000000
|
|
|
|
|
0.693147180559945
|
|
|
|
|
1.791759469228055
|
|
|
|
|
3.178053830347946
|
|
|
|
|
4.787491742782046
|
|
|
|
|
6.579251212010101
|
|
|
|
|
8.525161361065415
|
|
|
|
|
10.604602902745251
|
|
|
|
|
12.801827480081469
|
|
|
|
|
15.104412573075516
|
|
|
|
|
17.502307845873887
|
|
|
|
|
19.987214495661885
|
|
|
|
|
22.552163853123421
|
|
|
|
|
25.191221182738683
|
|
|
|
|
27.899271383840894
|
|
|
|
|
30.671860106080675
|
|
|
|
|
33.505073450136891
|
|
|
|
|
36.395445208033053
|
|
|
|
|
39.339884187199495
|
|
|
|
|
42.335616460753485
|
|
|
|
|
45.380138898476908
|
|
|
|
|
48.471181351835227
|
|
|
|
|
51.606675567764377
|
|
|
|
|
54.784729398112319
|
|
|
|
|
58.003605222980518
|
|
|
|
|
61.261701761002001
|
|
|
|
|
64.557538627006323
|
|
|
|
|
67.889743137181526
|
|
|
|
|
71.257038967168000
|
|
|
|
|
74.658236348830158
|
|
|
|
|
78.092223553315307
|
|
|
|
|
81.557959456115029
|
|
|
|
|
85.054467017581516
|
|
|
|
|
88.580827542197682
|
|
|
|
|
92.136175603687079
|
|
|
|
|
95.719694542143202
|
|
|
|
|
99.330612454787428
|
|
|
|
|
102.968198614513810
|
|
|
|
|
106.631760260643450
|
|
|
|
|
110.320639714757390
|
|
|
|
|
114.034211781461690
|
|
|
|
|
117.771881399745060
|
|
|
|
|
121.533081515438640
|
|
|
|
|
125.317271149356880
|
|
|
|
|
129.123933639127240
|
|
|
|
|
132.952575035616290
|
|
|
|
|
136.802722637326350
|
|
|
|
|
140.673923648234250
|
|
|
|
|
144.565743946344900
|
|
|
|
|
148.477766951773020
|
|
|
|
|
152.409592584497350
|
|
|
|
|
156.360836303078800
|
|
|
|
|
160.331128216630930
|
|
|
|
|
164.320112263195170
|
|
|
|
|
168.327445448427650
|
|
|
|
|
172.352797139162820
|
|
|
|
|
176.395848406997370
|
|
|
|
|
180.456291417543780
|
|
|
|
|
184.533828861449510
|
|
|
|
|
188.628173423671600
|
|
|
|
|
192.739047287844900
|
|
|
|
|
196.866181672889980
|
|
|
|
|
201.009316399281570
|
|
|
|
|
205.168199482641200
|
|
|
|
|
209.342586752536820
|
|
|
|
|
213.532241494563270
|
|
|
|
|
217.736934113954250
|
|
|
|
|
221.956441819130360
|
|
|
|
|
226.190548323727570
|
|
|
|
|
230.439043565776930
|
|
|
|
|
234.701723442818260
|
|
|
|
|
238.978389561834350
|
|
|
|
|
243.268849002982730
|
|
|
|
|
247.572914096186910
|
|
|
|
|
251.890402209723190
|
|
|
|
|
256.221135550009480
|
|
|
|
|
260.564940971863220
|
|
|
|
|
264.921649798552780
|
|
|
|
|
269.291097651019810
|
|
|
|
|
273.673124285693690
|
|
|
|
|
278.067573440366120
|
|
|
|
|
282.474292687630400
|
|
|
|
|
286.893133295426990
|
|
|
|
|
291.323950094270290
|
|
|
|
|
295.766601350760600
|
|
|
|
|
300.220948647014100
|
|
|
|
|
304.686856765668720
|
|
|
|
|
309.164193580146900
|
|
|
|
|
313.652829949878990
|
|
|
|
|
318.152639620209300
|
|
|
|
|
322.663499126726210
|
|
|
|
|
327.185287703775200
|
|
|
|
|
331.717887196928470
|
|
|
|
|
336.261181979198450
|
|
|
|
|
340.815058870798960
|
|
|
|
|
345.379407062266860
|
|
|
|
|
349.954118040770250
|
|
|
|
|
354.539085519440790
|
|
|
|
|
359.134205369575340
|
|
|
|
|
363.739375555563470
|
|
|
|
|
368.354496072404690
|
|
|
|
|
372.979468885689020
|
|
|
|
|
377.614197873918670
|
|
|
|
|
382.258588773060010
|
|
|
|
|
386.912549123217560
|
|
|
|
|
391.575988217329610
|
|
|
|
|
396.248817051791490
|
|
|
|
|
400.930948278915760
|
|
|
|
|
405.622296161144900
|
|
|
|
|
410.322776526937280
|
|
|
|
|
415.032306728249580
|
|
|
|
|
419.750805599544780
|
|
|
|
|
424.478193418257090
|
|
|
|
|
429.214391866651570
|
|
|
|
|
433.959323995014870
|
|
|
|
|
438.712914186121170
|
|
|
|
|
443.475088120918940
|
|
|
|
|
448.245772745384610
|
|
|
|
|
453.024896238496130
|
|
|
|
|
457.812387981278110
|
|
|
|
|
462.608178526874890
|
|
|
|
|
467.412199571608080
|
|
|
|
|
472.224383926980520
|
|
|
|
|
477.044665492585580
|
|
|
|
|
481.872979229887900
|
|
|
|
|
486.709261136839360
|
|
|
|
|
491.553448223298010
|
|
|
|
|
496.405478487217580
|
|
|
|
|
501.265290891579240
|
|
|
|
|
506.132825342034830
|
|
|
|
|
511.008022665236070
|
|
|
|
|
515.890824587822520
|
|
|
|
|
520.781173716044240
|
|
|
|
|
525.679013515995050
|
|
|
|
|
530.584288294433580
|
|
|
|
|
535.496943180169520
|
|
|
|
|
540.416924105997740
|
|
|
|
|
545.344177791154950
|
|
|
|
|
550.278651724285620
|
|
|
|
|
555.220294146894960
|
|
|
|
|
560.169054037273100
|
|
|
|
|
565.124881094874350
|
|
|
|
|
570.087725725134190
|
|
|
|
|
575.057539024710200
|
|
|
|
|
580.034272767130800
|
|
|
|
|
585.017879388839220
|
|
|
|
|
590.008311975617860
|
|
|
|
|
595.005524249382010
|
|
|
|
|
600.009470555327430
|
|
|
|
|
605.020105849423770
|
|
|
|
|
610.037385686238740
|
|
|
|
|
615.061266207084940
|
|
|
|
|
620.091704128477430
|
|
|
|
|
625.128656730891070
|
|
|
|
|
630.172081847810200
|
|
|
|
|
635.221937855059760
|
|
|
|
|
640.278183660408100
|
|
|
|
|
645.340778693435030
|
|
|
|
|
650.409682895655240
|
|
|
|
|
655.484856710889060
|
|
|
|
|
660.566261075873510
|
|
|
|
|
665.653857411105950
|
|
|
|
|
670.747607611912710
|
|
|
|
|
675.847474039736880
|
|
|
|
|
680.953419513637530
|
|
|
|
|
686.065407301994010
|
|
|
|
|
691.183401114410800
|
|
|
|
|
696.307365093814040
|
|
|
|
|
701.437263808737160
|
|
|
|
|
706.573062245787470
|
|
|
|
|
711.714725802289990
|
|
|
|
|
716.862220279103440
|
|
|
|
|
722.015511873601330
|
|
|
|
|
727.174567172815840
|
|
|
|
|
732.339353146739310
|
|
|
|
|
737.509837141777440
|
|
|
|
|
742.685986874351220
|
|
|
|
|
747.867770424643370
|
|
|
|
|
753.055156230484160
|
|
|
|
|
758.248113081374300
|
|
|
|
|
763.446610112640200
|
|
|
|
|
768.650616799717000
|
|
|
|
|
773.860102952558460
|
|
|
|
|
779.075038710167410
|
|
|
|
|
784.295394535245690
|
|
|
|
|
789.521141208958970
|
|
|
|
|
794.752249825813460
|
|
|
|
|
799.988691788643450
|
|
|
|
|
805.230438803703120
|
|
|
|
|
810.477462875863580
|
|
|
|
|
815.729736303910160
|
|
|
|
|
820.987231675937890
|
|
|
|
|
826.249921864842800
|
|
|
|
|
831.517780023906310
|
|
|
|
|
836.790779582469900
|
|
|
|
|
842.068894241700490
|
|
|
|
|
847.352097970438420
|
|
|
|
|
852.640365001133090
|
|
|
|
|
857.933669825857460
|
|
|
|
|
863.231987192405430
|
|
|
|
|
868.535292100464630
|
|
|
|
|
873.843559797865740
|
|
|
|
|
879.156765776907600
|
|
|
|
|
884.474885770751830
|
|
|
|
|
889.797895749890240
|
|
|
|
|
895.125771918679900
|
|
|
|
|
900.458490711945270
|
|
|
|
|
905.796028791646340
|
|
|
|
|
911.138363043611210
|
|
|
|
|
916.485470574328820
|
|
|
|
|
921.837328707804890
|
|
|
|
|
927.193914982476710
|
|
|
|
|
932.555207148186240
|
|
|
|
|
937.921183163208070
|
|
|
|
|
943.291821191335660
|
|
|
|
|
948.667099599019820
|
|
|
|
|
954.046996952560450
|
|
|
|
|
959.431492015349480
|
|
|
|
|
964.820563745165940
|
|
|
|
|
970.214191291518320
|
|
|
|
|
975.612353993036210
|
|
|
|
|
981.015031374908400
|
|
|
|
|
986.422203146368590
|
|
|
|
|
991.833849198223450
|
|
|
|
|
997.249949600427840
|
|
|
|
|
1002.670484599700300
|
|
|
|
|
1008.095434617181700
|
|
|
|
|
1013.524780246136200
|
|
|
|
|
1018.958502249690200
|
|
|
|
|
1024.396581558613400
|
|
|
|
|
1029.838999269135500
|
|
|
|
|
1035.285736640801600
|
|
|
|
|
1040.736775094367400
|
|
|
|
|
1046.192096209724900
|
|
|
|
|
1051.651681723869200
|
|
|
|
|
1057.115513528895000
|
|
|
|
|
1062.583573670030100
|
|
|
|
|
1068.055844343701400
|
|
|
|
|
1073.532307895632800
|
|
|
|
|
1079.012946818975000
|
|
|
|
|
1084.497743752465600
|
|
|
|
|
1089.986681478622400
|
|
|
|
|
1095.479742921962700
|
|
|
|
|
1100.976911147256000
|
|
|
|
|
1106.478169357800900
|
|
|
|
|
1111.983500893733000
|
|
|
|
|
1117.492889230361000
|
|
|
|
|
1123.006317976526100
|
|
|
|
|
1128.523770872990800
|
|
|
|
|
1134.045231790853000
|
|
|
|
|
1139.570684729984800
|
|
|
|
|
1145.100113817496100
|
|
|
|
|
1150.633503306223700
|
|
|
|
|
1156.170837573242400
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
PRIVATE>
|
|
|
|
|
|
|
|
|
|
: log-factorial ( n -- value )
|
|
|
|
|
{
|
|
|
|
|
{ [ dup 0 < ] [ "invalid input" throw ] }
|
|
|
|
|
{ [ dup 254 > ] [
|
|
|
|
|
1 + [| x |
|
|
|
|
|
x 0.5 - x log * x -
|
|
|
|
|
0.5 2 pi * log * +
|
|
|
|
|
1.0 12.0 x * / +
|
|
|
|
|
] call
|
|
|
|
|
] }
|
|
|
|
|
[ lf nth ]
|
|
|
|
|
} cond ;
|
|
|
|
|
|
|
|
|
|
:: log-one-plus-x ( x -- value )
|
|
|
|
|
x -1.0 <= [ "argument must be > -1" throw ] when
|
|
|
|
|
x abs 1e-4 > [ 1.0 x + log ] [ -0.5 x * 1.0 + x * ] if ;
|
|
|
|
|
|
|
|
|
|
<PRIVATE
|
|
|
|
|
|
|
|
|
|
CONSTANT: c0 2.515517
|
|
|
|
|
CONSTANT: c1 0.802853
|
|
|
|
|
CONSTANT: c2 0.010328
|
|
|
|
|
|
|
|
|
|
CONSTANT: d0 1.432788
|
|
|
|
|
CONSTANT: d1 0.189269
|
|
|
|
|
CONSTANT: d2 0.001308
|
|
|
|
|
|
|
|
|
|
PRIVATE>
|
|
|
|
|
|
|
|
|
|
:: rational-approximation ( t -- value )
|
|
|
|
|
c2 t * c1 + t * c0 + :> numerator
|
|
|
|
|
d2 t * d1 + t * d0 + t * 1.0 + :> denominator
|
|
|
|
|
t numerator denominator / - ;
|
|
|
|
|
|
|
|
|
|
:: normal-cdf-inverse ( p -- value )
|
|
|
|
|
p [ 0 > ] [ 1 < ] bi and [ p throw ] unless
|
|
|
|
|
p 0.5 <
|
|
|
|
|
[ p log -2.0 * sqrt rational-approximation neg ]
|
|
|
|
|
[ p 1.0 - log -2.0 * sqrt rational-approximation ] if ;
|
|
|
|
|
|
|
|
|
|
<PRIVATE
|
|
|
|
|
|
|
|
|
|
! Abramowitz and Stegun 6.1.41
|
|
|
|
|
! Asymptotic series should be good to at least 11 or 12 figures
|
|
|
|
|
! For error analysis, see Whittiker and Watson
|
|
|
|
|
! A Course in Modern Analysis (1927), page 252
|
|
|
|
|
CONSTANT: c {
|
|
|
|
|
1/12
|
|
|
|
|
-1/360
|
|
|
|
|
1/1260
|
|
|
|
|
-1/1680
|
|
|
|
|
1/1188
|
|
|
|
|
-691/360360
|
|
|
|
|
1/156
|
|
|
|
|
-3617/122400
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
CONSTANT: halfLogTwoPi 0.91893853320467274178032973640562
|
|
|
|
|
|
|
|
|
|
PRIVATE>
|
|
|
|
|
|
|
|
|
|
DEFER: gamma
|
|
|
|
|
|
|
|
|
|
:: log-gamma ( x -- value )
|
|
|
|
|
x 0 <= [ "Invalid input" throw ] when
|
|
|
|
|
x 12 < [ x gamma abs log ] [
|
2015-06-29 19:43:15 -04:00
|
|
|
1.0 x x * / :> z
|
2017-06-01 17:59:35 -04:00
|
|
|
7 c nth 7 <iota> reverse [ [ z * ] [ c nth ] bi* + ] each x / :> series
|
2015-04-01 17:03:53 -04:00
|
|
|
x 0.5 - x log * x - halfLogTwoPi + series +
|
|
|
|
|
] if ;
|
|
|
|
|
|
|
|
|
|
<PRIVATE
|
|
|
|
|
|
|
|
|
|
CONSTANT: GAMMA 0.577215664901532860606512090 ! Euler's gamma constant
|
|
|
|
|
|
|
|
|
|
! numerator coefficients for approximation over the interval (1,2)
|
|
|
|
|
CONSTANT: P {
|
|
|
|
|
-1.71618513886549492533811E+0
|
|
|
|
|
2.47656508055759199108314E+1
|
|
|
|
|
-3.79804256470945635097577E+2
|
|
|
|
|
6.29331155312818442661052E+2
|
|
|
|
|
8.66966202790413211295064E+2
|
|
|
|
|
-3.14512729688483675254357E+4
|
|
|
|
|
-3.61444134186911729807069E+4
|
|
|
|
|
6.64561438202405440627855E+4
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
! denominator coefficients for approximation over the interval (1,2)
|
|
|
|
|
CONSTANT: Q {
|
|
|
|
|
-3.08402300119738975254353E+1
|
|
|
|
|
3.15350626979604161529144E+2
|
|
|
|
|
-1.01515636749021914166146E+3
|
|
|
|
|
-3.10777167157231109440444E+3
|
|
|
|
|
2.25381184209801510330112E+4
|
|
|
|
|
4.75584627752788110767815E+3
|
|
|
|
|
-1.34659959864969306392456E+5
|
|
|
|
|
-1.15132259675553483497211E+5
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
:: (gamma) ( x -- value )
|
|
|
|
|
! The algorithm directly approximates gamma over (1,2) and uses
|
|
|
|
|
! reduction identities to reduce other arguments to this interval.
|
|
|
|
|
x :> y!
|
|
|
|
|
0 :> n!
|
|
|
|
|
y 1.0 < :> arg-was-less-than-one
|
|
|
|
|
arg-was-less-than-one
|
|
|
|
|
[ y 1.0 + y! ] [ y floor >integer 1 - n! y n - y! ] if
|
|
|
|
|
0.0 :> num!
|
|
|
|
|
1.0 :> den!
|
|
|
|
|
y 1 - :> z!
|
2017-06-01 17:59:35 -04:00
|
|
|
8 <iota> [
|
2015-04-01 17:03:53 -04:00
|
|
|
[ P nth num + z * num! ]
|
|
|
|
|
[ Q nth den z * + den! ] bi
|
|
|
|
|
] each
|
|
|
|
|
num den / 1.0 +
|
|
|
|
|
arg-was-less-than-one
|
|
|
|
|
[ y 1.0 - / ] [ n [ y * y 1.0 + y! ] times ] if ;
|
|
|
|
|
|
|
|
|
|
PRIVATE>
|
|
|
|
|
|
|
|
|
|
:: gamma ( x -- value )
|
|
|
|
|
x 0 <= [ "Invalid input" throw ] when
|
|
|
|
|
x {
|
|
|
|
|
{ [ dup 0.001 < ] [ GAMMA * 1.0 + x * 1.0 swap / ] }
|
|
|
|
|
{ [ dup 12.0 < ] [ (gamma) ] }
|
|
|
|
|
{ [ dup 171.624 > ] [ drop 1/0. ] }
|
|
|
|
|
[ log-gamma e^ ]
|
|
|
|
|
} cond ;
|
