Newer
Older
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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
;********************************************************
; file: polynomial.lisp
; description: Polynomials
; date: Tue Mar 21 2006 - 18:33
; author: Liam M. Healy
; modified: Tue Mar 21 2006 - 23:07
;********************************************************
;;; $Id: $
(in-package :gsl)
;;; To do:
;;; Awaits incorporation of arrays into CFFI.
;;; Divided difference needs to be checked.
;;; Quadratic solver always gives output range warning, except when it needs to.
;;; Cubic solver always gives output range warning.
;;; finish
;;;;****************************************************************************
;;;; Polynomial Evaluation
;;;;****************************************************************************
(export 'polynomial-eval)
(defun polynomial-eval (coefficients x)
"Evaluate the polyonomial with coefficients at the point x."
(let ((len (length coefficients)))
(cffi::with-foreign-array (coef coefficients :double (list len))
(cffi:foreign-funcall
"gsl_poly_eval"
:pointer coef
:int len
:double x
:double))))
;;;;****************************************************************************
;;;; Divided Difference Representation of Polynomials
;;;;****************************************************************************
;;; Use with-divided-difference to compute the divided difference, an
;;; opaque object that is then passed to
;;; get-divided-difference, eval-divided-difference
;;; in the body.
(export
'(with-divided-difference get-divided-difference
eval-divided-difference taylor-divided-difference))
(defmacro with-divided-difference ((dd xa ya) &body body)
"Compute the divided difference and bind dd to it.
This variable may then be passed to divided-difference
functions in the body."
`(let ((,dd (find-divided-difference ,xa ,ya)))
(unwind-protect
(progn ,@body)
(cffi::foreign-array-free (first ,dd))
(cffi::foreign-array-free (second ,dd)))))
;;; Do not call this function directly; use with-divided-difference.
(defun find-divided-difference (xa ya)
"Find the divided difference representation of the interpolating polynomial
for the points stored in the arrays @var{xa} and @var{ya}."
(let* ((size (length xa))
(dd (foreign-alloc :double :count size))
(xac (foreign-alloc :double :count size)))
(cffi::with-foreign-array (yac ya :double (list size))
;; return code not checked
(cffi:foreign-funcall
"gsl_poly_dd_init"
:pointer dd
:pointer xac
:pointer yac
:uint size ; "size_t is almost always an unsigned int"
:int))
(list dd xac size)))
(defun get-divided-difference (dd)
"Convert the divided difference into an array."
(cffi::foreign-array-to-lisp
(first dd) :double (list (third dd))))
(defun eval-divided-difference (dd x)
"Evaluate the polynomial stored in divided-difference form
at the point @var{x}."
(cffi:foreign-funcall
"gsl_poly_dd_eval"
:pointer (first dd)
:pointer (second dd)
:uint (third dd) ; "size_t is almost always an unsigned int"
:double x
:double))
(defun taylor-divided-difference (dd xp)
"Convert the divided-difference representation of a polynomial
to a Taylor expansion about the point xp."
(let ((cc (foreign-alloc :double :count (third dd)))
(workspace (foreign-alloc :double :count (third dd))))
(unwind-protect
(progn
;; Return value not checked.
(cffi:foreign-funcall
"gsl_poly_dd_taylor"
:pointer cc
:double xp
:pointer (first dd)
:pointer (second dd)
:uint (third dd) ; "size_t is almost always an unsigned int"
:pointer workspace
:int)
(cffi::foreign-array-to-lisp
cc :double (list (third dd))))
(cffi::foreign-array-free workspace)
(cffi::foreign-array-free cc))))
;;;;****************************************************************************
;;;; Quadratic Equations
;;;;****************************************************************************
;;; GSL seems to always give a output range error warning:
;;; WARNING:
;;; GSL condition ERANGE (2), Output range error in
;;; (SOLVE-QUADRATIC 1.0d0 2.0d0 1.0d0)
;;; Oddly, there is NO warning if a complex answer is generated
;;; from the real root solver.
(defun-sf solve-quadratic ((a :double) (b :double) (c :double))
"gsl_poly_solve_quadratic"
:documentation
"The real roots of the quadratic equation a x^2 + b x + c = 0.
Warning: if roots are not real, nonsense is produced and
no warning is given."
:return (:double :double))
(defun-sf solve-quadratic-complex ((a :double) (b :double) (c :double))
"gsl_poly_complex_solve_quadratic"
:documentation
"The complex roots of the quadratic equation a x^2 + b x + c = 0."
:return (gsl-complex gsl-complex))
;;;;****************************************************************************
;;;; Cubic Equations
;;;;****************************************************************************
;;; GSL seems to always give warning:
;;; WARNING:
;;; GSL condition EFAULT (3), Invalid pointer in
;;; (SOLVE-CUBIC -3.0d0 3.0d0 -1.0d0)
;;; GSL appears to give a different return code (EDOM) if there is
;;; only one real root.
;;; (solve-cubic -6.0d0 -13.0d0 42.0d0)
;;; -3.0d0
;;; 1.9999999999999996d0
;;; 7.0d0
;;; (solve-cubic -1.0d0 1.0d0 -1.0d0)
;;; 1.0d0
;;; 1.9999999999999996d0
;;; 7.0d0
(defun-sf solve-cubic ((a :double) (b :double) (c :double))
"gsl_poly_solve_cubic"
:documentation
"Find the real roots of the cubic equation, x^3 + a x^2 + b x + c = 0
with a leading coefficient of unity. The real roots are returned.
No special indication is given if there is only one real root."
:return (:double :double :double))
(defun-sf solve-cubic-complex ((a :double) (b :double) (c :double))
"gsl_poly_complex_solve_cubic"
:documentation
"Find the complex roots of the cubic equation, x^3 + a x^2 + b x + c = 0
with a leading coefficient of unity. The real roots are returned.
No special indication is given if there is only one real root."
:return (gsl-complex gsl-complex gsl-complex))