1 ;;; Airy functions Ai(z) and Bi(z) - A&S 10.4
3 ;;; airy_ai(z) - Airy function Ai(z)
4 ;;; airy_dai(z) - Derivative of Airy function Ai(z)
5 ;;; airy_bi(z) - Airy function Bi(z)
6 ;;; airy_dbi(z) - Derivative of Airy function Bi(z)
8 ;;;; Copyright (C) 2005 David Billinghurst
10 ;;;; airy.lisp is free software; you can redistribute it
11 ;;;; and/or modify it under the terms of the GNU General Public
12 ;;;; License as published by the Free Software Foundation; either
13 ;;;; version 2, or (at your option) any later version.
15 ;;;; airy.lisp is distributed in the hope that it will be
16 ;;;; useful, but WITHOUT ANY WARRANTY; without even the implied
17 ;;;; warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
18 ;;;; See the GNU General Public License for more details.
20 ;;;; You should have received a copy of the GNU General Public License
21 ;;;; along with command-line.lisp; see the file COPYING. If not,
22 ;;;; write to the Free Software Foundation, Inc., 59 Temple Place -
23 ;;;; Suite 330, Boston, MA 02111-1307, USA.
27 (declaim (special *flonum-op
*))
32 (simplify (list '(%airy_ai
) (resimplify z
))))
34 (defprop $airy_ai %airy_ai alias
)
35 (defprop $airy_ai %airy_ai verb
)
36 (defprop %airy_ai $airy_ai reversealias
)
37 (defprop %airy_ai $airy_ai noun
)
38 (defprop %airy_ai simp-%airy_ai operators
)
39 (defprop %airy_ai simplim%airy_ai simplim%function
)
40 (defprop %airy_ai
((z) ((%airy_dai
) z
)) grad
)
42 ;; airy_ai distributes over lists, matrices, and equations
43 (defprop %airy_ai
(mlist $matrix mequal
) distribute_over
)
45 ;; airy_ai has mirror symmetry
46 (defprop %airy_ai t commutes-with-conjugate
)
49 ;; http://functions.wolfram.com/03.05.21.0002.01
50 ;; (z/(3^(2/3)*gamma(2/3)))*hypergeometric([1/3],[2/3,4/3],z^3/9)
51 ;; - (3^(1/6)/(4*%pi))*z^2*gamma(2/3)*hypergeometric([2/3],[4/3,5/3],z^3/9);
56 ((mexpt) 3 ((rat) -
2 3))
57 ((mexpt) ((%gamma
) ((rat) 2 3)) -
1)
60 ((mlist) ((rat) 2 3) ((rat) 4 3))
61 ((mtimes) ((rat) 1 9) ((mexpt) z
3)))
64 ((rat) -
1 4) ((mexpt) 3 ((rat) 1 6)) ((mexpt) $%pi -
1) ((%gamma
) ((rat) 2 3))
67 ((mlist) ((rat) 4 3) ((rat) 5 3))
68 ((mtimes) ((rat) 1 9) ((mexpt) z
3)))
73 (cond ((floatp z
) (airy-ai-real z
))
74 ((complexp z
) (airy-ai-complex z
))))
76 (setf (gethash '%airy_ai
*flonum-op
*) #'airy-ai
)
78 (defun simplim%airy_ai
(expr var val
)
79 ; Look for the limit of the argument
80 (let ((z (limit (cadr expr
) var val
'think
)))
81 (cond ((or (eq z
'$inf
) ; A&S 10.4.59
82 (eq z
'$minf
)) ; A&S 10.4.60
85 ; Handle other cases with the function simplifier
86 (simplify (list '(%airy_ai
) z
))))))
88 (defun simp-%airy_ai
(form unused x
)
89 (declare (ignore unused
))
91 (let ((z (simpcheck (cadr form
) x
)))
92 (cond ((equal z
0) ; A&S 10.4.4: Ai(0) = 3^(-2/3)/gamma(2/3)
94 ((mexpt simp
) 3 ((rat simp
) -
2 3))
95 ((mexpt simp
) ((%gamma simp
) ((rat simp
) 2 3)) -
1)))
96 ((flonum-eval (mop form
) z
))
97 (t (eqtest (list '(%airy_ai
) z
) form
)))))
100 ;; Derivative dAi/dz of Airy function Ai(z)
101 (defmfun $airy_dai
(z)
102 "Derivative dAi/dz of Airy function Ai(z)"
103 (simplify (list '(%airy_dai
) (resimplify z
))))
105 (defprop $airy_dai %airy_dai alias
)
106 (defprop $airy_dai %airy_dai verb
)
107 (defprop %airy_dai $airy_dai reversealias
)
108 (defprop %airy_dai $airy_dai noun
)
109 (defprop %airy_dai simp-%airy_dai operators
)
110 (defprop %airy_dai simplim%airy_dai simplim%function
)
111 (defprop %airy_dai
((z) ((mtimes) z
((%airy_ai
) z
))) grad
)
112 (defprop %airy_dai
((z) ((%airy_ai
) z
)) integral
)
114 ;; airy_dai distributes over lists, matrices, and equations
115 (defprop %airy_dai
(mlist $matrix mequal
) distribute_over
)
117 ;; airy_dai has mirror symmetry
118 (defprop %airy_dai t commutes-with-conjugate
)
121 (cond ((floatp z
) (airy-dai-real z
))
122 ((complexp z
) (airy-dai-complex z
))))
124 (setf (gethash '%airy_dai
*flonum-op
*) #'airy-dai
)
126 (defun simplim%airy_dai
(expr var val
)
127 ; Look for the limit of the argument
128 (let ((z (limit (cadr expr
) var val
'think
)))
129 (cond ((eq z
'$inf
) ; A&S 10.4.61
131 ((eq z
'$minf
) ; A&S 10.4.62
134 ; Handle other cases with the function simplifier
135 (simplify (list '(%airy_dai
) z
))))))
137 (defun simp-%airy_dai
(form unused x
)
138 (declare (ignore unused
))
140 (let ((z (simpcheck (cadr form
) x
)))
141 (cond ((equal z
0) ; A&S 10.4.5: Ai'(0) = -3^(-1/3)/gamma(1/3)
143 ((mexpt simp
) 3 ((rat simp
) -
1 3))
144 ((mexpt simp
) ((%gamma simp
) ((rat simp
) 1 3)) -
1)))
145 ((flonum-eval (mop form
) z
))
146 (t (eqtest (list '(%airy_dai
) z
) form
)))))
149 (defmfun $airy_bi
(z)
150 "Airy function Bi(z)"
151 (simplify (list '(%airy_bi
) (resimplify z
))))
153 (defprop $airy_bi %airy_bi alias
)
154 (defprop $airy_bi %airy_bi verb
)
155 (defprop %airy_bi $airy_bi reversealias
)
156 (defprop %airy_bi $airy_bi noun
)
157 (defprop %airy_bi simp-%airy_bi operators
)
158 (defprop %airy_bi simplim%airy_bi simplim%function
)
159 (defprop %airy_bi
((z) ((%airy_dbi
) z
)) grad
)
161 ;; airy_bi distributes over lists, matrices, and equations
162 (defprop %airy_bi
(mlist $matrix mequal
) distribute_over
)
164 ;; airy_bi has mirror symmetry
165 (defprop %airy_bi t commutes-with-conjugate
)
168 ;; http://functions.wolfram.com/03.06.21.0002.01
169 ;; (z/(3^(1/6)*gamma(2/3)))*hypergeometric([1/3],[2/3,4/3],z^3/9)
170 ;; + (3^(2/3)/(4*%pi))*z^2*gamma(2/3)*hypergeometric([2/3],[4/3,5/3],z^3/9);
175 ((mexpt) 3 ((rat) -
1 6))
176 ((mexpt) ((%gamma
) ((rat) 2 3)) -
1)
178 ((mlist) ((rat) 1 3))
179 ((mlist) ((rat) 2 3) ((rat) 4 3))
180 ((mtimes) ((rat) 1 9) ((mexpt) z
3)))
183 ((rat) 1 4) ((mexpt) 3 ((rat) 2 3)) ((mexpt) $%pi -
1) ((%gamma
) ((rat) 2 3))
185 ((mlist) ((rat) 2 3))
186 ((mlist) ((rat) 4 3) ((rat) 5 3))
187 ((mtimes) ((rat) 1 9) ((mexpt) z
3)))
192 (cond ((floatp z
) (airy-bi-real z
))
193 ((complexp z
) (airy-bi-complex z
))))
195 (setf (gethash '%airy_bi
*flonum-op
*) #'airy-bi
)
197 (defun simplim%airy_bi
(expr var val
)
198 ; Look for the limit of the argument
199 (let ((z (limit (cadr expr
) var val
'think
)))
200 (cond ((eq z
'$inf
) ; A&S 10.4.63
202 ((eq z
'$minf
) ; A&S 10.4.64
205 ; Handle other cases with the function simplifier
206 (simplify (list '(%airy_bi
) z
))))))
208 (defun simp-%airy_bi
(form unused x
)
209 (declare (ignore unused
))
211 (let ((z (simpcheck (cadr form
) x
)))
212 (cond ((equal z
0) ; A&S 10.4.4: Bi(0) = sqrt(3) 3^(-2/3)/gamma(2/3)
214 ((mexpt simp
) 3 ((rat simp
) -
1 6))
215 ((mexpt simp
) ((%gamma simp
) ((rat simp
) 2 3)) -
1)))
216 ((flonum-eval (mop form
) z
))
217 (t (eqtest (list '(%airy_bi
) z
) form
)))))
219 ;; Derivative dBi/dz of Airy function Bi(z)
220 (defmfun $airy_dbi
(z)
221 "Derivative dBi/dz of Airy function Bi(z)"
222 (simplify (list '(%airy_dbi
) (resimplify z
))))
224 (defprop $airy_dbi %airy_dbi alias
)
225 (defprop $airy_dbi %airy_dbi verb
)
226 (defprop %airy_dbi $airy_dbi reversealias
)
227 (defprop %airy_dbi $airy_dbi noun
)
228 (defprop %airy_dbi simp-%airy_dbi operators
)
229 (defprop %airy_dbi simplim%airy_dbi simplim%function
)
230 (defprop %airy_dbi
((z) ((mtimes) z
((%airy_bi
) z
))) grad
)
231 (defprop %airy_dbi
((z) ((%airy_bi
) z
)) integral
)
233 ;; airy_dbi distributes over lists, matrices, and equations
234 (defprop %airy_dbi
(mlist $matrix mequal
) distribute_over
)
236 ;; airy_dbi has mirror symmetry
237 (defprop %airy_dbi t commutes-with-conjugate
)
240 (cond ((floatp z
) (airy-dbi-real z
))
241 ((complexp z
) (airy-dbi-complex z
))))
243 (setf (gethash '%airy_dbi
*flonum-op
*) #'airy-dbi
)
245 (defun simplim%airy_dbi
(expr var val
)
246 ; Look for the limit of the argument
247 (let ((z (limit (cadr expr
) var val
'think
)))
248 (cond ((eq z
'$inf
) ; A&S 10.4.66
250 ((eq z
'$minf
) ; A&S 10.4.67
253 ; Handle other cases with the function simplifier
254 (simplify (list '(%airy_dbi
) z
))))))
256 (defun simp-%airy_dbi
(form unused x
)
257 (declare (ignore unused
))
259 (let ((z (simpcheck (cadr form
) x
)))
260 (cond ((equal z
0) ; A&S 10.4.5: Bi'(0) = sqrt(3) 3^(-1/3)/gamma(1/3)
262 ((mexpt simp
) 3 ((rat simp
) 1 6))
263 ((mexpt simp
) ((%gamma simp
) ((rat simp
) 1 3)) -
1)))
264 ((flonum-eval (mop form
) z
))
265 (t (eqtest (list '(%airy_dbi
) z
) form
)))))
267 ;; Numerical routines using slatec functions
269 (defun airy-ai-real (z)
270 " Airy function Ai(z) for real z"
271 (declare (type flonum z
))
272 ;; slatec:dai issues underflow warning for z > zmax. See dai.{f,lisp}
273 ;; This value is correct for IEEE double precision
274 (let ((zmax 92.5747007268))
275 (declare (type flonum zmax
))
276 (if (< z zmax
) (slatec:dai z
) 0.0)))
278 (defun airy-ai-complex (z)
279 "Airy function Ai(z) for complex z"
280 (declare (type (complex flonum
) z
))
281 (multiple-value-bind (var-0 var-1 var-2 var-3 air aii nz ierr
)
282 (slatec:zairy
(realpart z
) (imagpart z
) 0 1 0.0 0.0 0 0)
283 (declare (type flonum air aii
)
284 (type f2cl-lib
:integer4 nz ierr
)
285 (ignore var-0 var-1 var-2 var-3
))
286 ;; Check nz and ierr for errors
287 (if (and (= nz
0) (= ierr
0)) (complex air aii
) nil
)))
289 (defun airy-dai-real (z)
290 "Derivative dAi/dz of Airy function Ai(z) for real z"
291 (declare (type flonum z
))
292 (let ((rz (sqrt (abs z
)))
293 (c (* 2/3 (expt (abs z
) 3/2))))
294 (declare (type flonum rz c
))
295 (multiple-value-bind (var-0 var-1 var-2 ai dai
)
296 (slatec:djairy z rz c
0.0 0.0)
297 (declare (ignore var-0 var-1 var-2 ai
))
300 (defun airy-dai-complex (z)
301 "Derivative dAi/dz of Airy function Ai(z) for complex z"
302 (declare (type (complex flonum
) z
))
303 (multiple-value-bind (var-0 var-1 var-2 var-3 air aii nz ierr
)
304 (slatec:zairy
(realpart z
) (imagpart z
) 1 1 0.0 0.0 0 0)
305 (declare (type flonum air aii
)
306 (type f2cl-lib
:integer4 nz ierr
)
307 (ignore var-0 var-1 var-2 var-3
))
308 ;; Check nz and ierr for errors
309 (if (and (= nz
0) (= ierr
0)) (complex air aii
) nil
)))
311 (defun airy-bi-real (z)
312 "Airy function Bi(z) for real z"
313 (declare (type flonum z
))
314 ;; slatec:dbi issues overflows for z > zmax. See dbi.{f,lisp}
315 ;; This value is correct for IEEE double precision
316 (let ((zmax 104.2179765192136))
317 (declare (type flonum zmax
))
318 (if (< z zmax
) (slatec:dbi z
) nil
)))
320 (defun airy-bi-complex (z)
321 "Airy function Bi(z) for complex z"
322 (declare (type (complex flonum
) z
))
323 (multiple-value-bind (var-0 var-1 var-2 var-3 bir bii ierr
)
324 (slatec:zbiry
(realpart z
) (imagpart z
) 0 1 0.0 0.0 0)
325 (declare (type flonum bir bii
)
326 (type f2cl-lib
:integer4 ierr
)
327 (ignore var-0 var-1 var-2 var-3
))
328 ;; Check ierr for errors
329 (if (= ierr
0) (complex bir bii
) nil
)))
331 (defun airy-dbi-real (z)
332 "Derivative dBi/dz of Airy function Bi(z) for real z"
333 (declare (type flonum z
))
334 ;; Overflows for z > zmax.
335 ;; This value is correct for IEEE double precision
336 (let ((zmax 104.1525))
337 (declare (type flonum zmax
))
339 (let ((rz (sqrt (abs z
)))
340 (c (* 2/3 (expt (abs z
) 3/2))))
341 (declare (type flonum rz c
))
342 (multiple-value-bind (var-0 var-1 var-2 bi dbi
)
343 (slatec:dyairy z rz c
0.0 0.0)
344 (declare (type flonum bi dbi
)
345 (ignore var-0 var-1 var-2 bi
))
347 ;; Will overflow. Return unevaluated.
350 (defun airy-dbi-complex (z)
351 "Derivative dBi/dz of Airy function Bi(z) for complex z"
352 (declare (type (complex flonum
) z
))
353 (multiple-value-bind (var-0 var-1 var-2 var-3 bir bii ierr
)
354 (slatec:zbiry
(realpart z
) (imagpart z
) 1 1 0.0 0.0 0)
355 (declare (type flonum bir bii
)
356 (type f2cl-lib
:integer4 ierr
)
357 (ignore var-0 var-1 var-2 var-3
))
358 ;; Check ierr for errors
359 (if (= ierr
0) (complex bir bii
) nil
)))
