1 ! Copyright (C) 2011 John Benediktsson
2 ! See http://factorcode.org/license.txt for BSD license
4 USING: combinators kernel locals math math.constants
5 math.functions sequences ;
11 CONSTANT: a1 0.254829592
12 CONSTANT: a2 -0.284496736
13 CONSTANT: a3 1.421413741
14 CONSTANT: a4 -1.453152027
15 CONSTANT: a5 1.061405429
20 ! Standalone error function erf(x)
21 ! http://www.johndcook.com/blog/2009/01/19/stand-alone-error-function-erf/
26 a5 t * a4 + t * a3 + t * a2 + t * a1 + t *
27 x x neg * e^ * 1 swap - :> y
30 :: expm1 ( x -- value )
31 x abs 1e-5 < [ x sq 0.5 * x + ] [ x e^ 1.0 - ] if ;
33 ! Standalone implementation of phi(x)
38 a5 t * a4 + t * a3 + t * a2 + t * a1 + t *
39 x x neg * e^ * 1 swap - :> y
304 : log-factorial ( n -- value )
306 { [ dup 0 < ] [ "invalid input" throw ] }
317 :: log-one-plus-x ( x -- value )
318 x -1.0 <= [ "argument must be > -1" throw ] when
319 x abs 1e-4 > [ 1.0 x + log ] [ -0.5 x * 1.0 + x * ] if ;
323 CONSTANT: c0 2.515517
324 CONSTANT: c1 0.802853
325 CONSTANT: c2 0.010328
327 CONSTANT: d0 1.432788
328 CONSTANT: d1 0.189269
329 CONSTANT: d2 0.001308
333 :: rational-approximation ( t -- value )
334 c2 t * c1 + t * c0 + :> numerator
335 d2 t * d1 + t * d0 + t * 1.0 + :> denominator
336 t numerator denominator / - ;
338 :: normal-cdf-inverse ( p -- value )
339 p [ 0 > ] [ 1 < ] bi and [ p throw ] unless
341 [ p log -2.0 * sqrt rational-approximation neg ]
342 [ p 1.0 - log -2.0 * sqrt rational-approximation ] if ;
346 ! Abramowitz and Stegun 6.1.41
347 ! Asymptotic series should be good to at least 11 or 12 figures
348 ! For error analysis, see Whittiker and Watson
349 ! A Course in Modern Analysis (1927), page 252
361 CONSTANT: halfLogTwoPi 0.91893853320467274178032973640562
367 :: log-gamma ( x -- value )
368 x 0 <= [ "Invalid input" throw ] when
369 x 12 < [ x gamma abs log ] [
371 7 c nth 7 <iota> reverse [ [ z * ] [ c nth ] bi* + ] each x / :> series
372 x 0.5 - x log * x - halfLogTwoPi + series +
377 CONSTANT: GAMMA 0.577215664901532860606512090 ! Euler's gamma constant
379 ! numerator coefficients for approximation over the interval (1,2)
381 -1.71618513886549492533811E+0
382 2.47656508055759199108314E+1
383 -3.79804256470945635097577E+2
384 6.29331155312818442661052E+2
385 8.66966202790413211295064E+2
386 -3.14512729688483675254357E+4
387 -3.61444134186911729807069E+4
388 6.64561438202405440627855E+4
391 ! denominator coefficients for approximation over the interval (1,2)
393 -3.08402300119738975254353E+1
394 3.15350626979604161529144E+2
395 -1.01515636749021914166146E+3
396 -3.10777167157231109440444E+3
397 2.25381184209801510330112E+4
398 4.75584627752788110767815E+3
399 -1.34659959864969306392456E+5
400 -1.15132259675553483497211E+5
403 :: (gamma) ( x -- value )
404 ! The algorithm directly approximates gamma over (1,2) and uses
405 ! reduction identities to reduce other arguments to this interval.
408 y 1.0 < :> arg-was-less-than-one
409 arg-was-less-than-one
410 [ y 1.0 + y! ] [ y floor >integer 1 - n! y n - y! ] if
415 [ P nth num + z * num! ]
416 [ Q nth den z * + den! ] bi
419 arg-was-less-than-one
420 [ y 1.0 - / ] [ n [ y * y 1.0 + y! ] times ] if ;
424 :: gamma ( x -- value )
425 x 0 <= [ "Invalid input" throw ] when
427 { [ dup 0.001 < ] [ GAMMA * 1.0 + x * 1.0 swap / ] }
428 { [ dup 12.0 < ] [ (gamma) ] }
429 { [ dup 171.624 > ] [ drop 1/0. ] }
projects / factor.git / blob