--- /dev/null
+! Copyright (C) 2011 John Benediktsson
+! See http://factorcode.org/license.txt for BSD license
+
+USING: kernel literals math math.functions math.ranges picomath
+sequences tools.test ;
+
+IN: picomath
+
+[ t ] [
+ {
+ { -3 -0.999977909503 }
+ { -1 -0.842700792950 }
+ { 0.0 0.0 }
+ { 0.5 0.520499877813 }
+ { 2.1 0.997020533344 }
+ } [ [ first erf ] [ second - ] bi abs ] map
+ supremum 1e-6 <
+] unit-test
+
+[ t ] [
+ {
+ { -1 -0.632120558828558 }
+ { 0.0 0.0 }
+ { $[ 1e-5 1e-8 - ] 0.000009990049900216168 }
+ { $[ 1e-5 1e-8 + ] 0.00001001005010021717 }
+ { 0.5 0.6487212707001282 }
+ } [ [ first expm1 ] [ second - ] bi abs ] map
+ supremum 1e-6 <
+] unit-test
+
+[ t ] [
+ {
+ { -3 0.00134989803163 }
+ { -1 0.158655253931 }
+ { 0.0 0.5 }
+ { 0.5 0.691462461274 }
+ { 2.1 0.982135579437 }
+ } [ [ first phi ] [ second - ] bi abs ] map
+ supremum 1e-3 <
+] unit-test
+
+: factorial ( n -- n! ) [ 1 ] [ [1,b] 1 [ * ] reduce ] if-zero ;
+
+[ t ] [
+ { 0 1 10 100 1000 10000 } [
+ [ factorial log ] [ log-factorial ] bi - abs
+ ] map supremum 1e-6 <
+] unit-test
+
+: relative-error ( approx value -- relative-error )
+ [ - abs ] keep / ;
+
+[ t ] [
+ {
+ { 1e-20 1e+20 }
+ { 2.19824158876e-16 4.5490905327e+15 } ! 0.99*DBL_EPSILON
+ { 2.24265050974e-16 4.45900953205e+15 } ! 1.01*DBL_EPSILON
+ { 0.00099 1009.52477271 }
+ { 0.00100 999.423772485 }
+ { 0.00101 989.522792258 }
+ { 6.1 142.451944066 }
+ { 11.999 39819417.4793 }
+ { 12 39916800.0 }
+ { 12.001 40014424.1571 }
+ { 15.2 149037380723.0 }
+ } [ [ first gamma ] [ second relative-error ] bi ] map
+ supremum 1e-6 <
+] unit-test
+
+[ t ] [
+ {
+ { 1e-12 27.6310211159 }
+ { 0.9999 5.77297915613e-05 }
+ { 1.0001 -5.77133422205e-05 }
+ { 3.1 0.787375083274 }
+ { 6.3 5.30734288962 }
+ { 11.9999 17.5020635801 }
+ { 12 17.5023078459 }
+ { 12.0001 17.5025521125 }
+ { 27.4 62.5755868211 }
+ } [ [ first log-gamma ] [ second relative-error ] bi ] map
+ supremum 1e-10 <
+] unit-test
--- /dev/null
+! 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!
+ p x * 1 + 1 swap / :> t
+ 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!
+ p x * 1 + 1 swap / :> t
+ 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 ] [
+ 1.0 x x * / :> z
+ 7 c nth 7 iota reverse [ [ z * ] [ c nth ] bi* + ] each x / :> series
+ 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!
+ 8 iota [
+ [ 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 ;
projects / factor.git / commitdiff