fixes typos and a missing reference.

[maxima.git] / share / orthopoly / orthopoly.lisp

;; 8/8 | 5/5 
;; Copyright (C) 2000, 2001, 2003, 2008, 2009 Barton Willis

#|
  This is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License,
  http://www.gnu.org/copyleft/gpl.html.

 This software has NO WARRANTY, not even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Maxima code for evaluating orthogonal polynomials listed in Chapter 22 of Abramowitz and Stegun (A & S). 
|#


(in-package :maxima)
($put '$orthopoly 1.0 '$version)

;; If the input can be evaluated to a floating point number (either
;; real or complex), convert the input to a Common Lisp complex number.
;; When the input doesn't evaluate to a float, return the input.

(defun maxima-to-lisp-complex-float (a)
  (let* (($ratprint nil)
         (b ($rectform a))
         (br ($float ($realpart b)))
         (bi ($float ($imagpart b))))
    (cond ((and (numberp br) (numberp bi))
           (if (= bi 0) (float br)
             (complex (float br) 
                      (float bi))))
          (t a))))

;; Return true iff a is a float or a complex number with either a
;; real or imaginary part that is a float.

(defun xcomplexp (a)
  (or (floatp a)
      (and (complexp a) (or (floatp (realpart a)) (floatp (imagpart a))))))

;; Convert a Lisp complex number into a Maxima complex number. When the
;; input isn't a Lisp complex number, return the argument.

