2
;; Copyright (C) 2000, 2001, 2003 Barton Willis
4
;; Maxima code for evaluating orthogonal polynomials listed in Chapter 22
5
;; of Abramowitz and Stegun (A & S).
7
;; Author: Barton Willis, University of Nebraska at Kearney (aka UNK).
8
;; License: see README. The user of this code assumes all risk
9
;; for its use. It has no warranty. If you don't know the meaning
10
;; of "no warranty," don't use this code. :-)
13
($put '$orthopoly 1.0 '$version)
15
;; A while loop taken from nset.lisp. Someday it should be placed
16
;; in a file with other Maxima macros and deleted from here and from nset.
18
(defmacro while (cond &rest body)
23
;; If the input can be evaluated to a floating point number (either
24
;; real or complex), convert the input to a Common Lisp complex number.
25
;; When the input doesn't evaluate to a float, return the input.
27
(defun maxima-to-lisp-complex-float (a)
28
(let* (($ratprint nil)
30
(br ($float ($realpart b)))
31
(bi ($float ($imagpart b))))
32
(cond ((and (numberp br) (numberp bi))
33
(if (= bi 0) (float br 1.0L0)
34
(complex (float br 1.0L0)
38
;; Return true iff a is a float or a complex number with either a
39
;; real or imaginary part that is a float.
43
(and (complexp a) (or (floatp (realpart a)) (floatp (imagpart a))))))
45
;; Convert a Lisp complex number into a Maxima complex number. When the
46
;; input isn't a Lisp complex number, return the argument.
48
(defun lisp-complex-to-maxima (x)
49
(if (complexp x) (add (realpart x) (mul '$%i (imagpart x))) x))
51
;; When orthopoly_returns_intervals is true, floating point evaluation
52
;; returns an interval using the form (($interval) c r), where c is the
53
;; center of the interval and r is its radius. We don't provide the user
54
;; with any tools for working with intervals; if a user wants
56
(defmvar $orthopoly_returns_intervals t)
58
(defun orthopoly-return-handler (d f e)
59
(cond ((or (floatp f) (complexp f))
60
(let ((df (maxima-to-lisp-complex-float d))
62
(setq f (lisp-complex-to-maxima
63
(if (numberp df) (* df f) (mul d f))))
64
(setq ef (if (and (numberp df) (numberp e)) (abs (* df e)) nil))
65
(cond ((and $orthopoly_returns_intervals (floatp ef))
66
`(($interval simp) ,f ,ef))
70
;; When a user requests the derivative of an a function in this package
71
;; with respect to the order or some other parameter, return a form
72
;; ((unk) input from user). We "simplify" this form by printing an error.
74
(defprop unk simp-unk operators)
76
(defun simp-unk (x y z)
77
(declare (ignore y z))
78
(merror "Maxima doesn't know the derivative of ~:M with respect the ~:M argument" (nth 2 x) (nth 1 x)))
80
;; A left continuous unit step function; thus
82
;; unit_step(x) = 0 for x <= 0 and 1 for x > 0.
84
;; This function differs from (1 + signum(x))/2 which isn't left or right
87
(defprop $unit_step simp-unit-step operators)
88
(defprop $unit_step "\\Theta" texword)
90
(defun simp-unit-step (a y z)
92
(setq y (simpcheck (cadr a) z))
94
(cond ((or (eq s '$nz) (eq s '$zero) (eq s '$neg)) 0)
96
(t `(($unit_step simp) ,y)))))
98
;; We barely support intervals; these functions aren't for user-level use.
99
;; The function intervalp returns true iff its argument is an interval.
101
(defun $intervalp (a)
102
(and (consp a) (consp (car a)) (eq (caar a) '$interval)))
104
;; Multiply x by a, where x isn't an interval and a might be an
105
;; interval. When a isn't an interval or x is an interval, return
106
;; x * a. The optional argument dx is an upper bound for the
107
;; relative error in x.
109
(defun interval-mult (x a &optional dx)
110
(if (null dx) (setq dx 0))
111
(cond ((and ($intervalp a) (not ($intervalp x)))
114
(if (or (floatp p) (floatp q))
115
(setq x ($float ($rectform x))))
116
(setq p ($expand (mult x p)))
117
(setq q ($expand (mult x (add q (simplify `((mabs) ,(mul p dx)))))))
118
(setq q (simplify `((mabs) ,q)))
119
`(($interval) ,p ,q)))
122
;; TeX a function with subscripts and superscripts. The string fn is the
123
;; function name, the list sub holds the positions of the subscripts, the list
124
;; sup holds the positions of the superscripts, and i is the position of the
125
;; function argument. When b1 is true, the subscript is surrounded by parens;
126
;; when b2 is true, the superscript is surrounded by parens. The lists sub and
127
;; sup may be nil, but the function must have at least on argument.
129
(defun tex-sub-and-super-scripted-function (fn sub b1 sup b2 i x l r)
131
(let ((lo) (hi) (s1) (s2))
132
(setq s1 (if b1 `("\\left(") nil))
133
(setq s2 (if b1 `("\\right)") nil))
135
(setq lo (cons (nth i x) lo)))
136
(setq lo (if lo (tex-list (nreverse lo) s1 s2 ",") nil))
137
(setq s1 (if b2 `("\\left(") nil))
138
(setq s2 (if b2 `("\\right)") nil))
140
(setq hi (cons (nth i x) hi)))
141
(setq hi (if hi (tex-list (nreverse hi) s1 s2 ",") nil))
143
(if lo `("_{",@lo "}") nil)
144
(if hi `("^{" ,@hi "}") nil)
145
`(,@(tex-list (nthcdr i x) `("\\left(") `("\\right)") ","))
148
(defun dimension-sub-and-super-scripted-function (fn sub sup b2 k x)
149
(let ((lo) (hi) (form))
151
(setq lo (cons (nth i x) lo)))
152
(setq lo (nreverse lo))
154
(setq hi (cons (nth i x) hi)))
155
(setq hi (nreverse hi))
157
(setq form `((,fn simp array) ,@lo)))
159
(setq form `((mexpt) ((,fn simp array) ,@lo) (("") ,@hi))))
161
(setq form `((mexpt) ((,fn simp array) ,@lo) ,@hi))))
162
`((,form simp) ,@(nthcdr k x))))
165
;; (1) each member of the list a evaluates to a floating
166
;; point number (using $float) and
167
;; (2) at least one member of a has a real or imaginary part that is a
168
;; floating point number *or* a bigfloat.
169
;; When we find an member of a that doesn't evaluate to a float,
170
;; immediately bail out and return nil.
172
(defun use-float (&rest a)
173
(let ((xr) (xi) (float-found nil) (okay t))
175
(setq x ($rectform x))
176
(setq xr ($realpart x))
177
(setq xi ($imagpart x))
178
(setq float-found (or float-found (floatp xr) (floatp xi)
179
($bfloatp xr) ($bfloatp xi)))
180
(setq okay (and (floatp ($float xr)) (floatp ($float xi))))
181
(if (not okay) (return)))
182
(and okay float-found)))
184
;; When n is a nonnegative integer, return 1F1(a,b,x); that is return
185
;; sum(pochhammer(a,k) * x^k /pochhammer(b,k),k,0,n). Regardless of the
186
;; value of n, when x = 0, immediately return 1.
188
;; Any orthopoly function that calls hypergeo should check that n is an
189
;; integer with n > -1; thus the summation code shouldn't get called for
190
;; any orthopoly function. It wouldn't be terrible if it happened -- I
191
;; think the summation form isn't useful and a user can get into trouble
192
;; with things like subst(0,x,sum(x^k,k,0,n)).
194
(defun $hypergeo11 (a b x n)
196
((and (integerp n) (> n -1))
197
(cond ((and (like a (- n)) (use-float b x))
199
(multiple-value-setq (f e)
200
(hypergeo11-float (- n) (maxima-to-lisp-complex-float b)
201
(maxima-to-lisp-complex-float x)))
204
(let ((sum 1) (cf 1))
206
(setq cf (div (mul cf (add i a) x) (mul (add i b) (+ 1 i))))
207
(setq sum (add cf sum)))))))
210
((mtimes) ((mexpt) ((mfactorial) ,$genindex) -1)
211
(($pochhammer) ,a ,$genindex)
212
((mexpt) (($pochhammer ) ,b ,$genindex) -1)
213
((mexpt) ,x ,$genindex)) ,$genindex 0 ,n))))
215
;; return the F[2,1] hypergeometic sum from 0 to n. Use genindex as the
216
;; sum index; genindex is defined in asum.lisp
218
(defun $hypergeo21 (a b c x n)
220
((and (integerp n) (> n -1))
221
(cond ((and (like a (- n)) (use-float b c x))
223
(multiple-value-setq (f e)
224
(hypergeo21-float (- n)
225
(maxima-to-lisp-complex-float b)
226
(maxima-to-lisp-complex-float c)
227
(maxima-to-lisp-complex-float x)))
230
(let ((sum 1) (cf 1))
232
(setq cf (div (mul cf (add i a) (add i b) x) (mul (+ 1 i)
234
(setq sum (add cf sum)))))))
238
((mtimes) (($pochhammer) ,a ,$genindex) (($pochhammer) ,b ,$genindex)
239
((mexpt) (($pochhammer) ,c ,$genindex) -1)
240
((mexpt) ((mfactorial) ,$genindex) -1)
241
((mexpt) ,x ,$genindex))
244
;; When n is an integer with n > -1, return the product
245
;; x (x + 1) (x + 2) ... (x + n - 1). This is the same as
246
;; gamma(x + n) / gamma(x). See A&S 6.1.22, page 256. When
247
;; n isn't an integer or n < 0, return the form (($pochhammer) x n).
248
;; Also notice that pochhammer(1,n) = n!.
250
;; We trap the case x is a complex float; otherwise, Maxima has to..
252
(defmvar $pochhammer_max_index 100)
254
(defun $pochhammer (x n)
256
(div (power -1 n) ($pochhammer (sub 1 x) (neg n))))
258
((and (integerp n) (use-float x))
260
(setq x (maxima-to-lisp-complex-float x))
261
(dotimes (i n (lisp-float-to-maxima-float acc))
262
(setq acc (* acc (+ i x))))))
264
((and (integerp n) (<= n $pochhammer_max_index))
267
(setq acc (mul acc (add i x))))))
269
(t `(($pochhammer simp) ,x ,n))))
271
(putprop '$pochhammer
273
((mtimes) (($pochhammer) x n)
274
((mplus) ((mtimes) -1 ((mqapply) (($psi array) 0) x))
275
((mqapply) (($psi array) 0) ((mplus) n x))))
276
((mtimes) (($pochhammer) x n)
277
((mqapply) (($psi array) 0) ((mplus) n x))))
280
(defprop $pochhammer tex-pochhammer tex)
282
(defun tex-pochhammer (x l r)
283
(setq x (mapcar #'(lambda (s) (tex s nil nil nil nil)) (cdr x)))
292
(setf (get '$pochhammer 'dimension) 'dimension-pochhammer)
294
(defun dimension-pochhammer (form result)
295
(setq form `(( (("") ,(nth 1 form)) simp array) ,(nth 2 form)))
296
(dimension-array form result))
298
;; pochhammer-quotient(a b x n) returns pochhammer( a,n) / pochhammer(b,n).
299
;; Only when one of a, b, or x is a floating point number and the others are
300
;; numeric types (that is when (use-float a b x) evaluates to true) does this
301
;; function differ from using (div ($pochhammer a n) ($pochhammer b n)).
302
;; In the floating point case, we arrange the operations to reduce rounding
303
;; errors and to reduce over and underflows.
305
(defun pochhammer-quotient (a b x n)
307
(pochhammer-quotient b a x (neg n)))
308
((and (integerp n) (use-float a b x))
310
(setq a (maxima-to-lisp-complex-float a))
311
(setq b (maxima-to-lisp-complex-float b))
312
(dotimes (i n (lisp-float-to-maxima-float acc))
313
(setq acc (* acc (/ (+ i a) (+ i b)))))))
314
(t (div ($pochhammer a n) ($pochhammer b n)))))
316
(defun use-hypergeo (n x)
318
(or (and (integerp n) (> n -1))
319
($featurep n '$integer)))
320
; (and ($featurep n '$integer) ($ratnump x))))
322
;; See A&S 22.5.42, page 779. For integer n, use the identity
323
;; binomial(n+a,n) = pochhammer(a+1,n)/pochhammer(1,n)
325
;; The condition number of the pochhammer function is
326
;; |1 + x/(x+1) + x/(x+2) + ... + x/(x+n-1)| <= n.
327
;; The relative floating point error in computing the pochhammer
328
;; function is bounded by 3 n eps. Putting these errors together,
329
;; the error in d is bounded 4 n |d| eps.
331
;; We're looking at (d +|- 4 n eps |d|)(f +|- e) = d (f +|- (e + 4 n eps |f|)) +
334
(defun $jacobi_p (n a b x)
335
(cond ((use-hypergeo n x)
337
;(setq d (div ($pochhammer (add a 1) n) ($pochhammer 1 n)))
338
(setq d (pochhammer-quotient (add a 1) 1 x n))
339
(multiple-value-setq (f e)
340
($hypergeo21 (mul -1 n) (add n a b 1) (add a 1)
341
(div (add 1 (mul -1 x)) 2) n))
342
(setq e (if e (+ e (* 4 n (abs f) double-float-epsilon)) nil))
343
(orthopoly-return-handler d f e)))
344
(t `(($jacobi_p simp) ,n ,a ,b ,x))))
348
((unk) "$first" "$jacobi_p")
349
((unk) "$second" "$jacobi_p")
350
((unk) "$third" "$jacobi_p")
353
((mexpt) ((mplus ) a b ((mtimes ) 2 n)) -1)
357
((mplus) a n) ((mplus) b n)
358
(($jacobi_p) ((mplus) -1 n) a b x))
359
((mtimes) n (($jacobi_p) n a b x)
360
((mplus) a ((mtimes ) -1 b)
362
((mplus) ((mtimes) -1 a) ((mtimes ) -1 b)
363
((mtimes) -2 n)) x))))
364
((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt ) x 2))) -1)))
367
(defprop $jacobi_p tex-jacobi-poly tex)
369
(defun tex-jacobi-poly (x l r)
370
(tex-sub-and-super-scripted-function "P" `(0) nil `(1 2) t 3 x l r))
372
(setf (get '$jacobi_p 'dimension) 'dimension-jacobi-p)
374
(defun dimension-jacobi-p (form result)
376
(dimension-sub-and-super-scripted-function "P" `(1) `(2 3) t 4 form)
379
;; See A&S 22.5.46, page 779.
381
(defun $ultraspherical (n a x)
382
(cond ((use-hypergeo n x)
384
;(setq d (div ($pochhammer (mul 2 a) n) ($pochhammer 1 n)))
385
(setq d (pochhammer-quotient (mul 2 a) 1 x n))
386
(multiple-value-setq (f e)
387
($hypergeo21 (mul -1 n) (add n (mul 2 a)) (add a (div 1 2))
388
(div (add 1 (mul -1 x)) 2) n))
389
(setq e (if e (+ e (* 4 n (abs f) double-float-epsilon)) nil))
390
(orthopoly-return-handler d f e)))
391
(t `(($ultraspherical simp) ,n ,a ,x))))
393
(putprop '$ultraspherical
395
((unk) "$first" "$ultraspherical")
396
((unk) "$second" "$ultrapsherical")
401
((mplus) -1 ((mtimes) 2 a) n)
402
(($ultraspherical) ((mplus) -1 n) a x))
403
((mtimes) -1 n (($ultraspherical) n a x) x))
404
((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))) -1)))
407
(defprop $ultraspherical tex-ultraspherical tex)
409
(defun tex-ultraspherical (x l r)
410
(tex-sub-and-super-scripted-function "C" `(0) nil `(1) t 2 x l r))
412
(setf (get '$ultraspherical 'dimension) 'dimension-ultraspherical)
414
(defun dimension-ultraspherical (form result)
416
(dimension-sub-and-super-scripted-function "C" `(1) `(2) t 3 form)
419
(defun $chebyshev_t (n x)
420
(cond ((use-hypergeo n x)
422
(multiple-value-setq (f e)
423
($hypergeo21 (mul -1 n) n (rat 1 2) (div (add 1 (mul -1 x)) 2) n))
424
(orthopoly-return-handler 1 f e)))
425
(t `(($chebyshev_t simp) ,n ,x))))
427
(putprop '$chebyshev_t
429
((unk) "$first" "$chebyshev_t")
432
((mtimes) n (($chebyshev_t) ((mplus ) -1 n) x))
433
((mtimes ) -1 n (($chebyshev_t) n x) x))
434
((mexpt) ((mplus ) 1 ((mtimes) -1 ((mexpt) x 2))) -1)))
437
(defprop $chebyshev_t tex-chebyshev-t tex)
439
(defun tex-chebyshev-t (x l r)
440
(tex-sub-and-super-scripted-function "T" `(0) nil nil nil 1 x l r))
442
(setf (get '$chebyshev_t 'dimension) 'dimension-chebyshev-t)
444
(defun dimension-chebyshev-t (form result)
446
(dimension-sub-and-super-scripted-function "T" `(1) nil nil 2 form)
449
;; See A & S 22.5.48, page 779.
451
(defun $chebyshev_u (n x)
452
(cond ((use-hypergeo n x)
455
(multiple-value-setq (f e)
456
($hypergeo21 (mul -1 n) (add 2 n) (rat 3 2)
457
(div (add 1 (mul -1 x)) 2) n))
458
(orthopoly-return-handler d f e)))
459
(t `(($chebyshev_u simp) ,n ,x))))
461
(putprop '$chebyshev_u
463
((unk) "$first" "$chebyshev_u")
468
((mplus) 1 n) (($chebyshev_u) ((mplus) -1 n) x))
469
((mtimes) -1 n (($chebyshev_u) n x) x))
470
((mexpt) ((mplus ) 1 ((mtimes) -1 ((mexpt) x 2))) -1)))
473
(defprop $chebyshev_u tex-chebyshev-u tex)
475
(defun tex-chebyshev-u (x l r)
476
(tex-sub-and-super-scripted-function "U" `(0) nil nil nil 1 x l r))
478
(setf (get '$chebyshev_u 'dimension) 'dimension-chebyshev-u)
480
(defun dimension-chebyshev-u (form result)
482
(dimension-sub-and-super-scripted-function "U" `(1) nil nil 2 form)
485
;; See A & S 22.5.35, page 779. We evaluate the legendre polynomials
486
;; as jacobi_p(n,0,0,x). Eat less exercise more.
488
(defun $legendre_p (n x)
489
(cond ((use-hypergeo n x)
491
(t `(($legendre_p simp) ,n ,x))))
493
(putprop '$legendre_p
495
((unk) "$first" "$legendre_p")
498
((mtimes) n (($legendre_p) ((mplus) -1 n) x))
499
((mtimes) -1 n (($legendre_p) n x) x))
500
((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))) -1)))
503
(defprop $legendre_p tex-legendre-p tex)
505
(defun tex-legendre-p (x l r)
506
(tex-sub-and-super-scripted-function "P" `(0) nil nil nil 1 x l r))
508
(setf (get '$legendre_p 'dimension) 'dimension-legendre-p)
510
(defun dimension-legendre-p (form result)
512
(dimension-sub-and-super-scripted-function "P" `(1) nil nil 2 form)
515
(defun $legendre_q (n x)
516
(if (and (integerp n) (> n -1))
517
($assoc_legendre_q n 0 x)
518
`(($legendre_q simp) ,n ,x)))
520
(putprop '$legendre_q
522
((unk) "$first" "$legendre_p")
524
((mtimes) -1 (($kron_delta) 0 n)
525
((mexpt) ((mplus) -1 ((mexpt) x 2)) -1))
528
((mtimes) n (($legendre_q) ((mplus) -1 n) x))
529
((mtimes) -1 n (($legendre_q) n x) x))
530
((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))) -1))))
533
(defprop $legendre_q tex-legendre-q tex)
535
(defun tex-legendre-q (x l r)
536
(tex-sub-and-super-scripted-function "Q" `(0) nil nil nil 1 x l r))
538
(setf (get '$legendre_q 'dimension) 'dimension-legendre-q)
540
(defun dimension-legendre-q (form result)
542
(dimension-sub-and-super-scripted-function "Q" `(1) nil nil 2 form)
545
;; See A & S 8.6.7 and 8.2.6 pages 333 and 334. I chose the
546
;; definition that is real valued on (-1,1).
548
;; For negative m, A & S 8.2.6 page 333 and G & R 8.706 page 1000
549
;; disagree; the factor of exp(i m pi) in A & S 8.1.6 suggests to me that
550
;; A & S 8.2.6 is bogus. As further evidence, Macsyma 2.4 seems to
551
;; agree with G & R 8.706. I'll use G & R.
553
;; Return assoc_legendre(0,m,x). This isn't a user-level function; we
554
;; don't check that m is a positive integer.
558
(merror "function q0m called with negative order. File a bug report"))
560
(div (simplify `((%log) ,(div (add 1 x) (sub 1 x)))) 2))
562
(mul (factorial (- m 1)) `((rat simp) 1 2)
563
(if (oddp m) -1 1) (power (sub 1 (mult x x)) (div m 2))
565
(mul (if (oddp m) 1 -1) (power (add 1 x) (neg m)))
566
(power (sub 1 x) (neg m)))))))
568
;; Return assoc_legendre(1,m,x). This isn't a user-level function; we
569
;; don't check that m is a positive integer; we don't check that m is
570
;; a positive integer.
574
(merror "function q1m called with negative order. File a bug report"))
576
(sub (mul x (q0m 0 x)) 1))
578
(mul -1 (power (sub 1 (mult x x)) `((rat simp) 1 2))
579
(sub (q0m 0 x) (div x (sub (mul x x) 1)))))
581
(mul (if (oddp m) -1 1) (power (sub 1 (mult x x)) (div m 2))
583
(mul (factorial (- m 2)) `((rat simp) 1 2) (if (oddp m) -1 1)
584
(sub (power (add x 1) (sub 1 m))
585
(power (sub x 1) (sub 1 m))))
587
(mul (factorial (- m 1)) `((rat simp) 1 2) (if (oddp m) -1 1)
588
(add (power (add x 1) (neg m))
589
(power (sub x 1) (neg m)))))))))
591
;; Return assoc_legendre_q(n,n,x). I don't have a reference that gives
592
;; a formula for assoc_legendre_q(n,n,x). To figure one out, I used
593
;; A&S 8.2.1 and a formula for assoc_legendre_p(n,n,x).
595
;; After finishing the while loop, q = int(1/(1-t^2)^(n+1),t,0,x).
597
(defun assoc-legendre-q-nn (n x)
599
(z (sub 1 (mul x x)))
601
(setq q (div (simplify `((%log) ,(div (add 1 x) (sub 1 x)))) 2))
603
(setq q (add (mul (sub 1 `((rat) 1 ,(* 2 i))) q)
604
(div x (mul 2 i (power z i)))))
606
(mul (^ -2 n) (factorial n) (power z (div n 2)) q)))
608
;; Use degree recursion to find the assoc_legendre_q function.
609
;; See A&S 8.5.3. When i = m in the while loop, we have a special
610
;; case. For negative order, see A&S 8.2.6.
612
(defun $assoc_legendre_q (n m x)
613
(cond ((and (integerp n) (> n -1) (integerp m) (<= (abs m) n))
615
(mul (div (factorial (+ n m)) (factorial (- n m)))
616
($assoc_legendre_q n (- m) x)))
618
(if (not (or (floatp x) ($bfloatp x))) (setq x ($rat x)))
619
(let* ((q0 (q0m m x))
620
(q1 (if (= n 0) q0 (q1m m x)))
622
(use-rat (or ($ratp x) (floatp x) ($bfloatp x))))
625
(setq q (if (= i m) (assoc-legendre-q-nn i x)
626
(div (sub (mul (- (* 2 i) 1) x q1)
627
(mul (+ i -1 m) q0)) (- i m))))
631
(if use-rat q1 ($ratsimp q1))))))
632
(t `(($assoc_legendre_q simp) ,n ,m ,x))))
634
;; See G & R, 8.733, page 1005 first equation.
636
(putprop '$assoc_legendre_q
638
((unk) "$first" "$assoc_legendre_q")
639
((unk) "$second" "$assoc_legendre_q")
644
((mtimes) -1 ((mplus) 1 ((mtimes) -1 m) n)
645
(($assoc_legendre_q ) ((mplus ) 1 n) m x))
646
((mtimes) ((mplus) 1 n)
647
(($assoc_legendre_q ) n m x) x))
648
((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))) -1))))
651
(defprop $assoc_legendre_q tex-assoc-legendre-q tex)
653
(defun tex-assoc-legendre-q (x l r)
654
(tex-sub-and-super-scripted-function "Q" `(0) nil `(1) nil 2 x l r))
656
(setf (get '$assoc_legendre_q 'dimension) 'dimension-assoc-legendre-q)
658
(defun dimension-assoc-legendre-q (form result)
660
(dimension-sub-and-super-scripted-function "Q" `(1) `(2) nil 3 form)
663
;; See A & S 22.5.37 page 779, A & S 8.6.6 (second equation) page 334, and
664
;; A & S 8.2.5 page 333. For n < 0, see A&S 8.2.1 page 333.
666
(defun $assoc_legendre_p (n m x)
667
(let ((f) (d) (dx 0))
668
(cond ((and (integerp n) (integerp m))
673
(setq f ($assoc_legendre_p n (neg m) x))
674
(setq d (div (factorial (+ n m)) (factorial (- n m))))
681
`((%genfact) ,(- (* 2 m) 1) ,(- m 1) 2)))
682
(setq d (mul d (if (oddp m) -1 1)))
683
(setq d (mul d (power (sub 1 (mul x x)) (div m 2))))))
686
($ultraspherical (- n m) (add m (rat 1 2)) x)))))
689
(setq f `(($assoc_legendre_p simp) ,n ,m ,x))))
690
(interval-mult d f (* double-float-epsilon dx))))
693
;; For the derivative of the associated legendre p function, see
694
;; A & S 8.5.4 page 334.
696
(putprop `$assoc_legendre_p
698
((unk) "$first" "$assoc_legendre_p")
699
((unk) "$second" "$assoc_legendre_p")
702
((mtimes simp) -1 ((mplus simp) m n) (($unit_step) n)
703
(($assoc_legendre_p simp) ((mplus simp) -1 n) m x))
704
((mtimes simp) n (($assoc_legendre_p simp) n m x) x))
705
((mexpt simp) ((mplus simp) -1 ((mexpt simp) x 2)) -1)))
708
(defprop $assoc_legendre_p tex-assoc-legendre-p tex)
710
(defun tex-assoc-legendre-p (x l r)
711
(tex-sub-and-super-scripted-function "P" `(0) nil `(1) nil 2 x l r))
713
(setf (get '$assoc_legendre_p 'dimension) 'dimension-assoc-legendre-p)
715
(defun dimension-assoc-legendre-p (form result)
717
(dimension-sub-and-super-scripted-function "P" `(1) `(2) nil 3 form)
720
;; See A&S 22.5.55 and 22.5.56, page 780.
722
(defun $hermite (n x)
727
(multiple-value-setq (f e)
728
($hypergeo11 (- n) (rat 1 2) (mul x x) n))
729
(setq d (mul (if (oddp n) -1 1) (factorial (* 2 n))
730
(div 1 (factorial n)))))
732
(setq n (/ (- n 1) 2))
733
(multiple-value-setq (f e)
734
($hypergeo11 (- n) (rat 3 2) (mul x x) n))
735
(setq d (mul (if (oddp n) -1 1) (factorial (+ 1 (* 2 n))) 2 x
736
(div 1 (factorial n))))))
737
(orthopoly-return-handler d f e)))
738
(t `(($hermite) ,n ,x))))
742
((unk) "$first" "$hermite")
743
((mtimes) 2 n (($hermite) ((mplus) -1 n) x)))
746
(defprop $hermite tex-hermite tex)
748
(defun tex-hermite (x l r)
749
(tex-sub-and-super-scripted-function "H" `(0) nil nil nil 1 x l r))
751
(setf (get '$hermite 'dimension) 'dimension-hermite)
753
(defun dimension-hermite (form result)
755
(dimension-sub-and-super-scripted-function "H" `(1) nil nil 2 form)
758
;; See A & S 22.5.54, page 780. For integer n, use the identity
759
;; binomial(n+a,n) = pochhammer(a+1,n)/pochhammer(1,n)
761
(defun $gen_laguerre (n a x)
762
(cond ((use-hypergeo n x)
764
;(setq d (div ($pochhammer (add a 1) n) ($pochhammer 1 n)))
765
(setq d (pochhammer-quotient (add a 1) 1 x n))
766
(multiple-value-setq (f e)
767
($hypergeo11 (mul -1 n) (add 1 a) x n))
768
(setq e (if e (+ e (* 4 (abs f) double-float-epsilon n)) nil))
769
(orthopoly-return-handler d f e)))
771
`(($gen_laguerre) ,n ,a ,x))))
773
(putprop '$gen_laguerre
775
((unk) "$first" "$gen_laguerre")
776
((unk) "$second" "$gen_laguerre")
779
((mtimes) -1 ((mplus) a n)
780
(($unit_step) n) (($gen_laguerre) ((mplus) -1 n) a x))
781
((mtimes) n (($gen_laguerre) n a x)))
785
(defprop $gen_laguerre tex-gen-laguerre tex)
787
(defun tex-gen-laguerre (x l r)
788
(tex-sub-and-super-scripted-function "L" `(0) nil `(1) t 1 x l r))
790
(setf (get '$gen_laguerre 'dimension) 'dimension-gen-laguerre)
792
(defun dimension-gen-laguerre (form result)
794
(dimension-sub-and-super-scripted-function "L" `(1) `(2) t 3 form)
797
;; See A & S 22.5.16, page 778.
799
(defun $laguerre (n x)
800
(cond ((use-hypergeo n x)
802
(multiple-value-setq (f e) ($hypergeo11 (mul -1 n) 1 x n))
803
(orthopoly-return-handler 1 f e)))
805
`(($laguerre) ,n ,x))))
809
((unk) "$first" "$laguerre")
812
((mtimes) -1 n (($laguerre) ((mplus) -1 n) x))
813
((mtimes) n (($laguerre) n x)))
817
(defprop $laguerre tex-laguerre tex)
819
(defun tex-laguerre (x l r)
820
(tex-sub-and-super-scripted-function "L" `(0) nil nil nil 1 x l r))
822
(setf (get '$laguerre 'dimension) 'dimension-laguerre)
824
(defun dimension-laguerre (form result)
826
(dimension-sub-and-super-scripted-function "L" `(1) nil nil 2 form)
829
(defun $spherical_hankel1 (n x)
831
(cond ((and (integerp n) (< n 0))
832
(setq d (mul '$%i (if (oddp n) 1 -1)))
833
(multiple-value-setq (f e)
834
($spherical_hankel1 (add -1 (mul -1 n)) x))
835
(orthopoly-return-handler d f e))
837
(multiple-value-setq (f e)
838
($hypergeo11 (mul -1 n) (mul -2 n) (mul -2 '$%i x) n))
839
(setq d (mul '$%i (if (= 0 n) 1
840
(simplify `((%genfact) ,(add (mul 2 n) -1)
841
,(add n (rat -1 2)) 2)))
842
(power '$%e (mul '$%i x)) (div -1 (power x (add 1 n)))))
843
(orthopoly-return-handler d f e))
845
`(($spherical_hankel1) ,n ,x)))))
847
(putprop '$spherical_hankel1
849
((unk) "$first" "$spherical_hankel1")
850
((mplus simp) (($spherical_hankel1) ((mplus) -1 n) x)
851
((mtimes simp) -1 ((mplus) 1 n)
852
(($spherical_hankel1) n x) ((mexpt) x -1))))
855
(defprop $spherical_hankel1 tex-spherical-hankel-1 tex)
857
(defun tex-spherical-hankel-1 (x l r)
858
(tex-sub-and-super-scripted-function "h^{(1)}" `(0) nil nil nil 1 x l r))
860
(setf (get '$spherical_hankel1 'dimension) 'dimension-spherical-hankel-1)
862
(defun dimension-spherical-hankel-1 (form result)
863
(let ((form1 `((mexpt) (($\h simp array) ,(nth 1 form))
865
(dimension-function `((,form1 simp) ,(nth 2 form)) result)))
867
;; See A & S 10.1.36.
869
(defun $spherical_hankel2 (n x)
871
(setq x (mul x (power '$%e (mul '$%i '$%pi (add (mul 2 n) 1)))))
873
(setq f ($spherical_hankel1 n x))
874
(if (oddp n) (interval-mult -1 f) f)))
875
(t `(($spherical_hankel2) ,n ,x))))
877
(putprop '$spherical_hankel2
879
((unk) "$first" "$spherical_hankel2")
880
((mplus simp) (($spherical_hankel2) ((mplus) -1 n) x)
881
((mtimes simp) -1 ((mplus) 1 n)
882
(($spherical_hankel2) n x) ((mexpt) x -1))))
885
(defprop $spherical_hankel2 tex-spherical-hankel-2 tex)
887
(defun tex-spherical-hankel-2 (x l r)
888
(tex-sub-and-super-scripted-function "h^{(2)}" `(0) nil nil nil 1 x l r))
890
(setf (get '$spherical_hankel2 'dimension) 'dimension-spherical-hankel-2)
892
(defun dimension-spherical-hankel-2 (form result)
893
(let ((form1 `((mexpt) (($\h simp array) ,(nth 1 form)) (2))))
894
(dimension-function `((,form1 simp) ,(nth 2 form)) result)))
896
;;---------------------------------------------------------------------
897
;; The spherical_bessel functions use the functions p-fun and q-fun.
898
;; See A&S 10.1.8 and 10.1.9 page 437.
906
(setq w (div (mul w `((rat) ,(* -1 (+ n m2 2) (+ n m2 1)
907
(- n m2) (- n (+ m2 1)))
908
,(* 4 (+ m2 1) (+ m2 2)))) x2))
909
(setq s (add s w)))))
912
(let ((s (if (= 0 n) 0 1))
913
(w 1) (m2) (x2 (mul x x))
914
(n1 (floor (/ (- n 1) 2))))
915
(dotimes (m n1 (div (mul n (+ n 1) s) (mul 2 x)))
917
(setq w (div (mul w `((rat) ,(* -1 (+ n m2 3) (+ n m2 2)
918
(- n (+ m2 1)) (- n (+ m2 2)))
919
,(* 4 (+ m2 3) (+ m2 2)))) x2))
920
(setq s (add s w)))))
923
;; See A&S 10.1.8 page 437 and A&S 10.1.15 page 439. When the order
924
;; is an integer and x is a float or a bigfloat, use the slatec code
925
;; for numerical evaluation. Yes, we coerce bigfloats to floats and
928
;; For numerical evaluation, we do our own analytic continuation -- otherwise
929
;; we get factors exp(%i n %pi / 2) that should evaluate to 1,%i,-1,-%, but
930
;; numerically have "fuzz" in their values. The fuzz can cause the spherical
931
;; bessel functions to have nonzero (but small) imaginary values on the
932
;; negative real axis. See A&S 10.1.34
934
(defun $spherical_bessel_j (n x)
935
(cond ((and (eq '$zero (csign ($ratdisrep x)))
936
(or (integerp n) ($featurep n '$integer)))
937
`(($kron_delta) 0 ,n))
939
((and (use-float x) (integerp n))
940
(let ((d 1) (xr) (xi) (z))
941
(setq x ($rectform ($float x)))
942
(setq xr ($realpart x))
943
(setq xi ($imagpart x))
944
(setq z (complex xr xi))
946
(setq d (if (oddp n) -1 1))
949
(setq n (+ 0.5L0 ($float n)))
950
(setq d (* d (sqrt (/ pi (* 2 z)))))
951
(setq d (lisp-float-to-maxima-float d))
952
($expand (mul ($rectform d) ($bessel x n)))))
954
((and (integerp n) (> n -1))
955
(let ((xt (sub x (div (mul n '$%pi) 2))))
957
(mul (p-fun n x) (simplify `((%sin) ,xt)))
958
(mul (q-fun n x) (simplify `((%cos) ,xt)))) x)))
961
(mul (if (oddp n) -1 1) ($spherical_bessel_y (- (+ n 1)) x)))
964
`(($spherical_bessel_j) ,n ,x))))
966
(putprop '$spherical_bessel_j
968
((unk) "$first" "$spherical_bessel_j")
969
((mtimes) ((mexpt) ((mplus) 1 ((mtimes) 2 n)) -1)
971
((mtimes) n (($spherical_bessel_j) ((mplus) -1 n) x))
972
((mtimes) -1 ((mplus) 1 n)
973
(($spherical_bessel_j) ((mplus) 1 n) x)))))
976
(defprop $spherical_bessel_j tex-spherical-bessel-j tex)
978
(defun tex-spherical-bessel-j (x l r)
979
(tex-sub-and-super-scripted-function "j^{(2)}" `(0) nil nil nil 1 x l r))
981
(setf (get '$spherical_bessel_j 'dimension) 'dimension-spherical-bessel-j)
983
(defun dimension-spherical-bessel-j (form result)
984
(let ((form1 `(($\j simp array) ,(nth 1 form))))
985
(dimension-function `((,form1 simp) ,(nth 2 form)) result)))
987
;; For analytic continuation, see A&S 10.1.35.
989
(defun $spherical_bessel_y (n x)
990
(cond ((and (use-float x) (integerp n))
991
(let ((d 1) (xr) (xi) (z))
992
(setq x ($rectform ($float x)))
993
(setq xr ($realpart x))
994
(setq xi ($imagpart x))
995
(setq z (complex xr xi))
997
(setq d (if (oddp n) 1 -1))
1000
(setq n (+ 0.5L0 ($float n)))
1001
(setq d (* d (sqrt (/ pi (* 2 z)))))
1002
(setq d (lisp-float-to-maxima-float d))
1003
($expand (mul ($rectform d) ($bessel_y n x)))))
1005
((and (integerp n) (> n -1))
1006
(let ((xt (add x (div (mul n '$%pi) 2))))
1007
(mul (if (oddp n) 1 -1)
1009
(mul (p-fun n x) (simplify `((%cos) ,xt)))
1010
(mul (q-fun n x) (simplify `((%sin) ,xt)))) x))))
1013
(mul (if (oddp n) 1 -1) ($spherical_bessel_j (- (+ n 1)) x)))
1014
(t `(($spherical_bessel_y) ,n ,x))))
1016
(putprop '$spherical_bessel_y
1018
((unk) "$first" "$spherical_bessel_y")
1019
((mtimes) ((mexpt) ((mplus) 1 ((mtimes) 2 n)) -1)
1021
((mtimes) n (($spherical_bessel_y) ((mplus) -1 n) x))
1022
((mtimes) -1 ((mplus) 1 n)
1023
(($spherical_bessel_y) ((mplus) 1 n) x)))))
1026
(defprop $spherical_bessel_y tex-spherical-bessel-y tex)
1028
(defun tex-spherical-bessel-y (x l r)
1029
(tex-sub-and-super-scripted-function "y^{(2)}" `(0) nil nil nil 1 x l r))
1031
(setf (get '$spherical_bessel_y 'dimension) 'dimension-spherical-bessel-y)
1033
(defun dimension-spherical-bessel-y (form result)
1034
(let ((form1 `(($\y simp array) ,(nth 1 form))))
1035
(dimension-function `((,form1 simp) ,(nth 2 form)) result)))
1037
;; Compute P_n^m(cos(theta)). See Merzbacher, 9.59 page 184
1038
;; and 9.64 page 185, and A & S 22.5.37 page 779. This function
1039
;; lacks error checking; it should only be called by spherical_harmonic.
1041
;; We need to be careful -- for the spherical harmonics we can't use
1042
;; assoc_legendre_p(n,m,cos(x)). If we did, we'd get factors
1043
;; (1 - cos^2(x))^(m/2) that simplify to |sin(x)|^m but we want them
1044
;; to simplify to sin^m(x). Oh my!
1046
(defun assoc-leg-cos (n m x)
1048
(if (= m 0) 1 (mul (simplify
1049
`((%genfact) ,(add (mul 2 m) -1) ,(add m (div 1 2)) 2))
1050
(power (simplify `((%sin) ,x)) m)))
1051
($ultraspherical (sub n m) (add m `((rat) 1 2))
1052
(simplify `((%cos) ,x)))))
1054
(defun $spherical_harmonic (n m th p)
1055
(cond ((and (integerp n) (integerp m) (> n -1))
1056
(cond ((> (abs m) n)
1059
(interval-mult (if (oddp m) -1 1)
1060
($spherical_harmonic n (- m) th (mul -1 p))))
1063
(mul (simplify ($exp (mul '$%i m p)))
1064
(power (div (* (+ (* 2 n) 1) (factorial (- n m)))
1065
(mul '$%pi (* 4 (factorial (+ n m)))))
1067
(assoc-leg-cos n m th)))))
1069
`(($spherical_harmonic) ,n ,m ,th ,p))))
1071
(defprop $spherical_harmonic tex-spherical-harmonic tex)
1073
(defun tex-spherical-harmonic (x l r)
1074
(tex-sub-and-super-scripted-function "Y" `(0) nil `(1) nil 2 x l r))
1076
(setf (get '$spherical_harmonic 'dimension) 'dimension-spherical-harmonic)
1078
(defun dimension-spherical-harmonic (form result)
1080
(dimension-sub-and-super-scripted-function "Y" `(1) `(2) nil 3 form)
1083
(putprop '$spherical_harmonic
1085
((unk) "$first" "$spherical_harmonic")
1086
((unk) "$second" "$spherical_harmonic")
1088
((mtimes) ((rat ) -1 2)
1090
((mtimes) ((mplus) ((mtimes) -1 m) n)
1093
(($spherical_harmonic) n ((mplus) 1 m) theta phi)
1094
((mexpt) $%e ((mtimes) -1 $%i phi)))
1095
((mtimes) ((rat) 1 2)
1097
((mtimes) ((mplus) 1 ((mtimes) -1 m) n)
1100
(($spherical_harmonic) n ((mplus) -1 m) theta phi)
1101
((mexpt) $%e ((mtimes) $%i phi))))
1103
((mtimes) $%i m (($spherical_harmonic) n m theta phi)))
1106
(defun maxima-float-to-lisp-float (y)
1107
(let* ((x ($rectform y))
1110
(cond ((or (floatp xr) (floatp xi))
1111
(setq xr (float xr 1.0d0)
1112
xi (float xi 1.0d0))
1113
(if (= 0.0 xi) xr (complex xr xi)))
1117
(defun lisp-float-to-maxima-float (x)
1119
(add (realpart x) (mul '$%i (imagpart x)))
1122
(defun hypergeo11-float (n b x)
1123
(let ((f0) (fm1) (f) (i 0) (k) (dk) (ak) (bk) (err)
1124
(as (make-array (- 1 n)
1125
:initial-element 0.0d0))
1126
(bs (make-array (- 1 n)
1127
:initial-element 0.0d0))
1128
(fs (make-array (- 1 n)
1129
:initial-element 0.0d0))
1137
(setf (aref fs i) f0)
1139
(setq ak (/ (+ (* 2 i) x) dk))
1140
(setq bk (/ (- i) dk))
1141
(setq f (+ (* ak f0) (* bk fm1)))
1142
(setf (aref as i) ak)
1143
(setf (aref bs i) bk)
1147
(setf (aref fs i) f0)
1154
(setq u (+ (* (aref as k) u0)
1155
(* (aref bs (+ 1 k)) um1)))
1158
(setq err (+ err (abs (* u0 (aref fs k)))))
1160
(values f0 (* 12 double-float-epsilon err))))
1162
(defun hypergeo21-float (n b c x)
1163
(let ((f0) (fm1) (f) (i 0) (k) (dk) (ak) (bk) (err)
1164
(as (make-array (- 1 n)
1165
:initial-element 0.0d0))
1166
(bs (make-array (- 1 n)
1167
:initial-element 0.0d0))
1168
(fs (make-array (- 1 n)
1169
:initial-element 0.0d0))
1176
(setf (aref fs i) f0)
1178
(setq ak (/ (+ (* 2 i) c (- (* x (+ b i)))) dk))
1179
(setq bk (/ (* i (- x 1)) dk))
1180
(setq f (+ (* ak f0) (* bk fm1)))
1181
(setf (aref as i) ak)
1182
(setf (aref bs i) bk)
1186
(setf (aref fs i) f0)
1193
(setq u (+ (* (aref as k) u0)
1194
(* (aref bs (+ 1 k)) um1)))
1198
(setq err (+ err (abs (* (aref fs k) u0)))))
1199
(values f0 (* 12 double-float-epsilon err))))
1201
;; For recursion relations, see A & S 22.7 page 782.
1203
;; jacobi_p(n+1,a,b,x) = (((2*n+a+b+1)*(a^2-b^2) +
1204
;; x*pochhammer(2*n+a+b,3)) * jacobi_p(n,a,b,x) -
1205
;; 2*(n+a)*(n+b)*(2*n+a+b+2)*jacobi_p(n-1,a,b,x))/(2*(n+1)*(n+a+b+1)*(2*n+a+b))
1207
;; ultraspherical(n+1,a,x) = (2*(n+a)*x * ultraspherical(n,a,x) -
1208
;; (n+2*a+1)*ultraspherical(n-1,a,x))/(n+1)
1210
;; chebyshev_t(n+1,x) = 2*x*chebyshev_t(n,x) -chebyshev_t(n-1,x)
1212
;; chebyshev_u(n+1,x) = 2*x*chebyshev_u(n,x) -chebyshev_u(n-1,x)
1214
;; laguerre(n+1,x) = (((2*n+1) - x)*laguerre(n,x) -(n)*laguerre(n-1,x))/(n+1)
1216
;; gen_laguerre(n+1,a,x) = (((2*n+a+1) - x)*gen_laguerre(n,a,x)
1217
;; -(n+a)*gen_laguerre(n-1,a,x))/(n+1)
1219
;; hermite(n+1,x) = 2*x*hermite(n,x) +2*n*hermite(n-1,x)
1221
;; See G & R 8.733.2; A & S 22.7.11 might be wrong -- or maybe I need
1224
;; (2*n+1)*x*assoc_legendre_p(n,m,x) = (n-m+1)*assoc_legendre_p(n+1,m,x)
1225
;; + (n+m)*assoc_legendre_p(n-1,m,x)
1227
;; For the half-integer bessel functions, See A & S 10.1.19
1229
;; fn(n-1,x) + fn(n+1,x) = (2*n+1)*fn(n,x)/x;
1231
(defun check-arg-length (fn n m)
1233
(merror "Function ~:M needs ~:M arguments, instead it received ~:M"
1236
(defun $orthopoly_recur (fn arg)
1237
(if (not ($listp arg))
1238
(merror "The second argument to orthopoly_recur must be a list"))
1239
(cond ((eq fn '$jacobi_p)
1240
(check-arg-length fn 4 (- (length arg) 1))
1241
(let ((n (nth 1 arg))
1246
`((mequal) (($jacobi_p ) ((mplus) 1 ,n) ,a ,b ,x)
1247
((mtimes) ((rat) 1 2) ((mexpt) ((mplus) 1 ,n) -1)
1248
((mexpt) ((mplus) 1 ,a ,b ,n) -1)
1249
((mexpt) ((mplus) ,a ,b ((mtimes) 2 ,n)) -1)
1251
((mtimes) -2 ((mplus) ,a ,n) ((mplus) ,b ,n)
1252
((mplus) 2 ,a ,b ((mtimes) 2 ,n))
1253
(($jacobi_p ) ((mplus) -1 ,n) ,a ,b ,x))
1254
((mtimes) (($jacobi_p ) ,n ,a ,b ,x)
1257
((mplus) ((mexpt) ,a 2)
1258
((mtimes) -1 ((mexpt) ,b 2)))
1259
((mplus) 1 ,a ,b ((mtimes) 2 ,n)))
1260
((mtimes) ((mplus) ,a ,b ((mtimes) 2 ,n))
1261
((mplus) 1 ,a ,b ((mtimes) 2 ,n))
1262
((mplus) 2 ,a ,b ((mtimes) 2 ,n)) ,x)))))))))
1264
((eq fn '$ultraspherical)
1265
(check-arg-length fn 3 (- (length arg) 1))
1266
(let ((n (nth 1 arg))
1270
`((mequal) (($ultraspherical) ((mplus) 1 ,n) ,a ,x)
1271
((mtimes) ((mexpt) ((mplus) 1 ,n) -1)
1273
((mtimes) -1 ((mplus) -1 ((mtimes) 2 ,a) ,n)
1274
(($ultraspherical) ((mplus) -1 ,n) ,a ,x))
1275
((mtimes) 2 ((mplus) ,a ,n)
1276
(($ultraspherical) ,n ,a ,x) ,x)))))))
1278
((member fn `($chebyshev_t $chebyshev_u) :test 'eq)
1279
(check-arg-length fn 2 (- (length arg) 1))
1280
(let ((n (nth 1 arg))
1283
`((mequal ) ((,fn) ((mplus ) 1 ,n) ,x)
1285
((mtimes ) -1 ((,fn) ((mplus ) -1 ,n) ,x))
1286
((mtimes ) 2 ((,fn) ,n ,x) ,x))))))
1288
((member fn '($legendre_p $legendre_q) :test 'eq)
1289
(check-arg-length fn 2 (- (length arg) 1))
1290
(let* ((n (nth 1 arg))
1292
(z (if (eq fn '$legendre_q)
1293
`((mtimes) -1 (($kron_delta) ,n 1)) 0)))
1295
`((mequal) ((,fn) ,n ,x)
1297
((mtimes) ((mexpt) ,n -1)
1299
((mtimes) ((mplus) 1 ((mtimes) -1 ,n))
1300
((,fn) ((mplus) -2 ,n) ,x))
1301
((mtimes) ((mplus) -1 ((mtimes) 2 ,n))
1302
((,fn) ((mplus) -1 ,n) ,x) ,x)))
1305
((member fn '($assoc_legendre_p $assoc_legendre_q) :test 'eq)
1306
(check-arg-length fn 3 (- (length arg) 1))
1307
(let ((n (nth 1 arg))
1311
`((mequal) ((,fn) ((mplus) 1 ,n) ,m ,x)
1313
((mexpt) ((mplus) 1 ((mtimes) -1 ,m) ,n) -1)
1316
((mplus) ((mtimes) -1 ,m)
1318
((,fn) ((mplus) -1 ,n) ,m ,x))
1319
((mtimes) ((mplus) 1 ((mtimes) 2 ,n))
1320
((,fn) ,n ,m ,x) ,x)))))))
1323
(check-arg-length fn 2 (- (length arg) 1))
1324
(let ((n (nth 1 arg))
1327
`((mequal ) (($laguerre ) ((mplus ) 1 ,n) ,x)
1328
((mtimes ) ((mexpt ) ((mplus ) 1 ,n) -1)
1330
((mtimes ) -1 ,n (($laguerre ) ((mplus ) -1 ,n) ,x))
1331
((mtimes ) (($laguerre ) ,n ,x)
1332
((mplus ) 1 ((mtimes ) 2 ,n) ((mtimes ) -1 ,x)))))))))
1334
((eq fn '$gen_laguerre)
1335
(check-arg-length fn 3 (- (length arg) 1))
1336
(let ((n (nth 1 arg))
1340
`((mequal) (($gen_laguerre) ((mplus) 1 ,n) ,a ,x)
1341
((mtimes) ((mexpt ) ((mplus) 1 ,n) -1)
1343
((mtimes) -1 ((mplus) ,a ,n)
1344
(($gen_laguerre) ((mplus) -1 ,n) ,a ,x))
1345
((mtimes) (($gen_laguerre) ,n ,a ,x)
1346
((mplus) 1 ,a ((mtimes) 2 ,n) ((mtimes ) -1 ,x)))))))))
1349
(check-arg-length fn 2 (- (length arg) 1))
1350
(let ((n (nth 1 arg))
1353
`((mequal) (($hermite) ((mplus) 1 ,n) ,x)
1355
((mtimes) -2 ,n (($hermite) ((mplus) -1 ,n) ,x))
1356
((mtimes) 2 (($hermite) ,n ,x) ,x))))))
1358
((member fn `($spherical_bessel_j $spherical_bessel_y
1359
$spherical_hankel1 $spherical_hankel2)
1361
(check-arg-length fn 2 (- (length arg) 1))
1362
(let ((n (nth 1 arg))
1365
`((mequal) ((,fn) ((mplus) 1 ,n) ,x)
1366
((mplus) ((mtimes) -1 ((,fn ) ((mplus) -1 ,n) ,x))
1367
((mtimes) ((,fn) ,n ,x) ((mexpt) ,x -1))
1368
((mtimes) 2 ,n ((,fn ) ,n ,x) ((mexpt) ,x -1)))))))
1370
(t (merror "A recursion relation for ~:M isn't known to Maxima" fn))))
1372
;; See A & S Table 2.2, page 774.
1374
(defun $orthopoly_weight (fn arg)
1375
(if (not ($listp arg))
1376
(merror "The second argument to orthopoly_weight must be a list"))
1378
(if (not (or ($symbolp (car (last arg))) ($subvarp (car (last arg)))))
1379
(merror "The last element of the second argument to orthopoly_weight must
1380
be a symbol or a subscripted symbol, instead found ~:M" (car (last arg))))
1382
(if (not (every #'(lambda (s)
1383
($freeof (car (last arg)) s)) (butlast (cdr arg))))
1384
(merror "Only the last element of ~:M may depend on the integration
1385
variable ~:M" arg (car (last arg))))
1387
(cond ((eq fn '$jacobi_p)
1388
(check-arg-length fn 4 (- (length arg) 1))
1389
(let ((a (nth 2 arg))
1394
((mtimes) ((mexpt) ((mplus) 1 ((mtimes) -1 ,x)) ,a)
1395
((mexpt) ((mplus ) 1 ,x) ,b))
1398
((eq fn '$ultraspherical)
1399
(check-arg-length fn 3 (- (length arg) 1))
1400
(let ((a (nth 2 arg))
1404
((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) ,x 2)))
1405
((mplus) ((rat) -1 2) ,a)) -1 1))))
1407
((eq fn '$chebyshev_t)
1408
(check-arg-length fn 2 (- (length arg) 1))
1409
(let ((x (nth 2 arg)))
1412
((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) ,x 2)))
1413
((rat) -1 2)) -1 1))))
1415
((eq fn '$chebyshev_u)
1416
(check-arg-length fn 2 (- (length arg) 1))
1417
(let ((x (nth 2 arg)))
1420
((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) ,x 2)))
1421
((rat) 1 2)) -1 1))))
1423
((eq fn '$legendre_p)
1424
(check-arg-length fn 2 (- (length arg) 1))
1428
(check-arg-length fn 2 (- (length arg) 1))
1429
(let ((x (nth 2 arg)))
1431
`((mlist) ((mexpt) $%e ((mtimes) -1 ,x)) 0 $inf))))
1433
((eq fn '$gen_laguerre)
1434
(check-arg-length fn 3 (- (length arg) 1))
1435
(let ((a (nth 2 arg))
1439
((mtimes) ((mexpt) ,x ,a)
1440
((mexpt) $%e ((mtimes) -1 ,x))) 0 $inf))))
1443
(check-arg-length fn 2 (- (length arg) 1))
1444
(let ((x (nth 2 arg)))
1446
`((mlist) ((mexpt) $%e ((mtimes) -1 ((mexpt) ,x 2)))
1447
((mtimes ) -1 $inf) $inf))))
1449
(t (merror "a weight for ~:m isn't known to maxima" fn))))
added
removed
share/orthopoly/orthopoly.lisp
