Fix bug #4307: partswitch affects op and operatorp
[maxima.git] / src / airy.lisp
blob041ba0b343851ecd2facd0bb1bae39748c929693
1 ;;; Airy functions Ai(z) and Bi(z) - A&S 10.4
2 ;;;
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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
24 (in-package :maxima)
26 ;; Airy Ai function
28 (defprop %airy_ai simplim%airy_ai simplim%function)
29 (defprop %airy_ai ((z) ((%airy_dai) z)) grad)
31 ;; airy_ai distributes over lists, matrices, and equations
32 (defprop %airy_ai (mlist $matrix mequal) distribute_over)
34 ;; airy_ai has mirror symmetry
35 (defprop %airy_ai t commutes-with-conjugate)
37 ;; Integral of Ai(z)
38 ;; http://functions.wolfram.com/03.05.21.0002.01
39 ;; (z/(3^(2/3)*gamma(2/3)))*hypergeometric([1/3],[2/3,4/3],z^3/9)
40 ;; - (3^(1/6)/(4*%pi))*z^2*gamma(2/3)*hypergeometric([2/3],[4/3,5/3],z^3/9);
41 (defprop %airy_ai
42 ((z)
43 ((mplus)
44 ((mtimes)
45 ((mexpt) 3 ((rat) -2 3))
46 ((mexpt) ((%gamma) ((rat) 2 3)) -1)
47 ((%hypergeometric)
48 ((mlist) ((rat) 1 3))
49 ((mlist) ((rat) 2 3) ((rat) 4 3))
50 ((mtimes) ((rat) 1 9) ((mexpt) z 3)))
52 ((mtimes)
53 ((rat) -1 4) ((mexpt) 3 ((rat) 1 6)) ((mexpt) $%pi -1) ((%gamma) ((rat) 2 3))
54 ((%hypergeometric)
55 ((mlist) ((rat) 2 3))
56 ((mlist) ((rat) 4 3) ((rat) 5 3))
57 ((mtimes) ((rat) 1 9) ((mexpt) z 3)))
58 ((mexpt) z 2))))
59 integral)
61 (defun airy-ai (z)
62 (cond ((floatp z) (airy-ai-real z))
63 ((complexp z) (airy-ai-complex z))))
65 (setf (gethash '%airy_ai *flonum-op*) 'airy-ai)
67 (defun simplim%airy_ai (expr var val)
68 ; Look for the limit of the argument
69 (let ((z (limit (cadr expr) var val 'think)))
70 (cond ((or (eq z '$inf) ; A&S 10.4.59
71 (eq z '$minf)) ; A&S 10.4.60
74 ; Handle other cases with the function simplifier
75 (simplify (list '(%airy_ai) z))))))
77 (defun airy-ai-hypergeometric (z)
78 "Returns the hypergeometric representation of Airy Ai"
79 ;; See http://functions.wolfram.com/03.05.26.0001.01 and
80 ;; https://fungrim.org/entry/01bbb6/:
82 ;; Ai(z) = Ai(0)*hypergeometric([],[2/3],z^3/9)
83 ;; + z*Ai'(0)*hypergeometric([],[4/3],z^3/9)
84 (add (mul (ftake '%airy_ai 0)
85 (ftake '%hypergeometric
86 (list '(mlist))
87 (list '(mlist) (div 2 3))
88 (div (power z 3)
89 9)))
90 (mul z
91 (ftake '%airy_dai 0)
92 (ftake '%hypergeometric
93 (list '(mlist))
94 (list '(mlist) (div 4 3))
95 (div (power z 3)
96 9)))))
98 (def-simplifier airy_ai (z)
99 (cond ((equal z 0) ; A&S 10.4.4: Ai(0) = 3^(-2/3)/gamma(2/3)
100 (div (power 3 (div -2 3))
101 (take '(%gamma) (div 2 3))))
102 ((flonum-eval (mop form) z))
103 ((or (bigfloat-numerical-eval-p z)
104 (complex-bigfloat-numerical-eval-p z))
105 ($rectform
106 ($bfloat (airy-ai-hypergeometric z))))
107 ($hypergeometric_representation
108 (airy-ai-hypergeometric z))
109 (t (give-up))))
112 ;; Derivative dAi/dz of Airy function Ai(z)
113 (defprop %airy_dai simplim%airy_dai simplim%function)
114 (defprop %airy_dai ((z) ((mtimes) z ((%airy_ai) z))) grad)
115 (defprop %airy_dai ((z) ((%airy_ai) z)) integral)
117 ;; airy_dai distributes over lists, matrices, and equations
118 (defprop %airy_dai (mlist $matrix mequal) distribute_over)
120 ;; airy_dai has mirror symmetry
121 (defprop %airy_dai t commutes-with-conjugate)
123 (defun airy-dai (z)
124 (cond ((floatp z) (airy-dai-real z))
125 ((complexp z) (airy-dai-complex z))))
127 (setf (gethash '%airy_dai *flonum-op*) 'airy-dai)
129 (defun simplim%airy_dai (expr var val)
130 ; Look for the limit of the argument
131 (let ((z (limit (cadr expr) var val 'think)))
132 (cond ((eq z '$inf) ; A&S 10.4.61
134 ((eq z '$minf) ; A&S 10.4.62
135 '$und)
137 ; Handle other cases with the function simplifier
138 (simplify (list '(%airy_dai) z))))))
140 (defun airy-dai-hypergeometric (z)
141 "Returns the hypergeometric representation of Ai'(x), the derivative
142 of the Airy function Ai(x)"
143 ;; See http://functions.wolfram.com/03.07.26.0001.01 and
144 ;; https://fungrim.org/entry/20e530/.
147 ;; Ai'(z) = Ai'(0)*hypergeometric([],[1/3],z^3/9)
148 ;; + z^2/2*Ai(0)*hypergeometric([],[5/3],z^3/9)
149 (add (mul (ftake '%airy_dai 0)
150 (ftake '%hypergeometric
151 (list '(mlist))
152 (list '(mlist) (div 1 3))
153 (div (power z 3)
154 9)))
155 (mul z z 1//2
156 (ftake '%airy_ai 0)
157 (ftake '%hypergeometric
158 (list '(mlist))
159 (list '(mlist) (div 5 3))
160 (div (power z 3)
161 9)))))
163 (def-simplifier airy_dai (z)
164 (cond ((equal z 0) ; A&S 10.4.5: Ai'(0) = -3^(-1/3)/gamma(1/3)
165 (div -1
166 (mul (power 3 (div 1 3))
167 (take '(%gamma) (div 1 3)))))
168 ((flonum-eval (mop form) z))
169 ((or (bigfloat-numerical-eval-p z)
170 (complex-bigfloat-numerical-eval-p z))
171 ($rectform
172 ($bfloat (airy-dai-hypergeometric z))))
173 ($hypergeometric_representation
174 (airy-dai-hypergeometric z))
175 (t (give-up))))
177 (defprop %airy_bi simplim%airy_bi simplim%function)
178 (defprop %airy_bi ((z) ((%airy_dbi) z)) grad)
180 ;; airy_bi distributes over lists, matrices, and equations
181 (defprop %airy_bi (mlist $matrix mequal) distribute_over)
183 ;; airy_bi has mirror symmetry
184 (defprop %airy_bi t commutes-with-conjugate)
186 ;; Integral of Bi(z)
187 ;; http://functions.wolfram.com/03.06.21.0002.01
188 ;; (z/(3^(1/6)*gamma(2/3)))*hypergeometric([1/3],[2/3,4/3],z^3/9)
189 ;; + (3^(2/3)/(4*%pi))*z^2*gamma(2/3)*hypergeometric([2/3],[4/3,5/3],z^3/9);
190 (defprop %airy_bi
191 ((z)
192 ((mplus)
193 ((mtimes)
194 ((mexpt) 3 ((rat) -1 6))
195 ((mexpt) ((%gamma) ((rat) 2 3)) -1)
196 ((%hypergeometric)
197 ((mlist) ((rat) 1 3))
198 ((mlist) ((rat) 2 3) ((rat) 4 3))
199 ((mtimes) ((rat) 1 9) ((mexpt) z 3)))
201 ((mtimes)
202 ((rat) 1 4) ((mexpt) 3 ((rat) 2 3)) ((mexpt) $%pi -1) ((%gamma) ((rat) 2 3))
203 ((%hypergeometric)
204 ((mlist) ((rat) 2 3))
205 ((mlist) ((rat) 4 3) ((rat) 5 3))
206 ((mtimes) ((rat) 1 9) ((mexpt) z 3)))
207 ((mexpt) z 2))))
208 integral)
210 (defun airy-bi (z)
211 (cond ((floatp z) (airy-bi-real z))
212 ((complexp z) (airy-bi-complex z))))
214 (setf (gethash '%airy_bi *flonum-op*) 'airy-bi)
216 (defun simplim%airy_bi (expr var val)
217 ; Look for the limit of the argument
218 (let ((z (limit (cadr expr) var val 'think)))
219 (cond ((eq z '$inf) ; A&S 10.4.63
220 '$inf)
221 ((eq z '$minf) ; A&S 10.4.64
224 ; Handle other cases with the function simplifier
225 (simplify (list '(%airy_bi) z))))))
227 (defun airy-bi-hypergeometric (z)
228 "Returns the hypergeometric representation of Airy Bi"
229 ;; See http://functions.wolfram.com/03.06.26.0001.01 and https://fungrim.org/entry/bd319e/
231 ;; Bi(z) = Bi(0)*hypergeometric([],[2/3],z^3/9)
232 ;; + z*Bi'(0)*hypergeometric([],[4/2],z^3/9)
233 (add (mul (ftake '%airy_bi 0)
234 (ftake '%hypergeometric
235 (list '(mlist))
236 (list '(mlist) (div 2 3))
237 (div (power z 3)
238 9)))
239 (mul z
240 (ftake '%airy_dbi 0)
241 (ftake '%hypergeometric
242 (list '(mlist))
243 (list '(mlist) (div 4 3))
244 (div (power z 3)
245 9)))))
247 (def-simplifier airy_bi (z)
248 (cond ((equal z 0) ; A&S 10.4.4: Bi(0) = sqrt(3) 3^(-2/3)/gamma(2/3)
249 (div (mul (power 3 1//2)
250 (power 3 (div -2 3)))
251 (take '(%gamma) (div 2 3))))
252 ((flonum-eval (mop form) z))
253 ((or (bigfloat-numerical-eval-p z)
254 (complex-bigfloat-numerical-eval-p z))
255 ($rectform
256 ($bfloat (airy-bi-hypergeometric z))))
257 ($hypergeometric_representation
258 (airy-bi-hypergeometric z))
259 (t (give-up))))
262 ;; Derivative dBi/dz of Airy function Bi(z)
263 (defprop %airy_dbi simplim%airy_dbi simplim%function)
264 (defprop %airy_dbi ((z) ((mtimes) z ((%airy_bi) z))) grad)
265 (defprop %airy_dbi ((z) ((%airy_bi) z)) integral)
267 ;; airy_dbi distributes over lists, matrices, and equations
268 (defprop %airy_dbi (mlist $matrix mequal) distribute_over)
270 ;; airy_dbi has mirror symmetry
271 (defprop %airy_dbi t commutes-with-conjugate)
273 (defun airy-dbi (z)
274 (cond ((floatp z) (airy-dbi-real z))
275 ((complexp z) (airy-dbi-complex z))))
277 (setf (gethash '%airy_dbi *flonum-op*) 'airy-dbi)
279 (defun simplim%airy_dbi (expr var val)
280 ; Look for the limit of the argument
281 (let ((z (limit (cadr expr) var val 'think)))
282 (cond ((eq z '$inf) ; A&S 10.4.66
283 '$inf)
284 ((eq z '$minf) ; A&S 10.4.67
285 '$und)
287 ; Handle other cases with the function simplifier
288 (simplify (list '(%airy_dbi) z))))))
290 (defun airy-dbi-hypergeometric (z)
291 "Returns the hypergeometric representation of Bi'(z), the derivative
292 of Airy Bi"
293 ;; See http://functions.wolfram.com/03.08.26.0001.01 and
294 ;; https://fungrim.org/entry/4d65e5/.
296 ;; Bi'(z) = Bi'(0)*hypergeometric([],[1/3],z^3/9)
297 ;; + z^2/2*Bi(0)*hypergeometric([],[5/3],z^3/9)
298 (add (mul (ftake '%airy_dbi 0)
299 (ftake '%hypergeometric
300 (list '(mlist))
301 (list '(mlist) (div 1 3))
302 (div (power z 3)
303 9)))
304 (mul z z 1//2
305 (ftake '%airy_bi 0)
306 (ftake '%hypergeometric
307 (list '(mlist))
308 (list '(mlist) (div 5 3))
309 (div (power z 3)
310 9)))))
312 (def-simplifier airy_dbi (z)
313 (cond ((equal z 0) ; A&S 10.4.5: Bi'(0) = sqrt(3) 3^(-1/3)/gamma(1/3)
314 (div (mul (power 3 1//2)
315 (power 3 (div -1 3)))
316 (take '(%gamma) (div 1 3))))
317 ((flonum-eval (mop form) z))
318 ((or (bigfloat-numerical-eval-p z)
319 (complex-bigfloat-numerical-eval-p z))
320 ($rectform
321 ($bfloat (airy-dbi-hypergeometric z))))
322 ($hypergeometric_representation
323 (airy-dbi-hypergeometric z))
324 (t (give-up))))
326 ;; Numerical routines using slatec functions
328 (defun airy-ai-real (z)
329 " Airy function Ai(z) for real z"
330 (declare (type flonum z))
331 ;; slatec:dai issues underflow warning for z > zmax. See dai.{f,lisp}
332 ;; This value is correct for IEEE double precision
333 (let ((zmax 92.5747007268))
334 (declare (type flonum zmax))
335 (if (< z zmax) (slatec:dai z) 0.0)))
337 (defun airy-ai-complex (z)
338 "Airy function Ai(z) for complex z"
339 (declare (type (complex flonum) z))
340 (multiple-value-bind (var-0 var-1 var-2 var-3 air aii nz ierr)
341 (slatec:zairy (realpart z) (imagpart z) 0 1 0.0 0.0 0 0)
342 (declare (type flonum air aii)
343 (type f2cl-lib:integer4 nz ierr)
344 (ignore var-0 var-1 var-2 var-3))
345 ;; Check nz and ierr for errors
346 (if (and (= nz 0) (= ierr 0)) (complex air aii) nil)))
348 (defun airy-dai-real (z)
349 "Derivative dAi/dz of Airy function Ai(z) for real z"
350 (declare (type flonum z))
351 (let ((rz (sqrt (abs z)))
352 (c (* 2/3 (expt (abs z) 3/2))))
353 (declare (type flonum rz c))
354 (multiple-value-bind (var-0 var-1 var-2 ai dai)
355 (slatec:djairy z rz c 0.0 0.0)
356 (declare (ignore var-0 var-1 var-2 ai))
357 dai)))
359 (defun airy-dai-complex (z)
360 "Derivative dAi/dz of Airy function Ai(z) for complex z"
361 (declare (type (complex flonum) z))
362 (multiple-value-bind (var-0 var-1 var-2 var-3 air aii nz ierr)
363 (slatec:zairy (realpart z) (imagpart z) 1 1 0.0 0.0 0 0)
364 (declare (type flonum air aii)
365 (type f2cl-lib:integer4 nz ierr)
366 (ignore var-0 var-1 var-2 var-3))
367 ;; Check nz and ierr for errors
368 (if (and (= nz 0) (= ierr 0)) (complex air aii) nil)))
370 (defun airy-bi-real (z)
371 "Airy function Bi(z) for real z"
372 (declare (type flonum z))
373 ;; slatec:dbi issues overflows for z > zmax. See dbi.{f,lisp}
374 ;; This value is correct for IEEE double precision
375 (let ((zmax 104.2179765192136))
376 (declare (type flonum zmax))
377 (if (< z zmax) (slatec:dbi z) nil)))
379 (defun airy-bi-complex (z)
380 "Airy function Bi(z) for complex z"
381 (declare (type (complex flonum) z))
382 (multiple-value-bind (var-0 var-1 var-2 var-3 bir bii ierr)
383 (slatec:zbiry (realpart z) (imagpart z) 0 1 0.0 0.0 0)
384 (declare (type flonum bir bii)
385 (type f2cl-lib:integer4 ierr)
386 (ignore var-0 var-1 var-2 var-3))
387 ;; Check ierr for errors
388 (if (= ierr 0) (complex bir bii) nil)))
390 (defun airy-dbi-real (z)
391 "Derivative dBi/dz of Airy function Bi(z) for real z"
392 (declare (type flonum z))
393 ;; Overflows for z > zmax.
394 ;; This value is correct for IEEE double precision
395 (let ((zmax 104.1525))
396 (declare (type flonum zmax))
397 (if (< z zmax)
398 (let ((rz (sqrt (abs z)))
399 (c (* 2/3 (expt (abs z) 3/2))))
400 (declare (type flonum rz c))
401 (multiple-value-bind (var-0 var-1 var-2 bi dbi)
402 (slatec:dyairy z rz c 0.0 0.0)
403 (declare (type flonum bi dbi)
404 (ignore var-0 var-1 var-2 bi))
405 dbi))
406 ;; Will overflow. Return unevaluated.
407 nil)))
409 (defun airy-dbi-complex (z)
410 "Derivative dBi/dz of Airy function Bi(z) for complex z"
411 (declare (type (complex flonum) z))
412 (multiple-value-bind (var-0 var-1 var-2 var-3 bir bii ierr)
413 (slatec:zbiry (realpart z) (imagpart z) 1 1 0.0 0.0 0)
414 (declare (type flonum bir bii)
415 (type f2cl-lib:integer4 ierr)
416 (ignore var-0 var-1 var-2 var-3))
417 ;; Check ierr for errors
418 (if (= ierr 0) (complex bir bii) nil)))