1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
|
(sc-include "macros/macros")
(pre-include "gsl/gsl_math.h")
(pre-include "gsl/gsl_blas.h")
(pre-include "gsl/gsl_randist.h")
(pre-include "gsl/gsl_roots.h")
(pre-include "gsl/gsl_sf.h")
(pre-define (SIGNUM x)
(if* (< x 0) -1 1))
;; This function exists so that guile code can access the value of
;; M_PI as provided by <gsl_math.h>.
(define (pi) (double)
"Return the value of pi."
(return M-PI))
(define (ln-volume-of-ball dimension) (double (unsigned int))
"Return the natural logarithm of the volume of the unit
ball. DIMENSION of 2 corresponds to a circle, 3 corresponds to a
sphere, etc."
(return (- (* 0.5 dimension M-LNPI)
(gsl-sf-lngamma (+ 1 (* 0.5 dimension))))))
(define (volume-of-ball dimension) (double (unsigned int))
"Return the volume of the unit ball of DIMENSION. DIMENSION of 2
corresponds to a circle, 3 corresponds to a sphere, etc."
(return (exp (ln-volume-of-ball dimension))))
(define (ln-surface-area-of-ball dimension) (double (unsigned int))
"Return the natural logarithm of the surface area of the unit
ball. DIMENSION of 2 corresponds to a circle, 3 corresponds to a
sphere, etc."
(return (+ (log dimension)
(ln-volume-of-ball dimension))))
(define (surface-area-of-ball dimension) (double (unsigned int))
"Return the surface area of the unit ball of DIMENSION. DIMENSION of
2 corresponds to a circle, 3 corresponds to a sphere, etc."
(return (* dimension (volume-of-ball dimension))))
(define (lower-incomplete-gamma s x) (double double double)
(return (* (gsl-sf-gamma s)
(gsl-sf-gamma-inc-P s x))))
(define (angle-between-vectors x y) (double (const gsl-vector*) (const gsl-vector*))
"Return the angle between vectors X and Y. The returned value is in
the range [0,pi]."
(declare dot-product double)
(gsl-blas-ddot x y (address-of dot-product))
;; TODO: Is this a valid floating point comparison?
(return (if* (= dot-product 0)
0
(acos (/ dot-product
(gsl-blas-dnrm2 x)
(gsl-blas-dnrm2 y))))))
(define (dot-product x y) (double (const gsl-vector*) (const gsl-vector*))
"Return the dot product of vectors X and Y."
(declare result double)
(gsl-blas-ddot x y (address-of result))
(return result))
(define (gaussian-pdf x) (double double)
"Return exp(-x^2/2) / sqrt(2*pi)"
(return (gsl-ran-gaussian-pdf x 1)))
(define (gaussian-cdf x) (double double)
"Return \\int_{-\\inf}^x gaussian-pdf(t) dt."
(return (* 0.5 (+ 1 (gsl-sf-erf (/ x M-SQRT2))))))
(define (rerror approx exact) (double double double)
"Return the relative error between approximate value APPROX and
exact value EXACT."
(return (fabs (- 1 (/ approx exact)))))
(sc-define-syntax* (with-root-fsolver solver solver-type function a b body ...)
(let ((solver-type (sc-gensym)))
`(begin
(let* ((,solver-type (const gsl-root-fsolver-type*) ,solver-type))
(with-alloc solver gsl-root-fsolver*
(gsl-root-fsolver-alloc ,solver-type)
gsl-root-fsolver-free
(gsl-root-fsolver-set solver ,function ,a ,b)
,@body)))))
(sc-define-syntax* (with-error-handler handler body ...)
(let ((old-handler (sc-gensym)))
`(begin
(let* ((,old-handler gsl-error-handler-t* (gsl-set-error-handler ,handler)))
,@body
(gsl-set-error-handler ,old-handler)))))
(pre-let* (BISECTION-EPSABS 0 BISECTION-EPSREL 1e-6 BISECTION-MAX-ITERATIONS 1000)
(define (bisection f a b) (double gsl-function* double double)
(declare solution double)
(define (error-handler reason file line gsl-errno) (void (const char*) (const char*) int int)
(fprintf stderr "Bisection error handler invoked.\n")
(fprintf stderr "f(%g) = %g\n" a (GSL-FN-EVAL f a))
(fprintf stderr "f(%g) = %g\n" b (GSL-FN-EVAL f b))
(fprintf stderr "gsl: %s:%d: ERROR: %s\n" file line reason)
(abort))
(with-root-fsolver solver gsl-root-fsolver-bisection f a b
(with-error-handler error-handler
(do-while (= (gsl-root-test-interval (gsl-root-fsolver-x-lower solver)
(gsl-root-fsolver-x-upper solver)
BISECTION-EPSABS
BISECTION-EPSREL)
GSL-CONTINUE)
(gsl-root-fsolver-iterate solver)))
(set solution (gsl-root-fsolver-root solver)))
(return solution))
(define (bisection-rlimit f a b) (double gsl-function* double double)
(let* ((sign int (SIGNUM (GSL-FN-EVAL f a))))
(for-i i BISECTION-MAX-ITERATIONS
(cond
((> (* sign (GSL-FN-EVAL f b)) 0)
(set* b 2))
(else (return b)))))
(fprintf stderr "Bisection bracketing failed\n")
(abort)))
|
