# ghewgill/picomath

1 parent 489164a commit 8162c0dc4dba345e11f2781ad2b7113bb2e42197 committed Jul 16, 2011
Showing with 481 additions and 0 deletions.
1. +18 −0 scheme/erf.scm
2. +4 −0 scheme/expm1.scm
3. +121 −0 scheme/gamma.scm
4. +263 −0 scheme/log-factorial.scm
5. +16 −0 scheme/normal-cdf-inverse.scm
6. +18 −0 scheme/phi.scm
7. +40 −0 scheme/test.scm
8. +1 −0 test.py
 @@ -0,0 +1,18 @@ +(define (erf x) + ; constants + (let* ((a1 0.254829592) + (a2 -0.284496736) + (a3 1.421413741) + (a4 -1.453152027) + (a5 1.061405429) + (p 0.3275911) + + ; Save the sign of x + (sign (if (< x 0) -1 1)) + (a (abs x)) + + ; A&S formula 7.1.26 + (t (/ 1.0 (+ 1.0 (* p a)))) + (y (- 1.0 (* (+ (* (+ (* (+ (* (+ (* a5 t) a4) t) a3) t) a2) t) a1) t (exp (* (- a) a)))))) + + (* sign y)))
 @@ -0,0 +1,4 @@ +(define (expm1 x) + (if (< (abs x) 1e-5) + (+ x (* 0.5 x x)) + (- (exp x) 1.0)))
 @@ -0,0 +1,121 @@ +; Visit http://www.johndcook.com/stand_alone_code.html for the source of this code and more like it. + +; Note that the functions Gamma and LogGamma are mutually dependent. + +(define (gamma x) + + (if (<= x 0) + (error "Invalid input")) + + ; Split the function domain into three intervals: + ; (0, 0.001), [0.001, 12), and (12, infinity) + + (cond + + ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; + ; First interval: (0, 0.001) + ; + ; For small x, 1/Gamma(x) has power series x + gamma x^2 - ... + ; So in this range, 1/Gamma(x) = x + gamma x^2 with error on the order of x^3. + ; The relative error over this interval is less than 6e-7. + + ((< x 0.001) + (let ((gamma 0.577215664901532860606512090)) ; Euler's gamma constant + (/ 1.0 (* x (+ 1.0 (* gamma x)))))) + + ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; + ; Second interval: [0.001, 12) + + ((< x 12.0) + ; The algorithm directly approximates gamma over (1,2) and uses + ; reduction identities to reduce other arguments to this interval. + + (letrec ((gamma-iter (lambda (z num den ps qs) + (if (null? ps) + (list num den) + (let ((new-num (* (+ num (car ps)) z)) + (new-den (+ (* den z) (car qs)))) + (gamma-iter z new-num new-den (cdr ps) (cdr qs)))))) + (gamma-z-n (lambda (result y n) + (if (= n 0) + result + (gamma-z-n (* result y) (+ y 1) (- n 1)))))) + + (let* ((arg-was-less-than-one (< x 1.0)) + ; Add or subtract integers as necessary to bring y into (1,2) + ; Will correct for this below + (n (if arg-was-less-than-one 0 (- (truncate x) 1))) + (y (if arg-was-less-than-one (+ x 1.0) (- x n))) + + ; numerator coefficients for approximation over the interval (1,2) + (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) + (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)) + + (z (- y 1)) + (a (gamma-iter z 0.0 1.0 p q)) + (num (car a)) + (den (cadr a)) + (result (+ (/ num den) 1.0))) + + ; Apply correction if argument was not initially in (1,2) + (if arg-was-less-than-one + (/ result (- y 1.0)) + (gamma-z-n result y n))))) + + ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; + ; Third interval: [12, infinity) + + ((<= x 171.624) + (exp (log-gamma x))) + + (else + ; Correct answer too large to display. + (inf)))) + +(define (log-gamma x) + (if (<= x 0) + (error "Invalid input")) + + (if (< x 12.0) + (log (abs (gamma x))) + + ; 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 + + (letrec ((log-gamma-iter (lambda (z sum c) + (if (null? c) + sum + (let ((s (+ (* sum z) (car c)))) + (log-gamma-iter z s (cdr c))))))) + + (let* ((c (reverse (list (/ 1.0 12.0) + (/ -1.0 360.0) + (/ 1.0 1260.0) + (/ -1.0 1680.0) + (/ 1.0 1188.0) + (/ -691.0 360360.0) + (/ 1.0 156.0) + (/ -3617.0 122400.0)))) + (z (/ 1.0 (* x x))) + (sum (log-gamma-iter z 0 c)) + (series (/ sum x)) + (half-log-two-pi 0.91893853320467274178032973640562)) + (+ (* (- x 0.5) (log x)) (- x) half-log-two-pi series)))))