(defun lisp-complex-to-maxima (x)
  (if (complexp x) (add (realpart x) (mul '$%i (imagpart x))) x))
   
;; When orthopoly_returns_intervals is true, floating point evaluation 
;; returns an interval using the form (($interval) c r), where c is the 
;; center of the interval and r is its radius.  We don't provide the user
;; with any tools for working with intervals; if a user wants

(defmvar $orthopoly_returns_intervals t)

(defun orthopoly-return-handler (d f e)
  (cond ((or (floatp f) (complexp f))
         (let ((df (maxima-to-lisp-complex-float d))
               (ef))
           (setq f (lisp-complex-to-maxima
                    (if (numberp df) (* df f) (mul d f))))
           (setq ef (if (and (numberp df) (numberp e)) (abs (* df e)) nil))
           (cond ((and $orthopoly_returns_intervals (floatp ef))
                  `(($interval simp) ,f ,ef))
                 (t f))))
        (t (mul d f))))
                             
;; When a user requests the derivative of an a function in this package
;; with respect to the order or some other parameter, return a form 
;; ((unk) input from user). We "simplify" this form by printing an error.

(defprop unk simp-unk operators)
(defun simp-unk (x y z)
  (declare (ignore y z))
  (merror "Maxima doesn't know the derivative of ~:M with respect the ~:M argument" (nth 2 x) (nth 1 x)))

;; A left continuous unit step function; thus 
;;
;;       unit_step(x) = 0 for x <= 0 and 1 for x > 0.  
;;
;; This function differs from (1 + signum(x))/2 which isn't left or right
;; continuous at 0.

(defprop $unit_step "\\Theta" texword)

(defun simp-unit-step (a y z)
  (oneargcheck a)
  (setq y (simpcheck (cadr a) z))
  (let ((s (csign y)))
    (cond ((or (eq s '$nz) (eq s '$zero) (eq s '$neg)) 0)
          ((eq s '$pos) 1)
          (t `(($unit_step simp) ,y)))))

;; We barely support intervals; these functions aren't for user-level use.
;; The function intervalp returns true iff its argument is an interval.

(defun $intervalp (a)
  (and (consp a) (consp (car a)) (eq (caar a) '$interval)))

;; Multiply x by a, where x isn't an interval and a might be an 
;; interval.  When a isn't an interval or x is an interval, return
;; x * a.  The optional argument dx is an upper bound for the 
;; relative error in x.

(defun interval-mult (x a &optional dx)
  (if (null dx) (setq dx 0))
  (cond ((and ($intervalp a) (not ($intervalp x)))
         (let ((p (nth 1 a))
               (q (nth 2 a)))
           (if (or (floatp p) (floatp q)) 
               (setq x ($float ($rectform x))))
           (setq p ($expand (mult x p)))
           (setq q ($expand (mult x (add q (simplify `((mabs) ,(mul p dx)))))))
           (setq q (simplify `((mabs) ,q)))
           `(($interval) ,p ,q)))
        (t (mult x a))))
           
;; TeX a function with subscripts and superscripts.  The string fn is the
;; function name, the list sub holds the positions of the subscripts, the list
;; sup holds the positions of the superscripts, and i is the position of the 
;; function argument.  When b1 is true, the subscript is surrounded by parens;
;; when b2 is true, the superscript is surrounded by parens.  The lists sub and
;; sup may be nil, but the function must have at least on argument.

(defun tex-sub-and-super-scripted-function (fn sub b1 sup b2 i x l r)
  (setq x (cdr x))
  (let ((lo) (hi) (s1) (s2))
    (setq s1 (if b1 `("\\left(")  nil))
    (setq s2 (if b1 `("\\right)") nil))
    (dolist (i sub)
      (setq lo (cons (nth i x) lo)))
    (setq lo (if lo (tex-list (nreverse lo) s1 s2 ",") nil))
    (setq s1 (if b2 `("\\left(")  nil))
    (setq s2 (if b2 `("\\right)") nil))
    (dolist (i sup)
      (setq hi (cons (nth i x) hi)))
    (setq hi (if hi (tex-list (nreverse hi) s1 s2 ",") nil))
    (append l `(,fn)
            (if lo `("_{",@lo "}") nil)
            (if hi `("^{" ,@hi "}") nil)           
            `(,@(tex-list (nthcdr i x) `("\\left(") `("\\right)") ","))
            r)))
            
(defun dimension-sub-and-super-scripted-function (fn sub sup b2 k x)
  (let ((lo) (hi) (form))
    (dolist (i sub)
      (setq lo (cons (nth i x) lo)))
    (setq lo (nreverse lo))
    (dolist (i sup)
      (setq hi (cons (nth i x) hi)))
    (setq hi (nreverse hi))
    (cond ((null hi)
           (setq form `((,fn simp array) ,@lo)))
          (b2
           (setq form `((mexpt) ((,fn simp array) ,@lo) (("") ,@hi))))
          (t
           (setq form `((mexpt) ((,fn simp array) ,@lo) ,@hi))))
    `((,form simp) ,@(nthcdr k x))))
 
;; Return true iff 
;;   (1) each member of the list a evaluates to a floating
;;       point number (using $float) and 
;;   (2) at least one member of a has a real or imaginary part that is a 
;;       floating point number *or* a bigfloat.
;; When we find an member of a that doesn't evaluate to a float, 
;; immediately bail out and return nil.

(defun use-float (&rest a)
  (let ((xr) (xi) (float-found nil) (okay t))
    (dolist (x a)
      (setq x ($rectform x))
      (setq xr ($realpart x))
      (setq xi ($imagpart x))
      (setq float-found (or float-found (floatp xr) (floatp xi)
                            ($bfloatp xr) ($bfloatp xi)))
      (setq okay (and (floatp ($float xr)) (floatp ($float xi))))
      (if (not okay) (return)))
    (and okay float-found)))
         
;; When n is a nonnegative integer, return 1F1(a,b,x); that is return
;; sum(pochhammer(a,k) * x^k /pochhammer(b,k),k,0,n). Regardless of the
;; value of n,  when x = 0, immediately return 1.

;; Any orthopoly function that calls hypergeo should check that n is an
;; integer with n > -1; thus the summation code shouldn't get called for
;; any orthopoly function.  It wouldn't be terrible if it happened -- I
;; think the summation form isn't useful and a user can get into trouble
;; with things like subst(0,x,sum(x^k,k,0,n)).

(defun $hypergeo11 (a b x n)
  (cond ((like x 0) 1)
        ((and (integerp n) (> n -1))
         (cond ((and (like a (- n)) (use-float b x))
                (let ((f) (e))
                  (multiple-value-setq (f e)
                    (hypergeo11-float (- n) (maxima-to-lisp-complex-float b)
                                      (maxima-to-lisp-complex-float x)))
                  (values f e)))
               (t
                (let ((sum 1) (cf 1))
                  (dotimes (i n sum)
                    (setq cf (div (mul cf (add i a) x) (mul (add i b) (+ 1 i))))
                    (setq sum (add cf sum)))))))
        (t
; The following is replaced with simplifying code.
;        `((%sum )
;          ((mtimes) ((mexpt) ((mfactorial) ,$genindex) -1)
;           (($pochhammer) ,a ,$genindex)
;           ((mexpt) (($pochhammer ) ,b ,$genindex) -1) 
;           ((mexpt) ,x ,$genindex)) ,$genindex 0 ,n))))
         (let ((index (gensumindex)))
           (simplify
             (list '(%sum)
                   (mul (inv (take '(mfactorial) index))
                        (take '($pochhammer) a index)
                        (inv (take '($pochhammer) b index))
                        (power x index))
                   index 0 n))))))

;; return the F[2,1] hypergeometic sum from 0 to n. Use genindex as the 
;; sum index; genindex is defined in asum.lisp

(defun $hypergeo21 (a b c x n)
  (cond ((like x 0) 1)
        ((and (integerp n) (> n -1))
         (cond ((and (like a (- n)) (use-float b c x))
                (let ((f) (e))
                  (multiple-value-setq (f e)
                    (hypergeo21-float (- n) 
                                      (maxima-to-lisp-complex-float b)
                                      (maxima-to-lisp-complex-float c)
                                      (maxima-to-lisp-complex-float x)))
                  (values f e)))
               (t
                (let ((sum 1) (cf 1))
                  (dotimes (i n sum)
                    (setq cf (div (mul cf (add i a) (add i b) x) (mul (+ 1 i) 
                                                                      (add i c))))
                    (setq sum (add cf sum)))))))
        
        (t
; The following is replaced with simplifying code.
;        `((%sum)
;          ((mtimes) (($pochhammer) ,a ,$genindex) (($pochhammer) ,b ,$genindex)
;           ((mexpt) (($pochhammer) ,c ,$genindex) -1)
;           ((mexpt) ((mfactorial) ,$genindex) -1)
;           ((mexpt) ,x ,$genindex)) 
;          ,$genindex 0 ,n))))
         (let ((index (gensumindex)))
           (simplify
             (list '(%sum)
                   (mul (take '($pochhammer) a index)
                        (take '($pochhammer) b index)
                        (inv (take '($pochhammer) c index))
                        (inv (take '(mfactorial) index))
                        (power x index))
                   index 0 n))))))

(eval-when
    #+gcl (load eval)
    #-gcl (:load-toplevel :execute)
    (let (($context '$global) (context '$global))
      (meval '(($declare) $pochhammer $complex))))

(defmvar $pochhammer_max_index 100)

;; This disallows noninteger assignments to $pochhammer_max_index.

(setf (get '$pochhammer_max_index 'assign)
      #'(lambda (a b) 
          (declare (ignore a))
          (if (not (and (atom b) (integerp b)))
              (progn
                (mtell "The value of `pochhammer_max_index' must be an integer.~%")
                'munbindp))))

(defun $pochhammer (x n)
  (take '($pochhammer) x n))

(in-package #-gcl #:bigfloat #+gcl "BIGFLOAT")

;; Numerical evaluation of pochhammer using the bigfloat package.
(defun pochhammer (x n)
  (let ((acc 1))
    (if (< n 0) (/ 1 (pochhammer (+ x n) (- n)))
      (dotimes (k n acc)
        (setq acc (* acc (+ k x)))))))

(in-package :maxima)

(defun simp-pochhammer (e y z)
  (declare (ignore y))
  (let ((x) (n) (return-a-rat))
    (twoargcheck e)
    (setq return-a-rat ($ratp (second e)))
    (setq x (simplifya (specrepcheck (second e)) z))
    (setq n (simplifya (specrepcheck (third e)) z))
 
    (cond ((or ($taylorp (second e)) ($taylorp (third e)))
           (setq x (simplifya (second e) z))
           (setq n (simplifya (third e) z))
           ($taylor (div (take '(%gamma) (add x n)) (take '(%gamma) x))))
           
          ((eql n 0) 1)
          
          ;; pochhammer(1,n) = n! (factorial is faster than bigfloat::pochhammer.)
          ((eql x 1) (take '(mfactorial) n))
     
          ;; pure numeric evaluation--use numeric package.
          ((and (integerp n) (complex-number-p x '$numberp))
           (maxima::to (bigfloat::pochhammer (bigfloat::to x) (bigfloat::to n))))

          ((and (integerp (mul 2 n)) (integerp (mul 2 x)) (> (mul 2 n) 0) (> (mul 2 x) 0))
           (div (take '(%gamma) (add x n)) (take '(%gamma) x)))

          ;; Use a reflection identity when (great (neg n) n) is true. Specifically,
          ;; use pochhammer(x,-n) * pochhammer(x-n,n) = 1; thus pochhammer(x,n) = 1/pochhammer(x+n,-n).
          ((great (neg n) n)
           (div 1 (take '($pochhammer) (add x n) (neg n))))
          
          ;; Expand when n is an integer and either abs(n) < $expop or abs(n) < $pochhammer_max_index.
          ;; Let's give $expand a bit of early help.
          ((and (integerp n) (or (< (abs n) $expop) (< (abs n) $pochhammer_max_index)))
           (if (< n 0) (div 1 (take '($pochhammer) (add x n) (neg n)))
             (let ((acc 1))
               (if (or (< (abs n) $expop) return-a-rat) (setq acc ($rat acc) x ($rat x)))
               (dotimes (k n (if return-a-rat acc ($ratdisrep acc)))
                 (setq acc (mul acc (add k x)))))))
           
          ;; return a nounform.
          (t `(($pochhammer simp) ,x ,n)))))

(putprop '$pochhammer
         '((x n)
           ((mtimes) (($pochhammer) x n)
            ((mplus) ((mtimes) -1 ((mqapply) (($psi array) 0) x))
             ((mqapply) (($psi array) 0) ((mplus) n x)))) 
           ((mtimes) (($pochhammer) x n)
            ((mqapply) (($psi array) 0) ((mplus) n x))))
         'grad)

(defprop $pochhammer tex-pochhammer tex)

(defun tex-pochhammer (x l r)
  (setq x (mapcar #'(lambda (s) (tex s nil nil nil nil)) (cdr x)))
  (append l 
          `("\\left(")
          (nth 0 x)
          `("\\right)_{")
          (nth 1 x)
          `("}")
          r))

(setf (get '$pochhammer 'dimension) 'dimension-pochhammer)

(defun dimension-pochhammer (form result)
  (setq form `(( (("") ,(nth 1 form)) simp array) ,(nth 2 form)))
  (dimension-array form result))

;; pochhammer-quotient(a b x n) returns pochhammer( a,n) / pochhammer(b,n).  
;; Only when one of a, b, or x  is a floating point number and the others are 
;; numeric types (that is when (use-float a b x) evaluates to true) does this 
;; function differ from using (div ($pochhammer a n) ($pochhammer b n)).  
;; In the floating point case, we arrange the operations to reduce rounding 
;; errors and to reduce over and underflows.

(defun pochhammer-quotient (a b x n)
  (cond ((mminusp n)
         (pochhammer-quotient b a x (neg n)))
        ((and (integerp n) (use-float a b x))
         (let ((acc 1.0))
           (setq a (maxima-to-lisp-complex-float a))
           (setq b (maxima-to-lisp-complex-float b))
           (dotimes (i n (lisp-float-to-maxima-float acc))
             (setq acc (* acc (/ (+ i a) (+ i b)))))))
        (t (div ($pochhammer a n) ($pochhammer b n)))))
  
(defun use-hypergeo (n x)
  (declare (ignore x))
  (or (and (integerp n) (> n -1))
      ($featurep n '$integer)))
;      (and ($featurep n '$integer) ($ratnump x))))

;; See A&S 22.5.42, page 779. For integer n, use the identity
;;     binomial(n+a,n) = pochhammer(a+1,n)/pochhammer(1,n)

;; The condition number of the pochhammer function is 
;;     |1 + x/(x+1) + x/(x+2) + ... + x/(x+n-1)| <= n.
;; The relative floating point error in computing the pochhammer 
;; function is bounded by 3 n eps.  Putting these errors together,
;; the error in d is bounded 4 n |d| eps.

;; We're looking at (d +|- 4 n eps |d|)(f +|- e) = d (f +|- (e + 4 n eps |f|)) + 
;; O(eps^2). 

(defun $jacobi_p (n a b x)
  (cond ((use-hypergeo n x)
         (let ((f) (d) (e))
           ;(setq d (div ($pochhammer (add a 1) n) ($pochhammer 1 n)))
           (setq d (pochhammer-quotient (add a 1) 1 x n))
           (multiple-value-setq (f e)
             ($hypergeo21 (mul -1 n) (add n a b 1) (add a 1)
                          (div (add 1 (mul -1 x)) 2) n))
           (setq e (if e (+ e (* 4 n (abs f) flonum-epsilon)) nil))
           (orthopoly-return-handler d f e)))
        (t `(($jacobi_p simp) ,n ,a ,b ,x))))

(putprop '$jacobi_p
         '((n a b x)
           ((unk) "$first" "$jacobi_p")
           ((unk) "$second" "$jacobi_p")
           ((unk) "$third" "$jacobi_p")

           ((mtimes)
            ((mexpt) ((mplus ) a b ((mtimes ) 2 n)) -1)
            ((mplus)
             ((mtimes) 2
              (($unit_step) n)
              ((mplus) a n) ((mplus) b n)
              (($jacobi_p) ((mplus) -1 n) a b x))
             ((mtimes) n (($jacobi_p) n a b x)
              ((mplus) a ((mtimes ) -1 b)
               ((mtimes)
                ((mplus) ((mtimes) -1 a) ((mtimes ) -1 b)
                 ((mtimes) -2 n)) x))))
            ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt ) x 2))) -1)))
         'grad)
           
(defprop $jacobi_p tex-jacobi-poly tex)

(defun tex-jacobi-poly (x l r)
  (tex-sub-and-super-scripted-function "P" `(0) nil `(1 2) t 3 x l r))

(setf (get '$jacobi_p 'dimension) 'dimension-jacobi-p)

(defun dimension-jacobi-p (form result)
  (dimension-function
   (dimension-sub-and-super-scripted-function "P" `(1) `(2 3) t 4 form)
   result))
          
;; See A&S 22.5.46, page 779.

(defun $ultraspherical (n a x)
  (cond ((use-hypergeo n x)
         (let ((f) (d) (e))
           ;(setq d (div ($pochhammer (mul 2 a) n) ($pochhammer 1 n)))
           (setq d (pochhammer-quotient (mul 2 a) 1 x n))
           (multiple-value-setq (f e)
             ($hypergeo21 (mul -1 n) (add n (mul 2 a)) (add a (div 1 2))
                          (div (add 1 (mul -1 x)) 2) n))
           (setq e (if e (+ e (* 4 n (abs f) flonum-epsilon)) nil))
           (orthopoly-return-handler d f e)))
        (t `(($ultraspherical simp) ,n ,a ,x))))

(putprop '$ultraspherical 
         '((n a x)
           ((unk) "$first" "$ultraspherical")
           ((unk) "$second" "$ultrapsherical")
           ((mtimes)
            ((mplus)
             ((mtimes)
              (($unit_step) n)
              ((mplus) -1 ((mtimes) 2 a) n)
              (($ultraspherical) ((mplus) -1 n) a x))
             ((mtimes) -1 n (($ultraspherical) n a x) x))
            ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))) -1)))
         'grad)            

(defprop $ultraspherical tex-ultraspherical tex)

(defun tex-ultraspherical (x l r)
  (tex-sub-and-super-scripted-function "C" `(0) nil `(1) t 2 x l r))

(setf (get '$ultraspherical 'dimension) 'dimension-ultraspherical)

(defun dimension-ultraspherical (form result)
    (dimension-function
     (dimension-sub-and-super-scripted-function "C" `(1) `(2) t 3 form)
     result))
  
(defun $chebyshev_t (n x)
  (cond ((use-hypergeo n x)
         (let ((f) (e))
           (multiple-value-setq (f e)
             ($hypergeo21 (mul -1 n) n (rat 1 2) (div (add 1 (mul -1 x)) 2) n))
           (orthopoly-return-handler 1 f e)))
        (t `(($chebyshev_t simp) ,n ,x))))

(putprop '$chebyshev_t 
         '((n x)
           ((unk) "$first" "$chebyshev_t")
           ((mtimes)
            ((mplus)
             ((mtimes) n (($chebyshev_t) ((mplus ) -1 n) x))
             ((mtimes ) -1 n (($chebyshev_t) n x) x))
            ((mexpt) ((mplus ) 1 ((mtimes) -1 ((mexpt) x 2))) -1)))
           'grad)

(defprop $chebyshev_t tex-chebyshev-t tex)

(defun tex-chebyshev-t (x l r)
  (tex-sub-and-super-scripted-function "T" `(0) nil nil nil 1 x l r))

(setf (get '$chebyshev_t 'dimension) 'dimension-chebyshev-t)

(defun dimension-chebyshev-t (form result)
  (dimension-function
   (dimension-sub-and-super-scripted-function "T" `(1) nil nil 2 form)
   result))

;; See A & S 22.5.48, page 779.

(defun $chebyshev_u (n x)
  (cond ((use-hypergeo n x)
         (let ((f) (d) (e))
           (setq d (add 1 n)) 
           (multiple-value-setq (f e)
             ($hypergeo21 (mul -1 n) (add 2 n) (rat 3 2)
                          (div (add 1 (mul -1 x)) 2) n))
           (orthopoly-return-handler d f e)))
        (t `(($chebyshev_u simp) ,n ,x))))

(putprop '$chebyshev_u
         '((n x)
           ((unk) "$first" "$chebyshev_u")
           ((mtimes)
            ((mplus)
             ((mtimes)
              (($unit_step) n)
              ((mplus) 1 n) (($chebyshev_u) ((mplus) -1 n) x))
             ((mtimes) -1 n (($chebyshev_u) n x) x))
            ((mexpt) ((mplus ) 1 ((mtimes) -1 ((mexpt) x 2))) -1)))
         'grad) 

(defprop $chebyshev_u tex-chebyshev-u tex)

(defun tex-chebyshev-u (x l r)
  (tex-sub-and-super-scripted-function "U" `(0) nil nil nil 1 x l r))

(setf (get '$chebyshev_u 'dimension) 'dimension-chebyshev-u)

(defun dimension-chebyshev-u (form result)
  (dimension-function
   (dimension-sub-and-super-scripted-function "U" `(1) nil nil 2 form)
   result))

;; See A&S 8.2.1 page 333 and 22.5.35 page 779.  We evaluate the legendre
;; polynomials as jacobi_p(n,0,0,x).  Eat less exercise more.

(defun $legendre_p (n x)
  (cond ((use-hypergeo n x) 
         (if (and (integerp n) (< n 0))
           ($legendre_p (- (abs n) 1) x)
           ($jacobi_p n 0 0 x)))
        (t `(($legendre_p simp) ,n ,x))))

(putprop '$legendre_p 
         '((n x) 
           ((unk) "$first" "$legendre_p")
           ((mtimes)
             ((mplus)
              ((mtimes) n (($legendre_p) ((mplus) -1 n) x))
              ((mtimes) -1 n (($legendre_p) n x) x))
             ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))) -1)))
         'grad)

(defprop $legendre_p tex-legendre-p tex)

(defun tex-legendre-p (x l r)
  (tex-sub-and-super-scripted-function "P" `(0) nil nil nil 1 x l r))

(setf (get '$legendre_p 'dimension) 'dimension-legendre-p)

(defun dimension-legendre-p (form result)
  (dimension-function
   (dimension-sub-and-super-scripted-function "P" `(1) nil nil 2 form)
   result))
  
(defun $legendre_q (n x)
  (if (and (integerp n) (> n -1)) 
      ($assoc_legendre_q n 0 x)
    `(($legendre_q simp) ,n ,x)))

(putprop '$legendre_q 
         '((n x) 
           ((unk) "$first" "$legendre_p")
           ((mplus)
            ((mtimes) -1 ((%kron_delta) 0 n)
             ((mexpt) ((mplus) -1 ((mexpt) x 2)) -1)) 
            ((mtimes)
             ((mplus)
              ((mtimes) n (($legendre_q) ((mplus) -1 n) x))
              ((mtimes) -1 n (($legendre_q) n x) x))
             ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))) -1))))
         'grad)

(defprop $legendre_q tex-legendre-q tex)

(defun tex-legendre-q (x l r)
  (tex-sub-and-super-scripted-function "Q" `(0) nil nil nil 1 x l r))

(setf (get '$legendre_q 'dimension) 'dimension-legendre-q)

(defun dimension-legendre-q (form result)
 (dimension-function
   (dimension-sub-and-super-scripted-function "Q" `(1) nil nil 2 form)
   result))
         
;; See A & S 8.6.7 and 8.2.6 pages 333 and 334. I chose the 
;; definition that is real valued on (-1,1).  

;; For negative m, A & S 8.2.6 page 333 and  G & R 8.706 page 1000
;; disagree; the factor of exp(i m pi) in A & S 8.1.6 suggests to me that
;; A & S  8.2.6 is bogus.  As further evidence, Macsyma 2.4 seems to 
;; agree with G & R 8.706. I'll use G & R.

;; Return assoc_legendre(0,m,x). This isn't a user-level function; we 
;; don't check that m is a positive integer.

(defun q0m (m x)
  (cond ((< m 0)
         (merror "function q0m called with negative order. File a bug report"))
        ((= m 0)
         (div (simplify `((%log) ,(div (add 1 x) (sub 1 x)))) 2))
        (t
         (mul (factorial (- m 1)) `((rat simp) 1 2)
              (if (oddp m) -1 1) (power (sub 1 (mult x x)) (div m 2))
              (add 
               (mul (if (oddp m) 1 -1) (power (add 1 x) (neg m)))
               (power (sub 1 x) (neg m)))))))
  
;; Return assoc_legendre(1,m,x).  This isn't a user-level function; we 
;; don't check that m is a positive integer; we don't check that m is 
;; a positive integer.

(defun q1m (m x)
  (cond ((< m 0)
         (merror "function q1m called with negative order. File a bug report"))
        ((= m 0)
         (sub (mul x (q0m 0 x)) 1))
        ((= m 1)
         (mul -1 (power (sub 1 (mult x x)) `((rat simp) 1 2))
              (sub (q0m 0 x) (div x (sub (mul x x) 1)))))
        (t
         (mul (if (oddp m) -1 1) (power (sub 1 (mult x x)) (div m 2))
              (add
               (mul (factorial (- m 2)) `((rat simp) 1 2) (if (oddp m) -1 1)
                    (sub (power (add x 1) (sub 1 m))
                         (power (sub x 1) (sub 1 m))))
               
               (mul (factorial (- m 1)) `((rat simp) 1 2) (if (oddp m) -1 1)
                    (add (power (add x 1) (neg m)) 
                         (power (sub x 1) (neg m)))))))))

;; Return assoc_legendre_q(n,n,x). I don't have a reference that gives
;; a formula for assoc_legendre_q(n,n,x). To figure one out,  I used
;; A&S 8.2.1 and a formula for assoc_legendre_p(n,n,x).  

;; After finishing the while loop, q = int(1/(1-t^2)^(n+1),t,0,x). 

(defun assoc-legendre-q-nn (n x)
  (let ((q) 
        (z (sub 1 (mul x x))) 
        (i 1))
    (setq q (div (simplify `((%log) ,(div (add 1 x) (sub 1 x)))) 2))
    (while (<= i n)
      (setq q (add (mul (sub 1 `((rat) 1 ,(* 2 i))) q)
                   (div x (mul 2 i (power z i)))))
      (incf i))
    (mul (expt -2 n) (factorial n) (power z (div n 2)) q)))

;; Use degree recursion to find the assoc_legendre_q function. 
;; See A&S 8.5.3. When i = m in the while loop, we have a special
;; case.  For negative order, see A&S 8.2.6.

(defun $assoc_legendre_q (n m x)
  (cond ((and (integerp n) (> n -1) (integerp m) (<= (abs m) n))
         (cond ((< m 0)
                (mul (div (factorial (+ n m)) (factorial (- n m)))
                     ($assoc_legendre_q n (- m) x)))
               (t
                (if (not (or (floatp x) ($bfloatp x))) (setq x ($rat x)))
                (let* ((q0 (q0m m x))
                       (q1 (if (= n 0) q0 (q1m m x)))
                       (q) (i 2)
                       (use-rat (or ($ratp x) (floatp x) ($bfloatp x))))
                  
                  (while (<= i n)
                    (setq q (if (= i m) (assoc-legendre-q-nn i x) 
                              (div (sub (mul (- (* 2 i) 1) x q1) 
                                        (mul (+ i -1 m) q0)) (- i m))))
                    (setq q0 q1)
                    (setq q1 q)
                    (incf i))
                  (if use-rat q1 ($ratsimp q1))))))
        (t `(($assoc_legendre_q simp) ,n ,m ,x))))

;; See G & R, 8.733, page 1005 first equation.

(putprop '$assoc_legendre_q
         '((n m x)
           ((unk) "$first" "$assoc_legendre_q")
           ((unk) "$second" "$assoc_legendre_q")
           
           ((mplus)
            ((mtimes)
             ((mplus)
              ((mtimes) -1 ((mplus) 1 ((mtimes) -1 m) n)
               (($assoc_legendre_q ) ((mplus ) 1 n) m x))
              ((mtimes) ((mplus) 1 n)
               (($assoc_legendre_q ) n m x) x))
             ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))) -1))))
         'grad) 
            
(defprop $assoc_legendre_q tex-assoc-legendre-q tex)

(defun tex-assoc-legendre-q (x l r)
  (tex-sub-and-super-scripted-function "Q" `(0) nil `(1) nil 2 x l r))

(setf (get '$assoc_legendre_q 'dimension) 'dimension-assoc-legendre-q)

(defun dimension-assoc-legendre-q (form result)
 (dimension-function
   (dimension-sub-and-super-scripted-function "Q" `(1) `(2) nil 3 form)
   result))

;; See A & S 22.5.37 page 779, A & S 8.6.6 (second equation) page 334, and 
;; A & S 8.2.5 page 333.  For n < 0, see A&S 8.2.1 page 333.

(defun $assoc_legendre_p (n m x)
  (let ((f) (d) (dx 0))
    (cond ((and (integerp n) (integerp m))
           (cond ((< n 0)
                  (setq f ($assoc_legendre_p (- (abs n) 1) m x))
                  (setq d 1)
                  (setq dx 1))
                 ((> (abs m) n)
                  (setq f 0)
                  (setq d 1))
                 ((< m 0)
                  (setq f ($assoc_legendre_p n (neg m) x))
                  ;; Adding a factor (-1)^m to the transformation to get the
                  ;; expected results for odd negative integers. DK 09/2009
                  (setq d (mul (power -1 m)
                               (div (factorial (+ n m)) (factorial (- n m)))))
                  (setq dx 1))
                 (t
                  (cond ((eql m 0)
                         (setq d 1))
                        (t
                         (setq d (simplify  
                                  `((%genfact) ,(- (* 2 m) 1) ,(- m 1) 2)))
                         (setq d (mul d (if (oddp m) -1 1)))
                         (setq d (mul d (power (sub 1 (mul x x)) (div m 2))))))
                  (setq dx 4)
                  (setq f 
                        ($ultraspherical (- n m) (add m (rat 1 2)) x)))))
          (t
           (setq d 1)
           (setq f `(($assoc_legendre_p simp) ,n ,m ,x))))
    (interval-mult d f (* flonum-epsilon dx))))


;; For the derivative of the associated legendre p function, see
;; A & S 8.5.4 page 334.

(putprop `$assoc_legendre_p
         '((n m x)
           ((unk) "$first" "$assoc_legendre_p")
           ((unk) "$second" "$assoc_legendre_p")
           ((mtimes simp)
            ((mplus simp)
             ((mtimes simp) -1 ((mplus simp) m n) (($unit_step) n)
              (($assoc_legendre_p simp) ((mplus simp) -1 n) m x))
             ((mtimes simp) n (($assoc_legendre_p simp) n m x) x))
            ((mexpt simp) ((mplus simp) -1 ((mexpt simp) x 2)) -1))) 
           'grad) 
           
(defprop $assoc_legendre_p tex-assoc-legendre-p tex)

(defun tex-assoc-legendre-p (x l r)
  (tex-sub-and-super-scripted-function "P" `(0) nil `(1) nil 2 x l r))

(setf (get '$assoc_legendre_p 'dimension) 'dimension-assoc-legendre-p)

(defun dimension-assoc-legendre-p (form result)
  (dimension-function
   (dimension-sub-and-super-scripted-function "P" `(1) `(2) nil 3 form)
   result))
                                
;; See A&S 22.5.55 and 22.5.56, page 780.

(defun $hermite (n x)
  (cond ((integerp n)
         (let ((f) (d) (e))
           (cond (($evenp n)
                  (setq n (/ n 2))
                  (multiple-value-setq (f e)  
                    ($hypergeo11 (- n) (rat 1 2) (mul x x) n))
                  (setq d (mul (if (oddp n) -1 1) (factorial (* 2 n))
                               (div 1 (factorial n)))))
                 (($oddp n)
                  (setq n (/ (- n 1) 2))
                  (multiple-value-setq (f e)
                    ($hypergeo11 (- n) (rat 3 2) (mul x x) n))
                  (setq d (mul (if (oddp n) -1 1) (factorial (+ 1 (* 2 n))) 2 x
                               (div 1 (factorial n))))))
           (orthopoly-return-handler d f e))) 
        (t `(($hermite) ,n ,x))))

(putprop '$hermite
         '((n x)
           ((unk) "$first" "$hermite")
           ((mtimes) 2 n (($hermite) ((mplus) -1 n) x)))
         'grad)

(defprop $hermite tex-hermite tex)

(defun tex-hermite (x l r)
  (tex-sub-and-super-scripted-function "H" `(0) nil nil nil 1 x l r))

(setf (get '$hermite 'dimension) 'dimension-hermite)

(defun dimension-hermite (form result)
 (dimension-function
   (dimension-sub-and-super-scripted-function "H" `(1) nil nil 2 form)
   result))

;; See A & S 22.5.54, page 780.  For integer n, use the identity
;;     binomial(n+a,n) = pochhammer(a+1,n)/pochhammer(1,n)

(defun $gen_laguerre (n a x)
  (cond ((use-hypergeo n x)
         (let ((f) (d) (e))
           ;(setq d (div ($pochhammer (add a 1) n) ($pochhammer 1 n)))
           (setq d (pochhammer-quotient (add a 1) 1 x n))
           (multiple-value-setq (f e)
             ($hypergeo11 (mul -1 n) (add 1 a) x n))
           (setq e (if e (+ e (* 4 (abs f) flonum-epsilon n)) nil))
           (orthopoly-return-handler d f e)))
        (t
         `(($gen_laguerre) ,n ,a ,x))))

(putprop '$gen_laguerre
         '((n a x)
           ((unk) "$first" "$gen_laguerre")
           ((unk) "$second" "$gen_laguerre")
           ((mtimes)
            ((mplus)
             ((mtimes) -1 ((mplus) a n)
              (($unit_step) n) (($gen_laguerre) ((mplus) -1 n) a x))
             ((mtimes) n (($gen_laguerre) n a x)))
            ((mexpt) x -1)))
         'grad)
                 
(defprop $gen_laguerre tex-gen-laguerre tex)

(defun tex-gen-laguerre (x l r)
  (tex-sub-and-super-scripted-function "L" `(0) nil `(1) t 1 x l r))

(setf (get '$gen_laguerre 'dimension) 'dimension-gen-laguerre)

(defun dimension-gen-laguerre (form result)
  (dimension-function
   (dimension-sub-and-super-scripted-function "L" `(1) `(2) t 3 form)
   result))

;; See A & S 22.5.16, page 778.

(defun $laguerre (n x)
  (cond ((use-hypergeo n x)
         (let ((f) (e))
           (multiple-value-setq (f e) ($hypergeo11 (mul -1 n) 1 x n))
           (orthopoly-return-handler 1 f e)))
        (t
         `(($laguerre) ,n ,x))))

(putprop '$laguerre
         '((n x)
           ((unk) "$first" "$laguerre")
           ((mtimes)
            ((mplus)
             ((mtimes) -1 n (($laguerre) ((mplus) -1 n) x))
             ((mtimes) n (($laguerre) n x)))
            ((mexpt) x -1)))
         'grad)
         
(defprop $laguerre tex-laguerre tex)

(defun tex-laguerre (x l r)
  (tex-sub-and-super-scripted-function "L" `(0) nil nil nil 1 x l r))

(setf (get '$laguerre 'dimension) 'dimension-laguerre)

(defun dimension-laguerre (form result)
  (dimension-function
   (dimension-sub-and-super-scripted-function "L" `(1) nil nil 2 form)
   result))

(defun $spherical_hankel1 (n x)
  (let ((f) (d) (e))
    (cond ((and (integerp n) (< n 0))
           (setq d (mul '$%i (if (oddp n) 1 -1)))
           (multiple-value-setq (f e)
             ($spherical_hankel1 (add -1 (mul -1 n)) x))
           (orthopoly-return-handler d f e))
          ((use-hypergeo n x)
           (multiple-value-setq (f e)
             ($hypergeo11 (mul -1 n) (mul -2 n) (mul -2 '$%i x) n))
           (setq d (mul '$%i (if (= 0 n) 1 
                               (simplify `((%genfact) ,(add (mul 2 n) -1) 
                                           ,(add n (rat -1 2)) 2)))
                        (power '$%e (mul '$%i x)) (div -1 (power x (add 1 n)))))
           (orthopoly-return-handler d f e))
          (t 
           `(($spherical_hankel1) ,n ,x)))))

(putprop '$spherical_hankel1
         '((n x)
           ((unk) "$first" "$spherical_hankel1")
           ((mplus simp) (($spherical_hankel1) ((mplus) -1 n) x)
            ((mtimes simp) -1 ((mplus) 1 n)
             (($spherical_hankel1) n x) ((mexpt) x -1))))
         'grad)

(defprop $spherical_hankel1 tex-spherical-hankel-1 tex)

(defun tex-spherical-hankel-1 (x l r)
  (tex-sub-and-super-scripted-function "h^{(1)}" `(0) nil nil nil 1 x l r))

(setf (get '$spherical_hankel1 'dimension) 'dimension-spherical-hankel-1)

(defun dimension-spherical-hankel-1 (form result)
  (let ((form1 `((mexpt) (($\h simp array) ,(nth 1 form)) 
                 (1))))
    (dimension-function `((,form1 simp) ,(nth 2 form)) result)))

;; See A & S 10.1.36.

(defun $spherical_hankel2 (n x)
  (cond ((integerp n)
         (setq x (mul x (power '$%e (mul '$%i '$%pi (add (mul 2 n) 1)))))
         (let ((f))
           (setq f ($spherical_hankel1 n x))
           (if (oddp n) (interval-mult -1 f) f)))
        (t `(($spherical_hankel2) ,n ,x))))

(putprop '$spherical_hankel2
         '((n x)
           ((unk) "$first" "$spherical_hankel2")
           ((mplus simp) (($spherical_hankel2) ((mplus) -1 n) x)
            ((mtimes simp) -1 ((mplus) 1 n)
             (($spherical_hankel2) n x) ((mexpt) x -1))))
         'grad)

(defprop $spherical_hankel2 tex-spherical-hankel-2 tex)

(defun tex-spherical-hankel-2 (x l r)
  (tex-sub-and-super-scripted-function "h^{(2)}" `(0) nil nil nil 1 x l r))

(setf (get '$spherical_hankel2 'dimension) 'dimension-spherical-hankel-2)

(defun dimension-spherical-hankel-2 (form result)
  (let ((form1 `((mexpt) (($\h simp array) ,(nth 1 form))  (2))))
    (dimension-function `((,form1 simp) ,(nth 2 form)) result)))
  
;;---------------------------------------------------------------------
;; The spherical_bessel functions use the functions p-fun and q-fun.
;; See A&S 10.1.8 and 10.1.9 page 437.

(defun p-fun (n x)
  (let ((s 1) (w 1) 
        (n1 (floor (/ n 2)))
        (x2 (mul x x)) (m2))
    (dotimes (m n1 s)
      (setq m2 (* 2 m))
      (setq w (div (mul w `((rat) ,(* -1 (+ n m2 2) (+ n m2 1) 
                                      (- n m2) (- n (+ m2 1)))
                            ,(* 4 (+ m2 1) (+ m2 2)))) x2))
      (setq s (add s w)))))

(defun q-fun (n x)
  (let ((s (if (= 0 n) 0 1))
        (w 1) (m2) (x2 (mul x x))
        (n1 (floor (/ (- n 1) 2))))
    (dotimes (m n1 (div (mul n (+ n 1) s) (mul 2 x)))
      (setq m2 (* 2 m))
      (setq w (div (mul w `((rat) ,(* -1 (+ n m2 3) (+ n m2 2) 
                                      (- n (+ m2 1)) (- n (+ m2 2)))
                            ,(* 4 (+ m2 3) (+ m2 2)))) x2))
      (setq s (add s w)))))


;; See A&S 10.1.8 page 437 and A&S 10.1.15 page 439.  When the order
;; is an integer and x is a float or a bigfloat, use the slatec code
;; for numerical evaluation.  Yes, we coerce bigfloats to floats and
;; return a float.

;; For numerical evaluation, we do our own analytic continuation -- otherwise
;; we get factors exp(%i n %pi / 2) that should evaluate to 1,%i,-1,-%, but
;; numerically have "fuzz" in their values.  The fuzz can cause the spherical
;; bessel functions to have nonzero (but small) imaginary values on the
;; negative real axis. See A&S 10.1.34

(defun $spherical_bessel_j (n x)
  (cond ((and (eq '$zero (csign ($ratdisrep x)))
              (or (integerp n) ($featurep n '$integer)))
         `((%kron_delta) 0 ,n))

        ((and (use-float x) (integerp n))
         (let ((d 1) (xr) (xi) (z))
           (setq x ($rectform ($float x)))
           (setq xr ($realpart x))
           (setq xi ($imagpart x))
           (setq z (complex xr xi))
           (cond ((< xr 0.0)
                  (setq d (if (oddp n) -1 1))
                  (setq x (mul -1 x))
                  (setq z (* -1 z))))
           (setq n (+ 0.5 ($float n)))
           (setq d (* d (sqrt (/ pi (* 2 z)))))
           (setq d (lisp-float-to-maxima-float d))
           ($expand (mul ($rectform d) ($bessel_j n x)))))

        ((and (integerp n) (> n -1))
         (let ((xt (sub x (div (mul n '$%pi) 2))))
           (div (add
                 (mul (p-fun n x) (simplify `((%sin) ,xt)))
                 (mul (q-fun n x) (simplify `((%cos) ,xt)))) x)))

        ((integerp n)
         (mul (if (oddp n) -1 1) ($spherical_bessel_y (- (+ n 1)) x)))

        (t 
         `(($spherical_bessel_j) ,n ,x))))
         
(putprop '$spherical_bessel_j
         '((n x)
           ((unk) "$first" "$spherical_bessel_j")
           ((mtimes) ((mexpt) ((mplus) 1 ((mtimes) 2 n)) -1)
            ((mplus)
             ((mtimes) n (($spherical_bessel_j) ((mplus) -1 n) x))
             ((mtimes) -1 ((mplus) 1 n)
              (($spherical_bessel_j) ((mplus) 1 n) x)))))
         'grad)
 
(defprop $spherical_bessel_j tex-spherical-bessel-j tex)

(defun tex-spherical-bessel-j (x l r)
  (tex-sub-and-super-scripted-function "j^{(2)}" `(0) nil nil nil 1 x l r))

(setf (get '$spherical_bessel_j 'dimension) 'dimension-spherical-bessel-j)

(defun dimension-spherical-bessel-j (form result)
  (let ((form1 `(($\j simp array) ,(nth 1 form)))) 
    (dimension-function `((,form1 simp) ,(nth 2 form)) result)))

;; For analytic continuation, see A&S 10.1.35.
 
(defun $spherical_bessel_y (n x)
  (cond ((and (use-float x) (integerp n))
         (let ((d 1) (xr) (xi) (z))
           (setq x ($rectform ($float x)))
           (setq xr ($realpart x))
           (setq xi ($imagpart x))
           (setq z (complex xr xi))
           (cond ((< xr 0.0)
                  (setq d (if (oddp n) 1 -1))
                  (setq x (mul -1 x))
                  (setq z (* -1 z))))
           (setq n (+ 0.5 ($float n)))
           (setq d (* d (sqrt (/ pi (* 2 z)))))
           (setq d (lisp-float-to-maxima-float d))
           ($expand (mul ($rectform d) ($bessel_y n x)))))

        ((and (integerp n) (> n -1))
         (let ((xt (add x (div (mul n '$%pi) 2))))
           (mul (if (oddp n) 1 -1)
                (div (sub
                      (mul (p-fun n x) (simplify `((%cos) ,xt)))
                      (mul (q-fun n x) (simplify `((%sin) ,xt)))) x))))

        ((integerp n)
         (mul (if (oddp n) 1 -1) ($spherical_bessel_j (- (+ n 1)) x)))
        (t  `(($spherical_bessel_y) ,n ,x))))

(putprop '$spherical_bessel_y
         '((n x)
           ((unk) "$first" "$spherical_bessel_y")
           ((mtimes) ((mexpt) ((mplus) 1 ((mtimes) 2 n)) -1)
            ((mplus)
             ((mtimes) n (($spherical_bessel_y) ((mplus) -1 n) x))
             ((mtimes) -1 ((mplus) 1 n)
              (($spherical_bessel_y) ((mplus) 1 n) x)))))
         'grad)
 
(defprop $spherical_bessel_y tex-spherical-bessel-y tex)

(defun tex-spherical-bessel-y (x l r)
  (tex-sub-and-super-scripted-function "y^{(2)}" `(0) nil nil nil 1 x l r))

 (setf (get '$spherical_bessel_y 'dimension) 'dimension-spherical-bessel-y)

(defun dimension-spherical-bessel-y (form result)
  (let ((form1 `(($\y simp array) ,(nth 1 form)))) 
    (dimension-function `((,form1 simp) ,(nth 2 form)) result)))

;; Compute P_n^m(cos(theta)).  See Merzbacher, 9.59 page 184
;; and 9.64 page 185, and A & S 22.5.37 page 779.  This function
;; lacks error checking; it should only be called by spherical_harmonic.

;; We need to be careful -- for the spherical harmonics we can't use
;; assoc_legendre_p(n,m,cos(x)).  If we did, we'd get factors 
;; (1 - cos^2(x))^(m/2) that simplify to |sin(x)|^m but we want them
;; to simplify to sin^m(x).  Oh my!

(defun assoc-leg-cos (n m x)
  (interval-mult
   (if (= m 0) 1 (mul (take '(%genfact) (sub (mul 2 m) 1) (sub m (div 1 2)) 2) (power (take '(%sin) x) m)))
   ($ultraspherical (sub n m) (add m (div 1 2)) (take  '(%cos) x))))

(defun $spherical_harmonic (n m th p)
  (cond ((and (integerp n) (integerp m) (> n -1))
         (cond ((> (abs m) n)
                0)
               ((< m 0)
                (interval-mult (if (oddp m) -1 1) 
                               ($spherical_harmonic n (- m) th (mul -1 p))))
               (t
                (interval-mult
                 (mul ($exp (mul '$%i m p))
                      (power (div (* (+ (* 2 n) 1) (factorial (- n m)))
                                  (mul '$%pi (* 4 (factorial (+ n m))))) 
                             `((rat) 1 2)))
                 (assoc-leg-cos n m th)))))
        (t
         `(($spherical_harmonic) ,n ,m ,th ,p))))

(defprop $spherical_harmonic tex-spherical-harmonic tex)

(defun tex-spherical-harmonic (x l r)
  (tex-sub-and-super-scripted-function "Y" `(0) nil `(1) nil 2 x l r))

(setf (get '$spherical_harmonic 'dimension) 'dimension-spherical-harmonic)

(defun dimension-spherical-harmonic (form result)
 (dimension-function
  (dimension-sub-and-super-scripted-function "Y" `(1) `(2) nil 3 form)
  result))

(putprop '$spherical_harmonic
         '((n m theta phi)
           ((unk) "$first" "$spherical_harmonic")
           ((unk) "$second" "$spherical_harmonic")
           ((mplus)
            ((mtimes) ((rat ) -1 2)
             ((mexpt)
              ((mtimes) ((mplus) ((mtimes) -1 m) n)
               ((mplus) 1 m n))
              ((rat) 1 2))
             (($spherical_harmonic) n ((mplus) 1 m) theta phi)
             ((mexpt) $%e ((mtimes) -1 $%i phi)))
            ((mtimes) ((rat) 1 2)
             ((mexpt)
              ((mtimes) ((mplus) 1 ((mtimes) -1 m) n)
               ((mplus) m n))
              ((rat) 1 2))
             (($spherical_harmonic) n ((mplus) -1 m) theta phi)
             ((mexpt) $%e ((mtimes) $%i phi)))) 
           
           ((mtimes) $%i m (($spherical_harmonic) n m theta phi)))
         'grad)
                                                                        
(defun maxima-float-to-lisp-float (y)
  (let* ((x ($rectform y))
         (xr ($realpart x))
         (xi ($imagpart x)))
    (cond ((or (floatp xr) (floatp xi))
           (setq xr (float xr)
                 xi (float xi))
           (if (= 0.0 xi) xr (complex xr xi)))
          (t
           y))))

(defun lisp-float-to-maxima-float (x)
  (if (complexp x)
      (add (realpart x) (mul '$%i (imagpart x)))
    x))
     
(defun hypergeo11-float (n b x)
  (let ((f0) (fm1) (f) (i 0) (k) (dk) (ak) (bk) (err)
        (as (make-array (- 1 n) 
                        :initial-element 0.0))
        (bs (make-array (- 1 n) 
                        :initial-element 0.0))
        (fs (make-array (- 1 n)
                        :initial-element 0.0))
        (u) (u0) (um1))

    (setq f0 1.0)
    (setq fm1 0.0)
    (setq x (- b x))
    (setq n (- n))
    (while (< i n)
      (setf (aref fs i) f0)
      (setq dk (+ b i))
      (setq ak (/ (+ (* 2 i) x) dk))
      (setq bk (/ (- i) dk))
      (setq f (+ (* ak f0) (* bk fm1)))
      (setf (aref as i) ak)
      (setf (aref bs i) bk)
      (setq fm1 f0)
      (setq f0 f)
      (incf i))
    (setf (aref fs i) f0)
    (setq i 1)
    (setq err 1.0)
    (setq u0 1.0)
    (setq um1 0.0)
    (while (< i n)
      (setq k (- n i))
      (setq u (+ (* (aref as k) u0)
                 (* (aref bs (+ 1 k)) um1)))
      (setq um1 u0)
      (setq u0 u)
      (setq err (+ err (abs (* u0 (aref fs k)))))
      (incf i))
    (values f0 (* 12 flonum-epsilon err))))
    
(defun hypergeo21-float (n b c x)
  (let ((f0) (fm1) (f) (i 0) (k) (dk) (ak) (bk) (err)
        (as (make-array (- 1 n) 
                        :initial-element 0.0))
        (bs (make-array (- 1 n) 
                        :initial-element 0.0))
        (fs (make-array (- 1 n)
                        :initial-element 0.0))
        (u) (u0) (um1))

    (setq f0 1.0)
    (setq fm1 0.0)
    (setq n (- n))
    (while (< i n)
      (setf (aref fs i) f0)
      (setq dk (+ c i))
      (setq ak (/ (+ (* 2 i) c (- (* x (+ b i)))) dk))
      (setq bk (/ (* i (- x 1)) dk))
      (setq f (+ (* ak f0) (* bk fm1)))
      (setf (aref as i) ak)
      (setf (aref bs i) bk)
      (setq fm1 f0)
      (setq f0 f)
      (incf i))
    (setf (aref fs i) f0)
    (setq i 1)
    (setq err 1.0)
    (setq u0 1.0)
    (setq um1 0.0)
    (while (< i n)
      (setq k (- n i))
      (setq u (+ (* (aref as k) u0)
                 (* (aref bs (+ 1 k)) um1)))
      (setq um1 u0)
      (setq u0 u)
      (incf i)
      (setq err (+ err (abs (* (aref fs k) u0)))))
    (values f0 (* 12 flonum-epsilon err))))
    
;; For recursion relations, see A & S 22.7 page 782. 

;; legendre_p(n+1,x) = ((2*n+1)*legendre_p(n,x)*x-n*legendre_p(n-1,x))/(n+1)

;;  jacobi_p(n+1,a,b,x) = (((2*n+a+b+1)*(a^2-b^2) + 
;;    x*pochhammer(2*n+a+b,3)) * jacobi_p(n,a,b,x) - 
;;    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))

;; ultraspherical(n+1,a,x) = (2*(n+a)*x * ultraspherical(n,a,x) - 
;;    (n+2*a-1)*ultraspherical(n-1,a,x))/(n+1)

;; chebyshev_t(n+1,x) = 2*x*chebyshev_t(n,x) -chebyshev_t(n-1,x)

;; chebyshev_u(n+1,x) = 2*x*chebyshev_u(n,x) -chebyshev_u(n-1,x)

;;  laguerre(n+1,x) = (((2*n+1) - x)*laguerre(n,x)  -(n)*laguerre(n-1,x))/(n+1)

;; gen_laguerre(n+1,a,x) = (((2*n+a+1) - x)*gen_laguerre(n,a,x)  
;;  -(n+a)*gen_laguerre(n-1,a,x))/(n+1)

;; hermite(n+1,x) = 2*x*hermite(n,x) -2*n*hermite(n-1,x)

;; See G & R 8.733.2; A & S 22.7.11 might be wrong -- or maybe I need
;; reading glasses.

;; (2*n+1)*x*assoc_legendre_p(n,m,x) = (n-m+1)*assoc_legendre_p(n+1,m,x) 
;; + (n+m)*assoc_legendre_p(n-1,m,x)

;; For the half-integer bessel functions, See A & S 10.1.19

;; fn(n-1,x) + fn(n+1,x) = (2*n+1)*fn(n,x)/x;

(defun check-arg-length (fn n m)
  (if (not (= n m))
      (merror "Function ~:M needs ~:M arguments, instead it received ~:M"
              fn n m)))

(defun $orthopoly_recur (fn arg)
  (if (not ($listp arg)) 
      (merror "The second argument to orthopoly_recur must be a list"))
  (cond ((eq fn '$jacobi_p)
         (check-arg-length fn 4 (- (length arg) 1))
         (let ((n (nth 1 arg))
               (a (nth 2 arg))
               (b (nth 3 arg))
               (x (nth 4 arg)))
           (simplify
            `((mequal) (($jacobi_p ) ((mplus) 1 ,n) ,a ,b ,x)
              ((mtimes) ((rat) 1 2) ((mexpt) ((mplus) 1 ,n) -1)
               ((mexpt) ((mplus) 1 ,a ,b ,n) -1)
               ((mexpt) ((mplus) ,a ,b ((mtimes) 2 ,n)) -1)
               ((mplus)
                ((mtimes) -2 ((mplus) ,a ,n) ((mplus) ,b ,n)
                 ((mplus) 2 ,a ,b ((mtimes) 2 ,n))
                 (($jacobi_p ) ((mplus) -1 ,n) ,a ,b ,x))
                ((mtimes) (($jacobi_p ) ,n ,a ,b ,x)
                 ((mplus)
                 ((mtimes)
                  ((mplus) ((mexpt) ,a 2)
                   ((mtimes) -1 ((mexpt) ,b 2)))
                  ((mplus) 1 ,a ,b ((mtimes) 2 ,n)))
                 ((mtimes) ((mplus) ,a ,b ((mtimes) 2 ,n))
                  ((mplus) 1 ,a ,b ((mtimes) 2 ,n))
                  ((mplus) 2 ,a ,b ((mtimes) 2 ,n)) ,x)))))))))

        ((eq fn '$ultraspherical)
         (check-arg-length fn 3 (- (length arg) 1))
         (let ((n (nth 1 arg))
               (a (nth 2 arg))
               (x (nth 3 arg)))
           (simplify
            `((mequal) (($ultraspherical) ((mplus) 1 ,n) ,a ,x)
             ((mtimes) ((mexpt) ((mplus) 1 ,n) -1)
              ((mplus)
               ((mtimes) -1 ((mplus) -1 ((mtimes) 2 ,a) ,n)
                (($ultraspherical) ((mplus) -1 ,n) ,a ,x))
               ((mtimes) 2 ((mplus) ,a ,n)
                (($ultraspherical) ,n ,a ,x) ,x)))))))

        ((member fn  `($chebyshev_t $chebyshev_u) :test 'eq)
         (check-arg-length fn 2 (- (length arg) 1))
         (let ((n (nth 1 arg))
               (x (nth 2 arg)))
          (simplify
           `((mequal ) ((,fn) ((mplus ) 1 ,n) ,x)
            ((mplus )
             ((mtimes ) -1 ((,fn) ((mplus ) -1 ,n) ,x))
             ((mtimes ) 2 ((,fn) ,n ,x) ,x))))))

        ((member fn '($legendre_p $legendre_q) :test 'eq)
         (check-arg-length fn 2 (- (length arg) 1))
         (let* ((n (nth 1 arg))
               (x (nth 2 arg))
               (z (if (eq fn '$legendre_q) 
                      `((mtimes) -1 ((%kron_delta) ,n 0)) 0))) 
           (simplify
             `((mequal) ((,fn) ((mplus) 1 ,n) ,x)
               ((mplus)
                ((mtimes) ((mexpt) ((mplus) 1 ,n) -1)
                 ((mplus)
                  ((mtimes) ((mtimes) -1 ,n)
                   ((,fn) ((mplus) -1 ,n) ,x))
                  ((mtimes) ((mplus) 1 ((mtimes) 2 ,n))
                   ((,fn) ,n ,x) ,x)))
                ,z)))))

        ((member fn '($assoc_legendre_p $assoc_legendre_q) :test 'eq)
         (check-arg-length fn 3 (- (length arg) 1))
         (let ((n (nth 1 arg))
               (m (nth 2 arg))
               (x (nth 3 arg)))
           (simplify
            `((mequal) ((,fn) ((mplus) 1 ,n) ,m ,x)
              ((mtimes)
               ((mexpt) ((mplus) 1 ((mtimes) -1 ,m) ,n) -1)
               ((mplus)
                ((mtimes)
                 ((mplus) ((mtimes) -1 ,m)
                  ((mtimes) -1 ,n))
                 ((,fn) ((mplus) -1 ,n) ,m ,x))
                ((mtimes) ((mplus) 1 ((mtimes) 2 ,n))
                 ((,fn) ,n ,m ,x) ,x))))))) 
        
        ((eq fn '$laguerre)
         (check-arg-length fn 2 (- (length arg) 1))
         (let ((n (nth 1 arg))
               (x (nth 2 arg)))
           (simplify
            `((mequal ) (($laguerre ) ((mplus ) 1 ,n) ,x)
              ((mtimes ) ((mexpt ) ((mplus ) 1 ,n) -1)
               ((mplus )
                ((mtimes ) -1 ,n (($laguerre ) ((mplus ) -1 ,n) ,x))
                ((mtimes ) (($laguerre ) ,n ,x)
                 ((mplus ) 1 ((mtimes ) 2 ,n) ((mtimes ) -1 ,x))))))))) 

        ((eq fn '$gen_laguerre)
         (check-arg-length fn 3 (- (length arg) 1))
         (let ((n (nth 1 arg))
               (a (nth 2 arg))
               (x (nth 3 arg)))
           (simplify
            `((mequal) (($gen_laguerre) ((mplus) 1 ,n) ,a ,x)
              ((mtimes) ((mexpt ) ((mplus) 1 ,n) -1)
               ((mplus)
                ((mtimes) -1 ((mplus) ,a ,n)
                 (($gen_laguerre) ((mplus) -1 ,n) ,a ,x))
                ((mtimes) (($gen_laguerre) ,n ,a ,x)
                 ((mplus) 1 ,a ((mtimes) 2 ,n) ((mtimes ) -1 ,x))))))))) 

        ((eq fn '$hermite)
         (check-arg-length fn 2 (- (length arg) 1))
         (let ((n (nth 1 arg))
               (x (nth 2 arg)))
           (simplify
            `((mequal) (($hermite) ((mplus) 1 ,n) ,x)
              ((mplus)
               ((mtimes) -2 ,n (($hermite) ((mplus) -1 ,n) ,x))
               ((mtimes) 2 (($hermite) ,n ,x) ,x))))))

        ((member fn `($spherical_bessel_j $spherical_bessel_y
                                          $spherical_hankel1 $spherical_hankel2)
                 :test 'eq)
         (check-arg-length fn 2 (- (length arg) 1))
         (let ((n (nth 1 arg))
               (x (nth 2 arg)))
           (simplify
            `((mequal) ((,fn) ((mplus) 1 ,n) ,x)
              ((mplus) ((mtimes) -1 ((,fn ) ((mplus) -1 ,n) ,x))
               ((mtimes) ((,fn) ,n ,x) ((mexpt) ,x -1))
               ((mtimes) 2 ,n ((,fn ) ,n ,x) ((mexpt) ,x -1)))))))
         
        (t (merror "A recursion relation for ~:M isn't known to Maxima" fn))))
    
;; See A & S Table 22.2, page 774.

(defun $orthopoly_weight (fn arg)
  (if (not ($listp arg)) 
      (merror "The second argument to orthopoly_weight must be a list"))

  (if (not (or ($symbolp (car (last arg))) ($subvarp (car (last arg)))))
      (merror "The last element of the second argument to orthopoly_weight must
be a symbol or a subscripted symbol, instead found ~:M" (car (last arg))))

  (if (not (every #'(lambda (s) 
                      ($freeof (car (last arg)) s)) (butlast (cdr arg))))
      (merror "Only the last element of ~:M may depend on the integration
variable ~:M" arg (car (last arg))))

  (cond ((eq fn '$jacobi_p)
         (check-arg-length fn 4 (- (length arg) 1))
         (let ((a (nth 2 arg))
               (b (nth 3 arg))
               (x (nth 4 arg)))
           (simplify
            `((mlist)
              ((mtimes) ((mexpt) ((mplus) 1 ((mtimes) -1 ,x)) ,a)
               ((mexpt) ((mplus ) 1 ,x) ,b))
              -1 1))))
        
        ((eq fn '$ultraspherical)
         (check-arg-length fn 3 (- (length arg) 1))
         (let ((a (nth 2 arg))
               (x (nth 3 arg)))
           (simplify
            `((mlist)
              ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) ,x 2)))
               ((mplus) ((rat) -1 2) ,a)) -1 1))))

        ((eq fn '$chebyshev_t)
         (check-arg-length fn 2 (- (length arg) 1))
         (let ((x (nth 2 arg)))
           (simplify
            `((mlist)
              ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) ,x 2)))
               ((rat) -1 2)) -1 1)))) 
          
        ((eq fn '$chebyshev_u)
         (check-arg-length fn 2 (- (length arg) 1))
         (let ((x (nth 2 arg)))
           (simplify
            `((mlist)
              ((mexpt) ((mplus) 1  ((mtimes) -1 ((mexpt) ,x 2)))
               ((rat) 1 2)) -1 1))))

        ((eq fn '$legendre_p)
         (check-arg-length fn 2 (- (length arg) 1))
         `((mlist) 1 -1 1))

        ; This is for a fixed order.  There is also an orthogonality
        ; condition for fixed degree with weight function 1/(1-x^2).
        ; See A & S 8.14.11 and 8.14.12.
        ((eq fn '$assoc_legendre_p)
         (check-arg-length fn 3 (- (length arg) 1))
         `((mlist) 1 -1 1))

        ((eq fn '$laguerre)
         (check-arg-length fn 2 (- (length arg) 1))
         (let ((x (nth 2 arg)))
           (simplify
            `((mlist) ((mexpt) $%e ((mtimes) -1 ,x)) 0 $inf))))

        ((eq fn '$gen_laguerre)
         (check-arg-length fn 3 (- (length arg) 1))
         (let ((a (nth 2 arg))
               (x (nth 3 arg)))
           (simplify
            `((mlist)
              ((mtimes) ((mexpt) ,x ,a)
               ((mexpt) $%e ((mtimes) -1 ,x))) 0 $inf))))

        ((eq fn '$hermite)
         (check-arg-length fn 2 (- (length arg) 1))
         (let ((x (nth 2 arg)))
           (simplify
            `((mlist) ((mexpt) $%e ((mtimes) -1 ((mexpt) ,x 2)))
              ((mtimes ) -1 $inf) $inf))))

        (t (merror "A weight for ~:M isn't known to Maxima" fn))))