;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10-*- ;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;     The data in this file contains enhancments.                    ;;;;;
;;;                                                                    ;;;;;
;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
;;;     All rights reserved                                            ;;;;;
;;;                                                                    ;;;;;
;;;  Copyright (c) 2001 by Raymond Toy.  Replaced everything and added ;;;;;
;;;     support for symbolic manipulation of all 12 Jacobian elliptic  ;;;;;
;;;     functions and the complete and incomplete elliptic integral    ;;;;;
;;;     of the first, second and third kinds.                          ;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(in-package :maxima)
;;(macsyma-module ellipt)

(defvar 3//2 '((rat simp) 3 2))
(defvar 1//2 '((rat simp) 1 2))
(defvar -1//2 '((rat simp) -1 2))

;;;
;;; Jacobian elliptic functions and elliptic integrals.
;;;
;;; References:
;;;
;;; [1] Abramowitz and Stegun
;;; [2] Lawden, Elliptic Functions and Applications, Springer-Verlag, 1989
;;;
;;; We use the definitions from Abramowitz and Stegun where our
;;; sn(u,m) is sn(u|m).  That is, the second arg is the parameter,
;;; instead of the modulus k or modular angle alpha.
;;;
;;; Note that m = k^2 and k = sin(alpha).
;;;

;;
;; Routines for computing the basic elliptic functions sn, cn, and dn.
;;
;;
;; A&S gives several methods for computing elliptic functions
;; including the AGM method (16.4) and ascending and descending Landen
;; transformations (16.12 and 16.14).  We use these latter because
;; they are actually quite fast, only requiring simple arithmetic and
;; square roots for the transformation until the last step.  The AGM
;; requires evaluation of several trignometric functions at each
;; stage.
;;
;; In addition, the Landen transformations are valid for all u and m.
;; Thus, we can compute the elliptic functions for complex u and
;; m. (The code below supports this, but we could make it run much
;; faster if we specialized it to for flonums.  However, if we
;; do that, we won't be able to handle the cases where m < 0 or m > 1.
;; We'll have to handle these specially via A&S 16.10 and 16.11.)
;;
;; See A&S 16.12 and 16.14.



(flet ((ascending-transform (u m)
	 ;; A&S 16.14.1
	 ;;
	 ;; Take care in computing this transform.  For the case where
	 ;; m is complex, we should compute sqrt(mu1) first as
	 ;; (1-sqrt(m))/(1+sqrt(m)), and then square this to get mu1.
	 ;; If not, we may choose the wrong branch when computing
	 ;; sqrt(mu1).
	 (let* ((root-m (cl:sqrt m))
		(mu (/ (* 4 root-m)
		       (cl:expt (1+ root-m) 2)))
		(root-mu1 (/ (- 1 root-m) (+ 1 root-m)))
		(v (/ u (1+ root-mu1))))
	   (values v mu root-mu1)))
       (descending-transform (u m)
	 ;; Note: Don't calculate mu first, as given in 16.12.1.  We
	 ;; should calculate sqrt(mu) = (1-sqrt(m1)/(1+sqrt(m1)), and
	 ;; then compute mu = sqrt(mu)^2.  If we calculate mu first,
	 ;; sqrt(mu) loses information when m or m1 is complex.
	 (let* ((root-m1 (cl:sqrt (- 1 m)))
		(root-mu (/ (- 1 root-m1) (+ 1 root-m1)))
		(mu (* root-mu root-mu))
		(v (/ u (1+ root-mu))))
	   (values v mu root-mu))))
  (declare (inline descending-transform ascending-transform))


  ;; Could use the descending transform, but some of my tests show
  ;; that it has problems with roundoff errors.
  (defun elliptic-dn-ascending (u m)
    (cond ((zerop m)
	   ;; A&S 16.6.3
	   1.0)
	  ((< (abs (- 1 m)) (* 4 flonum-epsilon))
	   ;; A&S 16.6.3
	   (/ (cl:cosh u)))
	  (t
	   (multiple-value-bind (v mu root-mu1)
	       (ascending-transform u m)
	     ;; A&S 16.14.4
	     (let* ((new-dn (elliptic-dn-ascending v mu)))
	       (* (/ (- 1 root-mu1) mu)
		  (/ (+ root-mu1 (* new-dn new-dn))
		     new-dn)))))))

  ;; Don't use the descending version because it requires cn, dn, and
  ;; sn.
  (defun elliptic-cn-ascending (u m)
    (cond ((zerop m)
	   ;; A&S 16.6.2
	   (cl:cos u))
	  ((< (abs (- 1 m)) (* 4 flonum-epsilon))
	   ;; A&S 16.6.2
	   (/ (cl:cosh u)))
	  (t
	   (multiple-value-bind (v mu root-mu1)
	       (ascending-transform u m)
	     ;; A&S 16.14.3
	     (let* ((new-dn (elliptic-dn-ascending v mu)))
	       (* (/ (+ 1 root-mu1) mu)
		  (/ (- (* new-dn new-dn) root-mu1)
		     new-dn)))))))

  ;; We don't use the ascending transform here because it requires
  ;; evaluating sn, cn, and dn.  The ascending transform only needs
  ;; sn.
  (defun elliptic-sn-descending (u m)
    (cond ((= m 1)
	   ;; A&S 16.6.1
	   (cl:tanh u))
	  ((< (abs m) flonum-epsilon)
	   ;; A&S 16.6.1
	   (cl:sin u))
	  (t
	   (multiple-value-bind (v mu root-mu)
	       (descending-transform u m)
	     (let* ((new-sn (elliptic-sn-descending v mu)))
	       (/ (* (1+ root-mu) new-sn)
		  (1+ (* root-mu new-sn new-sn))))))))
  #+nil
  (defun elliptic-sn-ascending (u m)
    (if (< (abs (- 1 m)) (* 4 flonum-epsilon))
	;; A&S 16.6.1
	(tanh u)
	(multiple-value-bind (v mu root-mu1)
	    (ascending-transform u m)
	  ;; A&S 16.14.2
	  (let* ((new-cn (elliptic-cn-ascending v mu))
		 (new-dn (elliptic-dn-ascending v mu))
		 (new-sn (elliptic-sn-ascending v mu)))
	    (/ (* (+ 1 root-mu1) new-sn new-cn)
	       new-dn)))))
  )

(defun sn (u m)
  (if (and (realp u) (realp m))
      (realpart (elliptic-sn-descending u m))
      (elliptic-sn-descending u m)))

(defun cn (u m)
  (if (and (realp u) (realp m))
      (realpart (elliptic-cn-ascending u m))
      (elliptic-cn-ascending u m)))

(defun dn (u m)
  (if (and (realp u) (realp m))
      (realpart (elliptic-dn-ascending u m))
      (elliptic-dn-ascending u m)))


;;
;; How this works, I think.
;;
;; $jacobi_sn is the user visible function JACOBI_SN.  We put
;; properties on this symbol so maxima can figure out what to do with
;; it.

;; Tell maxima how to simplify the functions $jacobi_sn, etc.  This
;; borrows heavily from trigi.lisp.
(defprop %jacobi_sn simp-%jacobi_sn operators)
(defprop %jacobi_cn simp-%jacobi_cn operators)
(defprop %jacobi_dn simp-%jacobi_dn operators)
(defprop %inverse_jacobi_sn simp-%inverse_jacobi_sn operators)
(defprop %inverse_jacobi_cn simp-%inverse_jacobi_cn operators)
(defprop %inverse_jacobi_dn simp-%inverse_jacobi_dn operators)

;; Tell maxima what the derivatives are.
;;
;; Lawden says the derivative wrt to k but that's not what we want.
;;
;; Here's the derivation we used, based on how Lawden get's his results.
;;
;; Let
;;
;; diff(sn(u,m),m) = s
;; diff(cn(u,m),m) = p
;; diff(dn(u,m),m) = q
;;
;; From the derivatives of sn, cn, dn wrt to u, we have
;;
;; diff(sn(u,m),u) = cn(u)*dn(u)
;; diff(cn(u,m),u) = -cn(u)*dn(u)
;; diff(dn(u,m),u) = -m*sn(u)*cn(u)
;;

;; Differentiate these wrt to m:
;;
;; diff(s,u) = p*dn + cn*q
;; diff(p,u) = -p*dn - q*dn
;; diff(q,u) = -sn*cn - m*s*cn - m*sn*q
;;
;; Also recall that
;;
;; sn(u)^2 + cn(u)^2 = 1
;; dn(u)^2 + m*sn(u)^2 = 1
;;
;; Differentiate these wrt to m:
;;
;; sn*s + cn*p = 0
;; 2*dn*q + sn^2 + 2*m*sn*s = 0
;;
;; Thus,
;;
;; p = -s*sn/cn
;; q = -m*s*sn/dn - sn^2/dn/2
;;
;; So
;; diff(s,u) = -s*sn*dn/cn - m*s*sn*cn/dn - sn^2*cn/dn/2
;;
;; or
;;
;; diff(s,u) + s*(sn*dn/cn + m*sn*cn/dn) = -1/2*sn^2*cn/dn
;;
;; diff(s,u) + s*sn/cn/dn*(dn^2 + m*cn^2) = -1/2*sn^2*cn/dn
;;
;; Multiply through by the integrating factor 1/cn/dn:
;;
;; diff(s/cn/dn, u) = -1/2*sn^2/dn^2 = -1/2*sd^2.
;;
;; Interate this to get
;;
;; s/cn/dn = C + -1/2*int sd^2
;;
;; It can be shown that C is zero.
;;
;; We know that (by differentiating this expression)
;;
;; int dn^2 = (1-m)*u+m*sn*cd + m*(1-m)*int sd^2
;;
;; or
;;
;; int sd^2 = 1/m/(1-m)*int dn^2 - u/m - sn*cd/(1-m)
;;
;; Thus, we get
;;
;; s/cn/dn = u/(2*m) + sn*cd/(2*(1-m)) - 1/2/m/(1-m)*int dn^2
;;
;; or
;;
;; s = 1/(2*m)*u*cn*dn + 1/(2*(1-m))*sn*cn^2 - 1/2/(m*(1-m))*cn*dn*E(u)
;;
;; where E(u) = int dn^2 = elliptic_e(am(u)) = elliptic_e(asin(sn(u)))
;;
;; This is our desired result:
;;
;; s = 1/(2*m)*cn*dn*[u - elliptic_e(asin(sn(u)),m)/(1-m)] + sn*cn^2/2/(1-m)
;;
;;
;; Since diff(cn(u,m),m) = p = -s*sn/cn, we have
;;
;; p = -1/(2*m)*sn*dn[u - elliptic_e(asin(sn(u)),m)/(1-m)] - sn^2*cn/2/(1-m)
;;
;; diff(dn(u,m),m) = q = -m*s*sn/dn - sn^2/dn/2
;;
;; q = -1/2*sn*cn*[u-elliptic_e(asin(sn),m)/(1-m)] - m*sn^2*cn^2/dn/2/(1-m)
;;
;;      - sn^2/dn/2
;;
;;   = -1/2*sn*cn*[u-elliptic_e(asin(sn),m)/(1-m)] + dn*sn^2/2/(m-1)
;;
(defprop %jacobi_sn
    ((u m)
     ((mtimes) ((%jacobi_cn) u m) ((%jacobi_dn) u m))
     ((mplus simp)
      ((mtimes simp) ((rat simp) 1 2)
       ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
       ((mexpt simp) ((%jacobi_cn simp) u m) 2) ((%jacobi_sn simp) u m))
      ((mtimes simp) ((rat simp) 1 2) ((mexpt simp) m -1)
       ((%jacobi_cn simp) u m) ((%jacobi_dn simp) u m)
       ((mplus simp) u
	((mtimes simp) -1 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
	 (($elliptic_e simp) ((%asin simp) ((%jacobi_sn simp) u m)) m))))))
  grad)

(defprop %jacobi_cn
    ((u m)
     ((mtimes simp) -1 ((%jacobi_sn simp) u m) ((%jacobi_dn simp) u m))
     ((mplus simp)
      ((mtimes simp) ((rat simp) -1 2)
       ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
       ((%jacobi_cn simp) u m) ((mexpt simp) ((%jacobi_sn simp) u m) 2))
      ((mtimes simp) ((rat simp) -1 2) ((mexpt simp) m -1)
       ((%jacobi_dn simp) u m) ((%jacobi_sn simp) u m)
       ((mplus simp) u
	((mtimes simp) -1 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
	 (($elliptic_e simp) ((%asin simp) ((%jacobi_sn simp) u m)) m))))))
  grad)

(defprop %jacobi_dn
    ((u m)
     ((mtimes) -1 m ((%jacobi_sn) u m) ((%jacobi_cn) u m))
     ((mplus simp)
      ((mtimes simp) ((rat simp) -1 2)
       ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
       ((%jacobi_dn simp) u m) ((mexpt simp) ((%jacobi_sn simp) u m) 2))
      ((mtimes simp) ((rat simp) -1 2) ((%jacobi_cn simp) u m)
       ((%jacobi_sn simp) u m)
       ((mplus simp) u
	((mtimes simp) -1
	 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
	 (($elliptic_e simp) ((%asin simp) ((%jacobi_sn simp) u m)) m))))))
  grad)

;; The inverse elliptic functions.
;;
;; F(phi|m) = asn(sin(phi),m)
;; 
;; so asn(u,m) = F(asin(u)|m)
(defprop %inverse_jacobi_sn
    ((x m)
     ;; 1/sqrt(1-x^2)/sqrt(1-m*x^2)
     ((mtimes simp)
      ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
       ((rat simp) -1 2))
      ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) x 2)))
       ((rat simp) -1 2)))
     ;; diff(F(asin(u)|m),m)
     ((mtimes simp) ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
      ((mplus simp)
       ((mtimes simp) -1 x
	((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
	 ((rat simp) 1 2))
	((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) x 2)))
	 ((rat simp) -1 2)))
       ((mtimes simp) ((mexpt simp) m -1)
	((mplus simp) ((%elliptic_e simp) ((%asin simp) x) m)
	 ((mtimes simp) -1 ((mplus simp) 1 ((mtimes simp) -1 m))
	  ((%elliptic_f simp) ((%asin simp) x) m)))))))
  grad)

;; Let u = inverse_jacobi_cn(x).  Then jacobi_cn(u) = x or
;; sqrt(1-jacobi_sn(u)^2) = x.  Or
;;
;; jacobi_sn(u) = sqrt(1-x^2)
;;
;; So u = inverse_jacobi_sn(sqrt(1-x^2),m) = inverse_jacob_cn(x,m)
;;
(defprop %inverse_jacobi_cn
    ((x m)
     ;; -1/sqrt(1-x^2)/sqrt(1-m+m*x^2)
     ((mtimes simp) -1
      ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
       ((rat simp) -1 2))
      ((mexpt simp)
       ((mplus simp) 1 ((mtimes simp) -1 m)
	             ((mtimes simp) m ((mexpt simp) x 2)))
       ((rat simp) -1 2)))
     ((mtimes simp) ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
      ((mplus simp)
       ((mtimes simp) -1
	((mexpt simp)
	 ((mplus simp) 1
	  ((mtimes simp) -1 m ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))))
	 ((rat simp) -1 2))
	((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2))) ((rat simp) 1 2))
	((mabs simp) x))
       ((mtimes simp) ((mexpt simp) m -1)
	((mplus simp)
	 ((%elliptic_e simp)
	  ((%asin simp)
	   ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2))) ((rat simp) 1 2)))
	  m)
	 ((mtimes simp) -1 ((mplus simp) 1 ((mtimes simp) -1 m))
	  ((%elliptic_f simp)
	   ((%asin simp)
	    ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2))) ((rat simp) 1 2)))
	   m)))))))
  grad)

;; Let u = inverse_jacobi_dn(x).  Then
;;
;; jacobi_dn(u) = x or
;;
;; x^2 = jacobi_dn(u)^2 = 1 - m*jacobi_sn(u)^2
;;
;; so jacobi_sn(u) = sqrt(1-x^2)/sqrt(m)
;;
;; or u = inverse_jacobi_sn(sqrt(1-x^2)/sqrt(m))
(defprop %inverse_jacobi_dn
    ((x m)
     ;; -1/sqrt(1-x^2)/sqrt(x^2+m-1)
     ((mtimes simp)
      ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
       ((rat simp) -1 2))
      ((mexpt simp) ((mplus simp) -1 m ((mexpt simp) x 2)) ((rat simp) -1 2)))
     ((mplus simp)
      ((mtimes simp) ((rat simp) -1 2) ((mexpt simp) m ((rat simp) -3 2))
       ((mexpt simp)
	((mplus simp) 1
	 ((mtimes simp) -1 ((mexpt simp) m -1)
	  ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))))
	((rat simp) -1 2))
       ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
	((rat simp) 1 2))
       ((mexpt simp) ((mabs simp) x) -1))
      ((mtimes simp) ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
       ((mplus simp)
	((mtimes simp) -1 ((mexpt simp) m ((rat simp) -1 2))
	 ((mexpt simp)
	  ((mplus simp) 1
	   ((mtimes simp) -1 ((mexpt simp) m -1)
	    ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))))
	  ((rat simp) 1 2))
	 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
	  ((rat simp) 1 2))
	 ((mexpt simp) ((mabs simp) x) -1))
	((mtimes simp) ((mexpt simp) m -1)
	 ((mplus simp)
	  ((%elliptic_e simp)
	   ((%asin simp)
	    ((mtimes simp) ((mexpt simp) m ((rat simp) -1 2))
	     ((mexpt simp) ((mplus simp) 1
			    ((mtimes simp) -1 ((mexpt simp) x 2)))
	      ((rat simp) 1 2))))
	   m)
	  ((mtimes simp) -1 ((mplus simp) 1 ((mtimes simp) -1 m))
	   ((%elliptic_f simp)
	    ((%asin simp)
	     ((mtimes simp) ((mexpt simp) m ((rat simp) -1 2))
	      ((mexpt simp) ((mplus simp) 1
			     ((mtimes simp) -1 ((mexpt simp) x 2)))
	       ((rat simp) 1 2))))
	    m))))))))
  grad)

;; Define the actual functions for the user
(defmfun $jacobi_sn (u m)
  (simplify (list '(%jacobi_sn) (resimplify u) (resimplify m))))
(defmfun $jacobi_cn (u m)
  (simplify (list '(%jacobi_cn) (resimplify u) (resimplify m))))
(defmfun $jacobi_dn (u m)
  (simplify (list '(%jacobi_dn) (resimplify u) (resimplify m))))

(defmfun $inverse_jacobi_sn (u m)
  (simplify (list '(%inverse_jacobi_sn) (resimplify u) (resimplify m))))

(defmfun $inverse_jacobi_cn (u m)
  (simplify (list '(%inverse_jacobi_cn) (resimplify u) (resimplify m))))

(defmfun $inverse_jacobi_dn (u m)
  (simplify (list '(%inverse_jacobi_dn) (resimplify u) (resimplify m))))

;; Possible forms of a complex number:
;;
;; 2.3
;; $%i
;; ((mplus simp) 2.3 ((mtimes simp) 2.3 $%i))
;; ((mplus simp) 2.3 $%i))
;; ((mtimes simp) 2.3 $%i)
;;


;; Is argument u a complex number with real and imagpart satisfying predicate ntypep?
(defun complex-number-p (u &optional (ntypep 'numberp))
  (labels ((a1 (x) (cadr x))
	   (a2 (x) (caddr x))
	   (a3+ (x) (cdddr x))
	   (N (x) (funcall ntypep x))	     ; N
	   (i (x) (and (eq x '$%i) (N 1)))   ; %i
	   (N+i (x) (and (null (a3+ x))	     ; mplus test is precondition
			 (N (a1 x))
			 (or (i (a2 x))
			     (and (mtimesp (a2 x)) (N*i (a2 x))))))
	   (N*i (x) (and (null (a3+ x))	     ; mtimes test is precondition
			 (N (a1 x))
			 (i (a2 x)))))
    (declare (inline a1 a2 a3+ N i N+i N*i))
    (cond ((N u))			     ;2.3
	  ((atom u) (i u))		     ;%i
	  ((mplusp u) (N+i u))		     ;N+%i, N+N*%i
	  ((mtimesp u) (N*i u))		     ;N*%i
	  (t nil))))
	
(defun complexify (x)
  ;; Convert a Lisp number to a maxima number
  (cond ((realp x) x)
	((complexp x) (add (realpart x) (mul '$%i (imagpart x))))
	(t (merror "Complexify called on ~:M" x))))
   
(defun kc-arg (exp m)
  ;; Replace elliptic_kc(m) in the expression with sym.  Check to see
  ;; if the resulting expression is linear in sym and the constant
  ;; term is zero.  If so, return the coefficient of sym, i.e, the
  ;; coefficient of elliptic_kc(m).
  (let* ((sym (gensym))
	 (arg (maxima-substitute sym `((%elliptic_kc) ,m) exp)))
    (if (and (not (equalp arg exp))
	     (linearp arg sym)
	     (zerop1 (coefficient arg sym 0)))
	(coefficient arg sym 1)
	nil)))

(defun kc-arg2 (exp m)
  ;; Replace elliptic_kc(m) in the expression with sym.  Check to see
  ;; if the resulting expression is linear in sym and the constant
  ;; term is zero.  If so, return the coefficient of sym, i.e, the
  ;; coefficient of elliptic_kc(m), and the constant term.  Otherwise,
  ;; return NIL.
  (let* ((sym (gensym))
	 (arg (maxima-substitute sym `((%elliptic_kc) ,m) exp)))
    (if (and (not (equalp arg exp))
	     (linearp arg sym))
	(list (coefficient arg sym 1)
	      (coefficient arg sym 0))
	nil)))

;; Tell maxima how to simplify the functions
;;
;; FORM is list containing the actual expression.  I don't really know
;; what Y and Z contain.  Most of this modeled after SIMP-%SIN.
(defmfun simp-%jacobi_sn (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z))
	coef)
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ;; Numerically evaluate sn
	   (sn (float u) (float m)))
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   ;; For complex values.  Should we really do this?
	   (let ((result (sn (complex ($realpart u) ($imagpart u))
			     (complex ($realpart m) ($imagpart m)))))
	     (complexify result)))
	  ((zerop1 u)
	   ;; A&S 16.5.1
	   0)
	  ((zerop1 m)
	   ;; A&S 16.6.1
	   `((%sin) ,u))
	  ((onep1 m)
	   ;; A&S 16.6.1
	   `((%tanh) ,u))
	  ((and $trigsign (mminusp* u))
	   (neg (cons-exp '%jacobi_sn (neg u) m)))
	  ((and $triginverses
		(listp u)
		(eq (caar u) '%inverse_jacobi_sn)
		(alike1 (third u) m))
	   ;; jacobi_sn(inverse_jacobi_sn(u,m), m) = u
	   (second u))
	  ;; A&S 16.20.1 (Jacobi's Imaginary transformation)
	  ((and $%iargs (multiplep u '$%i))
	   (mul '$%i
		(cons-exp '%jacobi_sc (coeff u '$%i 1) (add 1 (neg m)))))
	  ((setq coef (kc-arg2 u m))
	   ;; sn(m*K+u)
	   ;;
	   ;; A&S 16.8.1
	   (destructuring-bind (lin const)
	       coef
	     (cond ((integerp lin)
		    (ecase (mod lin 4)
		      (0
		       ;; sn(4*m*K + u) = sn(u), sn(0) = 0
		       (if (zerop1 const)
			   0
			   `((%jacobi_sn simp) ,const ,m)))
		      (1
		       ;; sn(4*m*K + K + u) = sn(K+u) = cd(u)
		       ;; sn(K) = 1
		       (if (zerop1 const)
			   1
			   `((%jacobi_cd simp) ,const ,m)))
		      (2
		       ;; sn(4*m*K+2*K + u) = sn(2*K+u) = -sn(u)
		       ;; sn(2*K) = 0
		       (if (zerop1 const)
			   0
			   (neg `((%jacobi_sn simp) ,const ,m))))
		      (3
		       ;; sn(4*m*K+3*K+u) = sn(2*K + K + u) = -sn(K+u) = -cd(u)
		       ;; sn(3*K) = -1
		       (if (zerop1 const)
			   -1
			   (neg `((%jacobi_cd simp) ,const ,m))))))
		   ((and (alike1 lin 1//2)
			 (zerop1 const))
		    ;; A&S 16.5.2
		    ;;
		    ;; sn(1/2*K) = 1/sqrt(1+sqrt(1-m))
		    `((mexpt simp)
		      ((mplus simp) 1
		       ((mexpt simp)
			((mplus simp) 1 ((mtimes simp) -1 ,m))
			((rat simp) 1 2)))
		      ((rat) -1 2)))
		   ((and (alike1 lin 3//2)
			 (zerop1 const))
		    ;; A&S 16.5.2
		    ;;
		    ;; sn(1/2*K + K) = cd(1/2*K,m)
		    (simplifya
		     `((%jacobi_cd) ((mtimes) ((rat) 1 2) ((%elliptic_kc) ,m))
		       ,m)
		     nil))
		   (t
		    (eqtest (list '(%jacobi_sn) u m) form)))))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%jacobi_sn) u m) form)))))

(defmfun simp-%jacobi_cn (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z))
	coef)
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   (cn (float u) (float m)))
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (let ((result (cn (complex ($realpart u) ($imagpart u))
			     (complex ($realpart m) ($imagpart m)))))
	     (complexify result)))
	  ((zerop1 u)
	   ;; A&S 16.5.1
	   1)
	  ((zerop1 m)
	   ;; A&S 16.6.2
	   `((%cos) ,u))
	  ((onep1 m)
	   ;; A&S 16.6.2
	   `((%sech) ,u))
	  ((and $trigsign (mminusp* u))
	   (cons-exp '%jacobi_cn (neg u) m))
	  ((and $triginverses
		(listp u)
		(eq (caar u) '%inverse_jacobi_cn)
		(alike1 (third u) m))
	   ;; jacobi_cn(inverse_jacobi_cn(u,m), m) = u
	   (second u))
	  ;; A&S 16.20.2 (Jacobi's Imaginary transformation)
	  ((and $%iargs (multiplep u '$%i))
	   (cons-exp '%jacobi_nc (coeff u '$%i 1) (add 1 (neg m))))
	  ((setq coef (kc-arg2 u m))
	   ;; cn(m*K+u)
	   ;;
	   ;; A&S 16.8.2
	   (destructuring-bind (lin const)
	       coef
	     (cond ((integerp lin)
		    (ecase (mod lin 4)
		      (0
		       ;; cn(4*m*K + u) = cn(u),
		       ;; cn(0) = 1
		       (if (zerop1 const)
			   1
			   `((%jacobi_cn simp) ,const ,m)))
		      (1
		       ;; cn(4*m*K + K + u) = cn(K+u) = -sqrt(m1)*sd(u)
		       ;; cn(K) = 0
		       (if (zerop1 const)
			   0
			   (neg `((mtimes simp)
				  ((mexpt simp)
				   ((mplus simp) 1 ((mtimes simp) -1 ,m))
				   ((rat simp) 1 2))
				  ((%jacobi_sd simp) ,const ,m)))))
		      (2
		       ;; cn(4*m*K + 2*K + u) = cn(2*K+u) = -cn(u)
		       ;; cn(2*K) = -1
		       (if (zerop1 const)
			   -1
			   (neg `((%jacobi_cn) ,const ,m))))
		      (3
		       ;; cn(4*m*K + 3*K + u) = cn(2*K + K + u) =
		       ;; -cn(K+u) = sqrt(m1)*sd(u)
		       ;;
		       ;; cn(3*K) = 0
		       (if (zerop1 const)
			   0
			   `((mtimes simp)
			     ((mexpt simp)
			      ((mplus simp) 1 ((mtimes simp) -1 ,m))
			      ((rat simp) 1 2))
			     ((%jacobi_sd simp) ,const ,m))))))
		   ((and (alike1 lin 1//2)
			 (zerop1 const))
		    ;; A&S 16.5.2
		    ;; cn(1/2*K) = (1-m)^(1/4)/sqrt(1+sqrt(1-m))
		    `((mtimes simp)
		      ((mexpt simp) ((mplus simp) 1
				     ((mtimes simp) -1 ,m))
		       ((rat simp) 1 4))
		      ((mexpt simp)
		       ((mplus simp) 1
			((mexpt simp)
			 ((mplus simp) 1
			  ((mtimes simp) -1 ,m))
			 ((rat simp) 1 2)))
		       ((rat simp) -1 2))))
		   (t
		    (eqtest (list '(%jacobi_cn) u m) form)))))
	  (t
	   (eqtest (list '(%jacobi_cn) u m) form)))))

(defmfun simp-%jacobi_dn (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z))
	coef)
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   (dn (float u) (float m)))
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (let ((result (dn (complex ($realpart u) ($imagpart u))
			     (complex ($realpart m) ($imagpart m)))))
	     (complexify result)))
	  ((zerop1 u)
	   ;; A&S 16.5.1
	   1)
	  ((zerop1 m)
	   ;; A&S 16.6.3
	   1)
	  ((onep1 m)
	   ;; A&S 16.6.3
	   `(($sech) ,u))
	  ((and $trigsign (mminusp* u))
	   (cons-exp '%jacobi_dn (neg u) m))
	  ((and $triginverses
		(listp u)
		(eq (caar u) '%inverse_jacobi_dn)
		(alike1 (third u) m))
	   ;; jacobi_dn(inverse_jacobi_dn(u,m), m) = u
	   (second u))
	  ((zerop1 ($ratsimp (sub u (power (sub 1 m) 1//2))))
	   ;; A&S 16.5.3
	   ;; dn(sqrt(1-m),m) = K(m)
	   ($elliptic_kc m))
	  ;; A&S 16.20.2 (Jacobi's Imaginary transformation)
	  ((and $%iargs (multiplep u '$%i))
	   (cons-exp '%jacobi_dc (coeff u '$%i 1)
		     (add 1 (neg m))))
	  ((setq coef (kc-arg2 u m))
	   ;; A&S 16.8.3
	   ;;
	   ;; dn(m*K+u) has period 2K
	   ;;
	   (destructuring-bind (lin const)
	       coef
	     (cond ((integerp lin)
		    (ecase (mod lin 2)
		      (0
		       ;; dn(2*m*K + u) = dn(u)
		       ;; dn(0) = 1
		       (if (zerop1 const)
			   1
			   ;; dn(4*m*K+2*K + u) = dn(2*K+u) = dn(u)
			   `((%jacobi_dn) ,const ,m)))
		      (1
		       ;; dn(2*m*K + K + u) = dn(K + u) = sqrt(1-m)*nd(u)
		       ;; dn(K) = sqrt(1-m)
		       (if (zerop1 const)
			   `((mexpt simp)
			     ((mplus simp) 1 ((mtimes simp) -1 ,m))
			     ((rat simp) 1 2))
			   `((mtimes simp)
			     ((mexpt simp)
			      ((mplus simp) 1 ((mtimes simp) -1 ,m))
			      ((rat simp) 1 2))
			     ((%jacobi_nd simp) ,const ,m))))))
		   ((and (alike1 lin 1//2)
			 (zerop1 const))
		    ;; A&S 16.5.2
		    ;; dn(1/2*K) = (1-m)^(1/4)
		    `((mexpt simp)
		      ((mplus simp) 1 ((mtimes simp) -1 ,m))
		      ((rat simp) 1 4)))
		   (t
		    (eqtest (list '(%jacobi_cn) u m) form)))))
	  (t (eqtest (list '(%jacobi_dn) u m) form)))))

;; Should we simplify the inverse elliptic functions into the
;; appropriate incomplete elliptic integral?  I think we should leave
;; it, but perhaps allow some way to do that transformation if
;; desired.

(defmfun simp-%inverse_jacobi_sn (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ;; Numerically evaluate asn
	   ;;
	   ;; asn(x,m) = F(asin(x),m)
	   (elliptic-f (cl:asin (float u)) m))
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (complexify (elliptic-f (cl:asin (complex ($realpart u) ($imagpart u)))
				   (complex ($realpart m) ($imagpart m)))))
	  ((zerop1 u)
	   ;; asn(0,m) = 0
	   0)
	  ((onep1 u)
	   ;; asn(1,m) = elliptic_kc(m)
	   `(($elliptic_kc) ,m))
	  ((and (numberp u) (onep1 (- u)))
	   ;; asn(-1,m) = -elliptic_kc(m)
	   `((mtimes) -1 ((%elliptic_kc) ,m)))
	  ((zerop1 m)
	   ;; asn(x,0) = F(asin(x),0) = asin(x)
	   `((%elliptic_f) ((%asin) ,u) 0))
	  ((onep1 m)
	   ;; asn(x,1) = F(asin(x),1) = log(tan(pi/2+asin(x)/2))
	   `((%elliptic_f) ((%asin) ,u) 1))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%inverse_jacobi_sn) u m) form)))))

(defmfun simp-%inverse_jacobi_cn (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ;; Numerically evaluate acn
	   ;;
	   ;; acn(x,m) = F(acos(x),m)
	   (elliptic-f (cl:acos (float u)) m))
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (complexify (elliptic-f (cl:acos (complex ($realpart u) ($imagpart u)))
				   (complex ($realpart m) ($imagpart m)))))
	  ((zerop1 m)
	   ;; asn(x,0) = F(acos(x),0) = acos(x)
	   `((%elliptic_f) ((%acos) ,u) 0))
	  ((onep1 m)
	   ;; asn(x,1) = F(asin(x),1) = log(tan(pi/2+asin(x)/2))
	   `((%elliptic_f) ((%acos) ,u) 1))
	  ((zerop1 u)
	   `((%elliptic_kc) ,m))
	  ((onep1 u)
	   0)
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%inverse_jacobi_cn) u m) form)))))

;;
;; adn(u) = x.
;; u = dn(x) = sqrt(1-m*sn(x)^2)
;; sn(x)^2 = (1-u^2)/m
;; sn(x) = sqrt((1-u^2)/m)
;; x = asn(sqrt((1-u^2)/m))
;; x = adn(u) = asn(sqrt((1-u^2)/m))

(defmfun simp-%inverse_jacobi_dn (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ;; Numerically evaluate adn
	   (let* ((u (float u))
		  (m (float m))
		  (phi (/ (* (sqrt (- 1 u)) (sqrt (+ 1 u)))
			  (sqrt m))))
	     (elliptic-f (cl:asin phi) m)))
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (let* ((u (complex ($realpart u) ($imagpart u)))
		  (phi (/ (* (sqrt (- 1 u)) (sqrt (+ 1 u)))
			 (sqrt m))))
	     (complexify (elliptic-f (cl:asin phi)
				     (complex ($realpart m) ($imagpart m))))))
	  ((onep1 m)
	   ;; x = dn(u,1) = sech(u).  so u = asech(x)
	   `((%asech) ,u))
	  ((onep1 u)
	   ;; jacobi_dn(0,m) = 1
	   0)
	  ((zerop1 ($ratsimp (sub u (power (sub 1 m) 1//2))))
	   ;; jacobi_dn(K(m),m) = sqrt(1-m) so
	   ;; inverse_jacobi_dn(sqrt(1-m),m) = K(m)
	   ($elliptic_kc m))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%inverse_jacobi_dn) u m) form)))))

;;;; Elliptic integrals

;; Carlson's elliptic integral of the first kind.
;;
;; ***PURPOSE  Compute the incomplete or complete elliptic integral of the
;;             1st kind.  For X, Y, and Z non-negative and at most one of
;;             them zero, RF(X,Y,Z) = Integral from zero to infinity of
;;                                 -1/2     -1/2     -1/2
;;                       (1/2)(t+X)    (t+Y)    (t+Z)    dt.
;;             If X, Y or Z is zero, the integral is complete.
;;
;;          Value of IER assigned by the DRF routine
;;
;;                   Value assigned         Error Message Printed
;;                   IER = 1                MIN(X,Y,Z) .LT. 0.0
;;                       = 2                MIN(X+Y,X+Z,Y+Z) .LT. LOLIM
;;                       = 3                MAX(X,Y,Z) .GT. UPLIM
;;
;;    Special double precision functions via DRF
;;
;;
;;
;;
;;                   Legendre form of ELLIPTIC INTEGRAL of 1st kind
;;
;;                   -----------------------------------------
;;
;;
;;
;;                                              2         2   2
;;                   F(PHI,K) = SIN(PHI) DRF(COS (PHI),1-K SIN (PHI),1)
;;
;;
;;                                   2
;;                   K(K) = DRF(0,1-K ,1)
;;
;;
;;                          PI/2     2   2      -1/2
;;                        = INT  (1-K SIN (PHI) )   D PHI
;;                           0
;;
;;
;;
;;                   Bulirsch form of ELLIPTIC INTEGRAL of 1st kind
;;
;;                   -----------------------------------------
;;
;;
;;                                           2 2    2
;;                   EL1(X,KC) = X DRF(1,1+KC X ,1+X )
;;
;;
;;                   Lemniscate constant A
;;
;;                   -----------------------------------------
;;
;;
;;                        1      4 -1/2
;;                   A = INT (1-S )    DS = DRF(0,1,2) = DRF(0,2,1)
;;                        0
;;
;;
;;
;;     -------------------------------------------------------------------
;;
;; ***REFERENCES  B. C. Carlson and E. M. Notis, Algorithms for incomplete
;;                  elliptic integrals, ACM Transactions on Mathematical
;;                  Software 7, 3 (September 1981), pp. 398-403.
;;                B. C. Carlson, Computing elliptic integrals by
;;                  duplication, Numerische Mathematik 33, (1979),
;;                  pp. 1-16.
;;                B. C. Carlson, Elliptic integrals of the first kind,
;;                  SIAM Journal of Mathematical Analysis 8, (1977),
;;                 pp. 231-242.

(let ((errtol (expt (* 4 flonum-epsilon) 1/6))
      (uplim (/ most-positive-flonum 5))
      (lolim (* #-gcl least-positive-normalized-flonum
		#+gcl least-positive-flonum
		5))
      (c1 (float 1/24))
      (c2 (float 3/44))
      (c3 (float 1/14)))
  (declare (type flonum errtol c1 c2 c3))
  (defun drf (x y z)
    "Compute Carlson's incomplete or complete elliptic integral of the
first kind:

                   INF
                  /
                  [                     1
  RF(x, y, z) =   I    ----------------------------------- dt
                  ]    SQRT(x + t) SQRT(y + t) SQRT(z + t)
                  /
                   0


where x >= 0, y >= 0, z >=0, and at most one of x, y, z is zero.
"
    (declare (type flonum x y z))
    ;; Check validity of input
    (assert (and (>= x 0) (>= y 0) (>= z 0)
		 (plusp (+ x y)) (plusp (+ x z)) (plusp (+ y z))))
    (assert (<= (max x y z) uplim))
    (assert (>= (min (+ x y) (+ x z) (+ y z)) lolim))
    (locally 
	(declare (type (flonum 0.0) x y)
		 (type (flonum (0.0)) z)
		 (optimize (speed 3)))
      (loop
       (let* ((mu (/ (+ x y z) 3))
	      (x-dev (- 2 (/ (+ mu x) mu)))
	      (y-dev (- 2 (/ (+ mu y) mu)))
	      (z-dev (- 2 (/ (+ mu z) mu))))
	 (when (< (max (abs x-dev) (abs y-dev) (abs z-dev)) errtol)
	   (let ((e2 (- (* x-dev y-dev) (* z-dev z-dev)))
		 (e3 (* x-dev y-dev z-dev)))
	     (return (/ (+ 1
			   (* e2 (- (* c1 e2)
				    1/10
				    (* c2 e3)))
			   (* c3 e3))
			(sqrt mu)))))
	 (let* ((x-root (sqrt x))
		(y-root (sqrt y))
		(z-root (sqrt z))
		(lam (+ (* x-root (+ y-root z-root)) (* y-root z-root))))
	   (setf x (* (+ x lam) 1/4))
	   (setf y (* (+ y lam) 1/4))
	   (setf z (* (+ z lam) 1/4))))))))

(let ((errtol (expt (* 4 flonum-epsilon) 1/6))
      (c1 (float 1/24))
      (c2 (float 3/44))
      (c3 (float 1/14)))
  (declare (type flonum errtol c1 c2 c3))
  (defun crf (x y z)
    "Compute Carlson's incomplete or complete elliptic integral of the
first kind:

                   INF
                  /
                  [                     1
  RF(x, y, z) =   I    ----------------------------------- dt
                  ]    SQRT(x + t) SQRT(y + t) SQRT(z + t)
                  /
                   0

  x, y, and z may be complex.
"
    (declare (number x y z))
    (let ((x (coerce x '(complex flonum)))
	  (y (coerce y '(complex flonum)))
	  (z (coerce z '(complex flonum))))
      (declare (type (complex flonum) x y z)
	       (optimize (speed 3)))
      (loop
	 (let* ((mu (/ (+ x y z) 3))
		(x-dev (- 2 (/ (+ mu x) mu)))
		(y-dev (- 2 (/ (+ mu y) mu)))
		(z-dev (- 2 (/ (+ mu z) mu))))
	   (when (< (max (abs x-dev) (abs y-dev) (abs z-dev)) errtol)
	     (let ((e2 (- (* x-dev y-dev) (* z-dev z-dev)))
		   (e3 (* x-dev y-dev z-dev)))
	       (return (/ (+ 1
			     (* e2 (- (* c1 e2)
				      1/10
				      (* c2 e3)))
			     (* c3 e3))
			  (sqrt mu)))))
	   (let* ((x-root (sqrt x))
		  (y-root (sqrt y))
		  (z-root (sqrt z))
		  (lam (+ (* x-root (+ y-root z-root)) (* y-root z-root))))
	     (setf x (* (+ x lam) 1/4))
	     (setf y (* (+ y lam) 1/4))
	     (setf z (* (+ z lam) 1/4))))))))

;; Elliptic integral of the first kind (Legendre's form):
;;
;;
;;      phi
;;     /
;;     [             1
;;     I    ------------------- ds
;;     ]                  2
;;     /    SQRT(1 - m SIN (s))
;;     0

(defun elliptic-f (phi-arg m-arg)
  (cond ((and (realp m-arg) (realp phi-arg))
	 (let ((phi (float phi-arg))
	       (m (float m-arg)))
	   (cond ((> m 1)
		  ;; A&S 17.4.15
		  (/ (elliptic-f (cl:asin (* (sqrt m) (sin phi))) (/ m))))
		 ((< m 0)
		  ;; A&S 17.4.17
		  (let* ((m (- m))
			 (m+1 (+ 1 m))
			 (root (sqrt m+1))
			 (m/m+1 (/ m m+1)))
		    (- (/ (elliptic-f (float (/ pi 2)) m/m+1)
			  root)
		       (/ (elliptic-f (- (float (/ pi 2)) phi) m/m+1)
			  root))))
		 ((= m 0)
		  ;; A&S 17.4.19
		  phi)
		 ((= m 1)
		  ;; A&S 17.4.21
		  (log (cl:tan (+ (/ phi 2) (float (/ pi 2))))))
		 ((minusp phi)
		  (- (elliptic-f (- phi) m)))
		 ((> phi pi)
		  ;; A&S 17.4.3
		  (multiple-value-bind (s phi-rem)
		      (truncate phi (float pi))
		    (+ (* 2 s (elliptic-k m))
		       (elliptic-f phi-rem m))))
		 ((<= phi (/ pi 2))
		  (let ((sin-phi (sin phi))
			(cos-phi (cos phi))
			(k (sqrt m)))
		    (* sin-phi
		       (drf (* cos-phi cos-phi)
			    (* (- 1 (* k sin-phi))
			       (+ 1 (* k sin-phi)))
			    1.0))))
		 ((< phi pi)
		  (+ (* 2 (elliptic-k m))
		     (elliptic-f (- phi (float pi)) m))))))
	(t
	 (let ((phi (coerce phi-arg '(complex flonum)))
	       (m (coerce m-arg '(complex flonum))))
	   (let ((sin-phi (sin phi))
		 (cos-phi (cos phi))
		 (k (sqrt m)))
	     (* sin-phi
		(crf (* cos-phi cos-phi)
		     (* (- 1 (* k sin-phi))
			(+ 1 (* k sin-phi)))
		     1.0)))))))

;; Complete elliptic integral of the first kind
(defun elliptic-k (m)
  (declare (type flonum m))
  (cond ((> m 1)
	 ;; A&S 17.4.15
	 (/ (elliptic-f (cl:asin (sqrt m)) (/ m))))
	((< m 0)
	 ;; A&S 17.4.17
	 (let* ((m (- m))
		(m+1 (+ 1 m))
		(root (sqrt m+1))
		(m/m+1 (/ m m+1)))
	   (- (/ (elliptic-k m/m+1)
		 root)
	      (/ (elliptic-f 0.0 m/m+1)
		 root))))
	((= m 0)
	 ;; A&S 17.4.19
	 (float (/ pi 2)))
	(t
	 (let ((k (sqrt m)))
	   (drf 0.0 (* (- 1 k)
		       (+ 1 k))
		1.0)))))
;;
;; Carlsons' elliptic integral of the second kind.
;;
;;   1.     DRD
;;          Evaluate an INCOMPLETE (or COMPLETE) ELLIPTIC INTEGRAL
;;          of the second kind
;;          Standard FORTRAN function routine
;;          Double precision version
;;          The routine calculates an approximation result to
;;          DRD(X,Y,Z) = Integral from zero to infinity of
;;                              -1/2     -1/2     -3/2
;;                    (3/2)(t+X)    (t+Y)    (t+Z)    dt,
;;          where X and Y are nonnegative, X + Y is positive, and Z is
;;          positive.  If X or Y is zero, the integral is COMPLETE.

;;    -------------------------------------------------------------------
;;
;;
;;   Special double precision functions via DRD and DRF
;;
;;
;;                  Legendre form of ELLIPTIC INTEGRAL of 2nd kind
;;
;;                  -----------------------------------------
;;
;;
;;                                             2         2   2
;;                  E(PHI,K) = SIN(PHI) DRF(COS (PHI),1-K SIN (PHI),1) -
;;
;;                     2      3             2         2   2
;;                  -(K/3) SIN (PHI) DRD(COS (PHI),1-K SIN (PHI),1)
;;
;;
;;                                  2        2            2
;;                  E(K) = DRF(0,1-K ,1) - (K/3) DRD(0,1-K ,1)
;;
;;                         PI/2     2   2      1/2
;;                       = INT  (1-K SIN (PHI) )  D PHI
;;                          0
;;
;;                  Bulirsch form of ELLIPTIC INTEGRAL of 2nd kind
;;
;;                  -----------------------------------------
;;
;;                                               2 2    2
;;                  EL2(X,KC,A,B) = AX DRF(1,1+KC X ,1+X ) +
;;
;;                                              3          2 2    2
;;                                 +(1/3)(B-A) X DRD(1,1+KC X ,1+X )
;;
;;
;;
;;
;;                  Legendre form of alternative ELLIPTIC INTEGRAL
;;                  of 2nd kind
;;
;;                  -----------------------------------------
;;
;;
;;
;;                            Q     2       2   2  -1/2
;;                  D(Q,K) = INT SIN P  (1-K SIN P)     DP
;;                            0
;;
;;
;;
;;                                     3          2     2   2
;;                  D(Q,K) = (1/3) (SIN Q) DRD(COS Q,1-K SIN Q,1)
;;
;;
;;
;;
;;                  Lemniscate constant  B
;;
;;                  -----------------------------------------
;;
;;
;;
;;
;;                       1    2    4 -1/2
;;                  B = INT  S (1-S )    DS
;;                       0
;;
;;
;;                  B = (1/3) DRD (0,2,1)
;;
;;
;;                  Heuman's LAMBDA function
;;
;;                  -----------------------------------------
;;
;;
;;
;;                  (PI/2) LAMBDA0(A,B) =
;;
;;                                    2                2
;;                 = SIN(B) (DRF(0,COS (A),1)-(1/3) SIN (A) *
;;
;;                            2               2         2       2
;;                  *DRD(0,COS (A),1)) DRF(COS (B),1-COS (A) SIN (B),1)
;;
;;                            2       3             2
;;                  -(1/3) COS (A) SIN (B) DRF(0,COS (A),1) *
;;
;;                           2         2       2
;;                   *DRD(COS (B),1-COS (A) SIN (B),1)
;;
;;
;;
;;                  Jacobi ZETA function
;;
;;                  -----------------------------------------
;;
;;                             2                 2       2   2
;;                  Z(B,K) = (K/3) SIN(B) DRF(COS (B),1-K SIN (B),1)
;;
;;
;;                                       2             2
;;                             *DRD(0,1-K ,1)/DRF(0,1-K ,1)
;;
;;                               2       3           2       2   2
;;                            -(K /3) SIN (B) DRD(COS (B),1-K SIN (B),1)
;;
;;

(let ((errtol (expt (/ flonum-epsilon 3) 1/6))
      (c1 (float 3/14))
      (c2 (float 1/6))
      (c3 (float 9/22))
      (c4 (float 3/26)))
  (declare (type flonum errtol c1 c2 c3 c4))
  (defun drd (x y z)
    ;; Check validity of input

    (assert (and (>= x 0) (>= y 0) (>= z 0) (plusp (+ x y))) (x y z))
	    
    (let ((x (float x))
	  (y (float y))
	  (z (float z))
	  (sigma 0.0)
	  (power4 1.0))
      (declare (type (flonum 0.0) x y power4 sigma)
	       (type (flonum (0.0)) z)
	       (optimize (speed 3)))
      (loop
       (let* ((mu (* 1/5 (+ x y (* 3 z))))
	      (x-dev (/ (- mu x) mu))
	      (y-dev (/ (- mu y) mu))
	      (z-dev (/ (- mu z) mu)))
	 (when (< (max (abs x-dev) (abs y-dev) (abs z-dev)) errtol)
	   (let* ((ea (* x-dev y-dev))
		  (eb (* z-dev z-dev))
		  (ec (- ea eb))
		  (ed (- ea (* 6 eb)))
		  (ef (+ ed ec ec))
		  (s1 (* ed (+ (- c1)
			       (* 1/4 c3 ed)
			       (* -3/2 c4 z-dev ef))))
		  (s2 (* z-dev (+ (* c2 ef)
				  (* z-dev (+ (* (- c3) ec)
					      (* z-dev c4 ea)))))))
	     (return (+ (* 3 sigma)
			(/ (* power4 (+ 1 s1 s2))
			   (* mu (sqrt mu)))))))
	 (let* ((x-root (sqrt x))
		(y-root (sqrt y))
		(z-root (sqrt z))
		(lam (+ (* x-root (+ y-root z-root)) (* y-root z-root))))
	   (incf sigma (/ power4 z-root (+ z lam)))
	   (setf power4 (/ power4 4))
	   (setf x (/ (+ x lam) 4))
	   (setf y (/ (+ y lam) 4))
	   (setf z (/ (+ z lam) 4))))))))

;; Elliptic integral of the second kind (Legendre's form):
;;
;;
;;      phi
;;     /
;;     [                  2
;;     I    SQRT(1 - m SIN (s)) ds
;;     ]
;;     /
;;      0

(defun elliptic-e (phi m)
  (declare (type flonum phi m))
  (cond ((= m 0)
	 ;; A&S 17.4.23
	 phi)
	((= m 1)
	 ;; A&S 17.4.25
	 (sin phi))
	(t
	 (let* ((sin-phi (sin phi))
		(cos-phi (cos phi))
		(k (sqrt m))
		(y (* (- 1 (* k sin-phi))
		      (+ 1 (* k sin-phi)))))
	   (- (* sin-phi
		 (drf (* cos-phi cos-phi) y 1.0))
	      (* (/ m 3)
		 (expt sin-phi 3)
		 (drd (* cos-phi cos-phi) y 1.0)))))))

;; Complete version
(defun elliptic-ec (m)
  (declare (type flonum m))
  (cond ((= m 0)
	 ;; A&S 17.4.23
	 (float (/ pi 2)))
	((= m 1)
	 ;; A&S 17.4.25
	 1.0)
	(t
	 (let* ((k (sqrt m))
		(y (* (- 1 k)
		      (+ 1 k))))
	   (- (drf 0.0 y 1.0)
	      (* (/ m 3)
		 (drd 0.0 y 1.0)))))))


;; Define the elliptic integrals for maxima
;;
;; We use the definitions given in A&S 17.2.6 and 17.2.8.  In particular:
;;
;;                 phi
;;                /
;;                [             1
;; F(phi|m)  =    I    ------------------- ds
;;                ]                  2
;;                /    SQRT(1 - m SIN (s))
;;                 0
;;
;; and
;;
;;              phi
;;             /
;;             [                  2
;; E(phi|m) =  I    SQRT(1 - m SIN (s)) ds
;;             ]
;;             /
;;              0
;;
;; That is, we do not use the modular angle, alpha, as the second arg;
;; the parameter m = sin(alpha)^2 is used.
;;


(defprop $elliptic_f simp-$elliptic_f operators)
(defprop $elliptic_e simp-$elliptic_e operators)

;; The derivative of F(phi|m) wrt to phi is easy.  The derivative wrt
;; to m is harder.  Here is a derivation.  Hope I got it right.
;;
;; diff(integrate(1/sqrt(1-m*sin(x)^2),x,0,phi), m);
;;
;; 			   PHI
;; 			  /	       2
;; 			  [	    SIN (x)
;; 			  I    ------------------ dx
;; 			  ]		 2    3/2
;; 			  /    (1 - m SIN (x))
;; 			   0
;;  			  --------------------------
;; 				      2
;;
;; 
;; Now use the following relationship that is easily verified:
;;
;;               2                 2
;;    (1 - m) SIN (x)           COS (x)                 COS(x) SIN(x)
;;  ------------------- = ------------------- - DIFF(-------------------, x)
;;                 2                     2                          2
;;  SQRT(1 - m SIN (x))   SQRT(1 - m SIN (x))         SQRT(1 - m SIN (x))
;;
;;
;; Now integrate this to get:
;;
;; 
;; 	       PHI
;; 	      /		    2
;; 	      [		 SIN (x)
;;    (1 - m) I	   ------------------- dx =
;; 	      ]			 2
;; 	      /	   SQRT(1 - m SIN (x))
;; 	       0

;;
;; 			   PHI
;; 			  /	        2
;; 			  [	     COS (x)
;;  			+ I    ------------------- dx
;; 			  ]		     2
;; 			  /    SQRT(1 - m SIN (x))
;; 			   0
;; 			       COS(PHI) SIN(PHI)
;;  			  -  ---------------------
;; 					   2
;; 			     SQRT(1 - m SIN (PHI))
;;
;; Use the fact that cos(x)^2 = 1 - sin(x)^2 to show that this
;; integral on the RHS is:
;;
;;
;;		  (1 - m) elliptic_F(PHI, m) + elliptic_E(PHI, m)
;; 		  -------------------------------------------
;;				       m
;; So, finally, we have
;;
;;
;; 
;;   d			    
;; 2 -- (elliptic_F(PHI, m)) = 
;;   dm				
;;
;;  elliptic_E(PHI, m) - (1 - m) elliptic_F(PHI, m)     COS(PHI) SIN(PHI)
;;  ---------------------------------------------- - ---------------------
;; 			   m					  2
;; 						     SQRT(1 - m SIN (PHI))
;;   ----------------------------------------------------------------------
;; 				     1 - m

(defprop $elliptic_f
    ((phi m)
     ;; diff wrt phi
     ;; 1/sqrt(1-m*sin(phi)^2)
     ((mexpt simp)
      ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) ((%sin simp) phi) 2)))
      ((rat simp) -1 2))
     ;; diff wrt m
     ((mtimes simp) ((rat simp) 1 2)
      ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
      ((mplus simp)
       ((mtimes simp) ((mexpt simp) m -1)
	((mplus simp) (($elliptic_e simp) phi m)
	 ((mtimes simp) -1 ((mplus simp) 1 ((mtimes simp) -1 m))
	  (($elliptic_f simp) phi m))))
       ((mtimes simp) -1 ((%cos simp) phi) ((%sin simp) phi)
	((mexpt simp)
	 ((mplus simp) 1
	  ((mtimes simp) -1 m ((mexpt simp) ((%sin simp) phi) 2)))
	 ((rat simp) -1 2))))))
  grad)

;;
;; The derivative of E(phi|m) wrt to m is much simpler to derive than F(phi|m).
;;
;; Take the derivative of the definition to get
;;
;; 	    PHI
;; 	   /		 2
;; 	   [	      SIN (x)
;; 	   I    ------------------- dx
;; 	   ]		      2
;; 	   /    SQRT(1 - m SIN (x))
;; 	    0
;; 	 - ---------------------------
;; 			2
;;
;; It is easy to see that
;;
;; 			    PHI
;; 			   /		 2
;; 			   [	      SIN (x)
;;  elliptic_F(PHI, m) - m I    ------------------- dx = elliptic_E(PHI, m)
;; 			   ]		      2
;; 			   /    SQRT(1 - m SIN (x))
;; 			    0
;;
;; So we finally have
;;
;;   d			       elliptic_E(PHI, m) - elliptic_F(PHI, m)
;;   -- (elliptic_E(PHI, m)) = ---------------------------------------
;;   dm					        2 m

(defprop $elliptic_e
    ((phi m)
     ;; sqrt(1-m*sin(phi)^2)
     ((mexpt simp)
      ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) ((%sin simp) phi) 2)))
      ((rat simp) 1 2))
     ;; diff wrt m
     ((mtimes simp) ((rat simp) 1 2) ((mexpt simp) m -1)
      ((mplus simp) (($elliptic_e simp) phi m)
       ((mtimes simp) -1 (($elliptic_f simp) phi m)))))
  grad)
		    
(defmfun $elliptic_f (phi m)
  (simplify (list '($elliptic_f) (resimplify phi) (resimplify m))))
(defmfun $elliptic_e (phi m)
  (simplify (list '($elliptic_e) (resimplify phi) (resimplify m))))

(defmfun simp-$elliptic_f (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((phi (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp phi) (floatp m))
	       (and $numer (numberp phi) (numberp m)))
	   ;; Numerically evaluate it
	   (elliptic-f (float phi) (float m)))
	  ((and $numer (complex-number-p phi)
		(complex-number-p m))
	   (complexify (elliptic-f (complex ($realpart phi) ($imagpart phi))
				   (complex ($realpart m) ($imagpart m)))))
	  ((zerop1 phi)
	   0)
	  ((zerop1 m)
	   ;; A&S 17.4.19
	   phi)
	  ((onep1 m)
	   ;; A&S 17.4.21.  Let's pick the log tan form.
	   `((%log) ((%tan)
		     ((mplus) ((mtimes) $%pi ((rat) 1 4))
		      ((mtimes) ((rat) 1 2) ,phi)))))
	  ((alike1 phi '((mtimes) ((rat) 1 2) $%pi))
	   ;; Complete elliptic integral
	   `((%elliptic_kc) ,m))
	  (t
	   ;; Nothing to do
	   (eqtest (list '($elliptic_f) phi m) form)))))

(defmfun simp-$elliptic_e (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((phi (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp phi) (floatp m))
	       (and $numer (numberp phi) (numberp m)))
	   ;; Numerically evaluate it
	   (elliptic-e (float phi) (float m)))
	  ((zerop1 phi)
	   0)
	  ((zerop1 m)
	   ;; A&S 17.4.23
	   phi)
	  ((onep1 m)
	   ;; A&S 17.4.25
	   `((%sin) ,phi))
	  ((alike1 phi '((mtimes) ((rat) 1 2) $%pi))
	   ;; Complete elliptic integral
	   `((%elliptic_ec) ,m))
	  (t
	   ;; Nothing to do
	   (eqtest (list '($elliptic_e) phi m) form)))))
    

;; Complete elliptic integrals
;;
;; elliptic_kc(m) = elliptic_f(%pi/2, m)
;;
;; elliptic_ec(m) = elliptic_e(%pi/2, m)
;;
(defmfun $elliptic_kc (m)
  (simplify (list '(%elliptic_kc) (resimplify m))))
(defmfun $elliptic_ec (m)
  (simplify (list '(%elliptic_ec) (resimplify m))))


(defprop %elliptic_kc simp-%elliptic_kc operators)
(defprop %elliptic_ec simp-%elliptic_ec operators)

(defmfun simp-%elliptic_kc (form yy z)
  (declare (ignore yy))
  (oneargcheck form)
  (let ((m (simpcheck (cadr form) z)))
    (cond ((or (and (floatp m))
	       (and $numer (numberp m)))
	   ;; Numerically evaluate it
	   (elliptic-k (float m)))
	  ((zerop1 m)
	   '((mtimes) ((rat) 1 2) $%pi))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%elliptic_kc) m) form)))))

(defprop %elliptic_kc
    ((m)
     ;; diff wrt m
     ((mtimes)
      ((rat) 1 2)
      ((mplus) ((%elliptic_ec) m)
       ((mtimes) -1
	((%elliptic_kc) m)
	((mplus) 1 ((mtimes) -1 m))))
      ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
      ((mexpt) m -1)))
  grad)

(defmfun simp-%elliptic_ec (form yy z)
  (declare (ignore yy))
  (oneargcheck form)
  (let ((m (simpcheck (cadr form) z)))
    (cond ((or (and (floatp m))
	       (and $numer (numberp m)))
	   ;; Numerically evaluate it
	   (elliptic-ec (float m)))
	  ((zerop1 m)
	   '((mtimes) ((rat) 1 2) $%pi))
	  ;; Some special cases we know about.
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%elliptic_ec) m) form)))))

(defprop %elliptic_ec
    ((m)
     ((mtimes) ((rat) 1 2)
      ((mplus) ((%elliptic_ec) m)
       ((mtimes) -1 ((%elliptic_kc)
		     m)))
      ((mexpt) m -1)))
  grad)

;;
;; Elliptic integral of the third kind:
;;
;; (A&S 17.2.14)
;;
;;                 phi
;;                /
;;                [                     1
;; PI(n;phi|m) =  I    ----------------------------------- ds
;;                ]                  2               2
;;                /    SQRT(1 - m SIN (s)) (1 - n SIN (s))
;;                 0
;;
;; As with E and F, we do not use the modular angle alpha but the
;; parameter m = sin(alpha)^2.
;;
(defprop $elliptic_pi simp-$elliptic_pi operators)

(defmfun $elliptic_pi (n phi m)
  (simplify (list '($elliptic_pi)
		  (resimplify n) (resimplify phi) (resimplify m))))

(defmfun simp-$elliptic_pi (form yy z)
  (declare (ignore yy))
  ;;(threeargcheck form)
  (let ((n (simpcheck (cadr form) z))
	(phi (simpcheck (caddr form) z))
	(m (simpcheck (cadddr form) z)))
    (cond ((or (and (floatp n) (floatp phi) (floatp m))
	       (and $numer (numberp n) (numberp phi) (numberp m)))
	   ;; Numerically evaluate it
	   (elliptic-pi (float n) (float phi) (float m)))
	  ((zerop1 n)
	   `(($elliptic_f) ,phi ,m))
	  ((zerop1 m)
	   ;; 3 cases depending on n < 1, n > 1, or n = 1.
	   (let ((s (asksign `((mplus) -1 ,n))))
	     (case s
	       ($positive
		`((mtimes)
		  ((mexpt) ((mplus) -1 ,n) ((rat) -1 2))
		  ((%atanh)
		   ((mtimes) ((mexpt) ((mplus) -1 ,n) ((rat) 1 2))
		    ((%tan) ,phi)))))
	       ($negative
		`((mtimes)
		  ((mexpt) ((mplus) 1 ((mtimes) -1 ,n)) ((rat) -1 2))
		  ((%atan)
		   ((mtimes) ((mexpt)
			      ((mplus) 1 ((mtimes) -1 ,n))
			      ((rat) 1 2))
		    ((%tan) ,phi)))))
	       ($zero
		`((%tan) ,phi)))))
	  (t
	   ;; Nothing to do
	   (eqtest (list '($elliptic_pi) n phi m) form)))))

(defun elliptic-pi (n phi m)
  ;; Note: Carlson's DRJ has n defined as the negative of the n given
  ;; in A&S.
  (let* ((nn (- n))
	 (sin-phi (sin phi))
	 (cos-phi (cos phi))
	 (k (sqrt m))
	 (k2sin (* (- 1 (* k sin-phi))
		   (+ 1 (* k sin-phi)))))
    (- (* sin-phi (drf (expt cos-phi 2) k2sin 1.0))
       (* (/ nn 3) (expt sin-phi 3)
	  (drj (expt cos-phi 2) k2sin 1.0
	       (- 1 (* n (expt sin-phi 2))))))))
    
;;***PURPOSE  Calculate a double precision approximation to
;;             DRC(X,Y) = Integral from zero to infinity of
;;                              -1/2     -1
;;                    (1/2)(t+X)    (t+Y)  dt,
;;            where X is nonnegative and Y is positive.
;;   1.     DRC
;;          Standard FORTRAN function routine
;;          Double precision version
;;          The routine calculates an approximation result to
;;          DRC(X,Y) = integral from zero to infinity of
;;
;;                              -1/2     -1
;;                    (1/2)(t+X)    (t+Y)  dt,
;;
;;          where X is nonnegative and Y is positive.  The duplication
;;          theorem is iterated until the variables are nearly equal,
;;          and the function is then expanded in Taylor series to fifth
;;          order.  Logarithmic, inverse circular, and inverse hyper-
;;          bolic functions can be expressed in terms of DRC.
;;
;;   --------------------------------------------------------------------
;;
;;   Special functions via DRC
;;
;;
;;
;;                  LN X                X .GT. 0
;;
;;                                             2
;;                  LN(X) = (X-1) DRC(((1+X)/2)  , X )
;;
;;
;;   --------------------------------------------------------------------
;;
;;                  ARCSIN X            -1 .LE. X .LE. 1
;;
;;                                       2
;;                  ARCSIN X = X DRC (1-X  ,1 )
;;
;;   --------------------------------------------------------------------
;;
;;                  ARCCOS X            0 .LE. X .LE. 1
;;
;;
;;                                     2       2
;;                  ARCCOS X = SQRT(1-X ) DRC(X  ,1 )
;;
;;   --------------------------------------------------------------------
;;
;;                  ARCTAN X            -INF .LT. X .LT. +INF
;;
;;                                        2
;;                  ARCTAN X = X DRC(1,1+X  )
;;
;;   --------------------------------------------------------------------
;;
;;                  ARCCOT X            0 .LE. X .LT. INF
;;
;;                                  2   2
;;                  ARCCOT X = DRC(X  ,X +1 )
;;
;;   --------------------------------------------------------------------
;;
;;                  ARCSINH X           -INF .LT. X .LT. +INF
;;
;;                                       2
;;                  ARCSINH X = X DRC(1+X  ,1 )
;;
;;   --------------------------------------------------------------------
;;
;;                  ARCCOSH X           X .GE. 1
;;
;;                                    2         2
;;                  ARCCOSH X = SQRT(X -1) DRC(X  ,1 )
;;
;;   --------------------------------------------------------------------
;;
;;                  ARCTANH X           -1 .LT. X .LT. 1
;;
;;                                         2
;;                  ARCTANH X = X DRC(1,1-X  )
;;
;;   --------------------------------------------------------------------
;;
;;                  ARCCOTH X           X .GT. 1
;;
;;                                   2   2
;;                  ARCCOTH X = DRC(X  ,X -1 )
;;
;;   --------------------------------------------------------------------
;;
;;***REFERENCES  B. C. Carlson and E. M. Notis, Algorithms for incomplete
;;                 elliptic integrals, ACM Transactions on Mathematical
;;                 Software 7, 3 (September 1981), pp. 398-403.
;;               B. C. Carlson, Computing elliptic integrals by
;;                 duplication, Numerische Mathematik 33, (1979),
;;                 pp. 1-16.
;;               B. C. Carlson, Elliptic integrals of the first kind,
;;                 SIAM Journal of Mathematical Analysis 8, (1977),
;;                 pp. 231-242.


(let ((errtol (expt (/ flonum-epsilon 3) 1/6))
      (c1 (float 1/7))
      (c2 (float 9/22)))
  (declare (type flonum errtol c1 c2))
  (defun drc (x y)
    (declare (type (flonum 0.0) x)
	     (type (flonum (0.0)) y)
	     (optimize (speed 3)))
    (let ((xn x)
	  (yn y))
      (declare (type (flonum 0.0) xn)
	       (type (flonum (0.0)) yn))
      (loop
       (let* ((mu (/ (+ xn yn yn) 3))
	      (sn (- (/ (+ yn mu) mu) 2)))
	 (declare (type flonum sn))
	 (when (< (abs sn) errtol)
	   (let ((s (* sn sn
		       (+ 0.3
			  (* sn (+ c1 (* sn (+ 0.375 (* sn c2)))))))))
	     (return-from drc (/ (+ 1 s) (sqrt mu)))))
	 (let ((lam (+ (* 2 (sqrt xn) (sqrt yn)) yn)))
	   (setf xn (* (+ xn lam) 0.25))
	   (setf yn (* (+ yn lam) 0.25))))))))

;;   1.     DRJ
;;          Standard FORTRAN function routine
;;          Double precision version
;;          The routine calculates an approximation result to
;;          DRJ(X,Y,Z,P) = Integral from zero to infinity of
;;
;;                                -1/2     -1/2     -1/2     -1
;;                      (3/2)(t+X)    (t+Y)    (t+Z)    (t+P)  dt,
;;
;;          where X, Y, and Z are nonnegative, at most one of them is
;;          zero, and P is positive.  If X or Y or Z is zero, the
;;          integral is COMPLETE.  The duplication theorem is iterated
;;          until the variables are nearly equal, and the function is
;;          then expanded in Taylor series to fifth order.
;;
;;
;;    -------------------------------------------------------------------
;;
;;
;;   Special double precision functions via DRJ and DRF
;;
;;
;;                  Legendre form of ELLIPTIC INTEGRAL of 3rd kind
;;                  -----------------------------------------
;;
;;
;;                          PHI         2         -1
;;             P(PHI,K,N) = INT (1+N SIN (THETA) )   *
;;                           0
;;
;;
;;                                  2    2         -1/2
;;                             *(1-K  SIN (THETA) )     D THETA
;;
;;
;;                                           2          2   2
;;                        = SIN (PHI) DRF(COS (PHI), 1-K SIN (PHI),1)
;;
;;                                   3             2         2   2
;;                         -(N/3) SIN (PHI) DRJ(COS (PHI),1-K SIN (PHI),
;;
;;                                  2
;;                         1,1+N SIN (PHI))
;;
;;
;;
;;                  Bulirsch form of ELLIPTIC INTEGRAL of 3rd kind
;;                  -----------------------------------------
;;
;;
;;                                            2 2    2
;;                  EL3(X,KC,P) = X DRF(1,1+KC X ,1+X ) +
;;
;;                                            3           2 2    2     2
;;                               +(1/3)(1-P) X  DRJ(1,1+KC X ,1+X ,1+PX )
;;
;;
;;                                           2
;;                  CEL(KC,P,A,B) = A RF(0,KC ,1) +
;;
;;
;;                                                      2
;;                                 +(1/3)(B-PA) DRJ(0,KC ,1,P)
;;
;;
;;                  Heuman's LAMBDA function
;;                  -----------------------------------------
;;
;;
;;                                2                      2      2    1/2
;;                  L(A,B,P) =(COS (A)SIN(B)COS(B)/(1-COS (A)SIN (B))   )
;;
;;                                            2         2       2
;;                            *(SIN(P) DRF(COS (P),1-SIN (A) SIN (P),1)
;;
;;                                 2       3            2       2
;;                            +(SIN (A) SIN (P)/(3(1-COS (A) SIN (B))))
;;
;;                                    2         2       2
;;                            *DRJ(COS (P),1-SIN (A) SIN (P),1,1-
;;
;;                                2       2          2       2
;;                            -SIN (A) SIN (P)/(1-COS (A) SIN (B))))
;;
;;
;;
;;                  (PI/2) LAMBDA0(A,B) =L(A,B,PI/2) =
;;
;;                   2                         2       2    -1/2
;;              = COS (A)  SIN(B) COS(B) (1-COS (A) SIN (B))
;;
;;                           2                  2       2
;;                 *DRF(0,COS (A),1) + (1/3) SIN (A) COS (A)
;;
;;                                      2       2    -3/2
;;                 *SIN(B) COS(B) (1-COS (A) SIN (B))
;;
;;                           2         2       2          2       2
;;                 *DRJ(0,COS (A),1,COS (A) COS (B)/(1-COS (A) SIN (B)))
;;
;;
;;                  Jacobi ZETA function
;;                  -----------------------------------------
;;
;;                        2                     2   2    1/2
;;             Z(B,K) = (K/3) SIN(B) COS(B) (1-K SIN (B))
;;
;;
;;                                  2      2   2                 2
;;                        *DRJ(0,1-K ,1,1-K SIN (B)) / DRF (0,1-K ,1)
;;
;;
;;
;;***REFERENCES  B. C. Carlson and E. M. Notis, Algorithms for incomplete
;;                 elliptic integrals, ACM Transactions on Mathematical
;;                 Software 7, 3 (September 1981), pp. 398-403.
;;               B. C. Carlson, Computing elliptic integrals by
;;                 duplication, Numerische Mathematik 33, (1979),
;;                 pp. 1-16.
;;               B. C. Carlson, Elliptic integrals of the first kind,
;;                 SIAM Journal of Mathematical Analysis 8, (1977),
;;                 pp. 231-242.
(let ((errtol (expt (/ flonum-epsilon 3) 1/6))
      (c1 3/14)
      (c2 1/3)
      (c3 3/22)
      (c4 3/26))
  (defun drj (x y z p)
    (declare (type (flonum 0.0) x y z)
	     (type (flonum (0.0)) p))
    (let ((xn x)
	  (yn y)
	  (zn z)
	  (pn p)
	  (sigma 0.0)
	  (power4 1.0))
      (declare (type (flonum 0.0) xn yn zn)
	       (type (flonum (0.0)) pn)
	       (type flonum sigma power4))
      (loop
       (let* ((mu (* 0.2 (+ xn yn zn pn pn)))
	      (xndev (/ (- mu xn) mu))
	      (yndev (/ (- mu yn) mu))
	      (zndev (/ (- mu zn) mu))
	      (pndev (/ (- mu pn) mu))
	      (eps (max (abs xndev) (abs yndev) (abs zndev) (abs pndev))))
	 (when (< eps errtol)
	   (let* ((ea (+ (* xndev (+ yndev zndev))
			 (* yndev zndev)))
		  (eb (* xndev yndev zndev))
		  (ec (* pndev pndev))
		  (e2 (- ea (* 3 ec)))
		  (e3 (+ eb (* 2 pndev (- ea ec))))
		  (s1 (+ 1 (* e2 (+ (- c1)
				    (* 3/4 c3 e2)
				    (* -3/2 c4 e3)))))
		  (s2 (* eb (+ (* 1/2 c2)
			       (* pndev (+ (- c3) (- c3)
					   (* pndev c4))))))
		  (s3 (- (* pndev ea (- c2 (* pndev c3)))
			 (* c2 pndev ec))))
	     (return-from drj (+ (* 3 sigma)
				 (/ (* power4 (+ s1 s2 s3))
				    (* mu (sqrt mu)))))))
	 (let* ((xnroot (sqrt xn))
		(ynroot (sqrt yn))
		(znroot (sqrt zn))
		(lam (+ (* xnroot (+ ynroot znroot))
			(* ynroot znroot)))
		(alfa (expt (+ (* pn (+ xnroot ynroot znroot))
			       (* xnroot ynroot znroot))
			    2))
		(beta (* pn (expt (+ pn lam) 2))))
	   (incf sigma (* power4 (drc alfa beta)))
	   (setf power4 (* power4 0.25))
	   (setf xn (* 0.25 (+ xn lam)))
	   (setf yn (* 0.25 (+ yn lam)))
	   (setf zn (* 0.25 (+ zn lam)))
	   (setf pn (* 0.25 (+ pn lam)))))))))

(defun check-drj (x y z p)
  (let* ((w (/ (* x y) z))
	 (a (* p p (+ x y z w)))
	 (b (* p (+ p x) (+ p y)))
	 (drj-1 (drj x (+ x z) (+ x w) (+ x p)))
	 (drj-2 (drj y (+ y z) (+ y w) (+ y p)))
	 (drj-3 (drj a b b a))
	 (drj-4 (drj 0.0 z w p)))
    ;; Both values should be equal.
    (values (+ drj-1 drj-2 (* (- a b) drj-3) (/ 3 (sqrt a)))
	    drj-4)))

;;; Other Jacobian elliptic functions

;; jacobi_ns(u,m) = 1/jacobi_sn(u,m)
(defmfun  $jacobi_ns (u m)
  (simplify (list '(%jacobi_ns) (resimplify u) (resimplify m))))

(defprop %jacobi_ns simp-%jacobi_ns operators)

(defprop %jacobi_ns
    ((u m)
     ;; diff wrt u
     ((mtimes) -1 ((%jacobi_cn) u m) ((%jacobi_dn) u m)
      ((mexpt) ((%jacobi_sn) u m) -2))
     ;; diff wrt m
     ((mtimes) -1 ((mexpt) ((%jacobi_sn) u m) -2)
      ((mplus)
       ((mtimes) ((rat) 1 2)
	((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	((mexpt) ((%jacobi_cn) u m) 2)
	((%jacobi_sn) u m))
       ((mtimes) ((rat) 1 2) ((mexpt) m -1)
	((%jacobi_cn) u m) ((%jacobi_dn) u m)
	((mplus) u
	 ((mtimes) -1
	  ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	  (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
	   m)))))))
  grad)

(defmfun simp-%jacobi_ns (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z))
	coef)
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ;; Numerically evaluate sn
	   (/ (sn (float u) (float m))))
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (let ((u-r ($realpart u))
		 (u-i ($imagpart u))
		 (m-r ($realpart m))
		 (m-i ($imagpart m)))
	     (complexify (/ (sn (complex u-r u-i) (complex m-r m-i))))))
	  ((zerop1 m)
	   ;; A&S 16.6.10
	   `(($csc) ,u))
	  ((onep1 m)
	   ;; A&S 16.6.10
	   `(($coth) ,u))
	  ((zerop1 u)
	   (dbz-err1 'jacobi_ns))
	  ((and $trigsign (mminusp* u))
	   ;; ns is odd
	   (neg (cons-exp '%jacobi_ns (neg u) m)))
	  ;; A&S 16.20 (Jacobi's Imaginary transformation)
	  ((and $%iargs (multiplep u '$%i))
	   ;; ns(i*u) = 1/sn(i*u) = -i/sc(u,m1) = -i*cs(u,m1)
	   (neg (mul '$%i
		     (cons-exp '%jacobi_cs (coeff u '$%i 1) (add 1 (neg m))))))
	  ((setq coef (kc-arg2 u m))
	   ;; A&S 16.8.10
	   ;;
	   ;; ns(m*K+u) = 1/sn(m*K+u)
	   ;;
	   (destructuring-bind (lin const)
	       coef
	     (cond ((integerp lin)
		    (ecase (mod lin 4)
		      (0
		       ;; ns(4*m*K+u) = ns(u)
		       ;; ns(0) = infinity
		       (if (zerop1 const)
			   (dbz-err1 'jacobi_ns)
			   `((%jacobi_ns simp) ,const ,m)))
		      (1
		       ;; ns(4*m*K + K + u) = ns(K+u) = dc(u)
		       ;; ns(K) = 1
		       (if (zerop1 const)
			   1
			   `((%jacobi_dc simp) ,const ,m)))
		      (2
		       ;; ns(4*m*K+2*K + u) = ns(2*K+u) = -ns(u)
		       ;; ns(2*K) = infinity
		       (if (zerop1 const)
			   (dbz-err1 'jacobi_ns)
			   (neg `((%jacobi_ns simp) ,const ,m))))
		      (3
		       ;; ns(4*m*K+3*K+u) = ns(2*K + K + u) = -ns(K+u) = -dc(u)
		       ;; ns(3*K) = -1
		       (if (zerop1 const)
			   -1
			   (neg `((%jacobi_dc simp) ,const ,m))))))
		   ((and (alike1 lin 1//2)
			 (zerop1 const))
		    `((mexpt) ((%jacobi_sn) ,u ,m) -1))
		   (t
		    (eqtest (list '(%jacobi_ns) u m) form)))))	  
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%jacobi_ns) u m) form)))))

;; jacobi_nc(u,m) = 1/jacobi_cn(u,m)

(defmfun  $jacobi_nc (u m)
  (simplify (list '(%jacobi_nc) (resimplify u) (resimplify m))))

(defprop %jacobi_nc simp-%jacobi_nc operators)

(defprop %jacobi_nc
    ((u m)
     ;; wrt u
     ((mtimes) ((mexpt) ((%jacobi_cn) u m) -2)
      ((%jacobi_dn) u m) ((%jacobi_sn) u m))
     ;; wrt m
     ((mtimes) -1 ((mexpt) ((%jacobi_cn) u m) -2)
      ((mplus)
       ((mtimes) ((rat) -1 2)
	((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	((%jacobi_cn) u m) ((mexpt) ((%jacobi_sn) u m) 2))
       ((mtimes) ((rat) -1 2) ((mexpt) m -1)
	((%jacobi_dn) u m) ((%jacobi_sn) u m)
	((mplus) u
	 ((mtimes) -1 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	  (($elliptic_e) ((%asin) ((%jacobi_sn) u m)) m)))))))
  grad)

(defmfun simp-%jacobi_nc (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z))
	coef)
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   (/ (cn (float u) (float m))))
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (let ((u-r ($realpart u))
		 (u-i ($imagpart u))
		 (m-r ($realpart m))
		 (m-i ($imagpart m)))
	     (complexify (/ (cn (complex u-r u-i) (complex m-r m-i))))))
	  ((zerop1 u)
	   1)
	  ((zerop1 m)
	   ;; A&S 16.6.8
	   `(($sec) ,u))
	  ((onep1 m)
	   ;; A&S 16.6.8
	   `((%cosh) ,u))
	  ((and $trigsign (mminusp* u))
	   ;; nc is even
	   (cons-exp '%jacobi_nc (neg u) m))
	  ;; A&S 16.20 (Jacobi's Imaginary transformation)
	  ((and $%iargs (multiplep u '$%i))
	   ;; nc(i*u) = 1/cn(i*u) = 1/nc(u,1-m) = cn(u,1-m)
	   (cons-exp '%jacobi_cn (coeff u '$%i 1) (add 1 (neg m))))
	  ((setq coef (kc-arg2 u m))
	   ;; A&S 16.8.8
	   ;;
	   ;; nc(u) = 1/cn(u)
	   ;;
	   (destructuring-bind (lin const)
	       coef
	     (cond ((integerp lin)
		    (ecase (mod lin 4)
		      (0
		       ;; nc(4*m*K+u) = nc(u)
		       ;; nc(0) = 1
		       (if (zerop1 const)
			   1
			   `((%jacobi_nc simp) ,const ,m)))
		      (1
		       ;; nc(4*m*K+K+u) = nc(K+u) = -ds(u)/sqrt(1-m)
		       ;; nc(K) = infinity
		       (if (zerop1 const)
			   (dbz-err1 'jacobi_nc)
			   (neg `((mtimes simp)
				  ((mexpt simp)
				   ((mplus simp) 1 ((mtimes simp) -1 ,m))
				   ((rat simp) -1 2))
				  ((%jacobi_ds simp) ,const ,m)))))
		      (2
		       ;; nc(4*m*K+2*K+u) = nc(2*K+u) = -nc(u)
		       ;; nc(2*K) = -1
		       (if (zerop1 const)
			   -1
			   (neg `((%jacobi_nc) ,const ,m))))
		      (3
		       ;; nc(4*m*K+3*K+u) = nc(3*K+u) = nc(2*K+K+u) =
		       ;; -nc(K+u) = ds(u)/sqrt(1-m)
		       ;;
		       ;; nc(3*K) = infinity
		       (if (zerop1 const)
			   (dbz-err1 'jacobi_nc)
			   `((mtimes simp)
			     ((mexpt simp)
			      ((mplus simp) 1 ((mtimes simp) -1 ,m))
			      ((rat simp) -1 2))
			     ((%jacobi_ds simp) ,const ,m))))))
		   ((and (alike1 1//2 lin)
			 (zerop1 const))
		    `((mexpt) ((%jacobi_cn) ,u ,m) -1))
		   (t
		    (eqtest (list '(%jacobi_cn) u m) form)))))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%jacobi_nc) u m) form)))))

;; jacobi_nd(u,m) = 1/jacobi_dn(u,m)
(defmfun  $jacobi_nd (u m)
  (simplify (list '(%jacobi_nd) (resimplify u) (resimplify m))))

(defprop %jacobi_nd simp-%jacobi_nd operators)

(defprop %jacobi_nd
    ((u m)
     ;; wrt u
     ((mtimes) m ((%jacobi_cn) u m)
      ((mexpt) ((%jacobi_dn) u m) -2) ((%jacobi_sn) u m))
     ;; wrt m
     ((mtimes) -1 ((mexpt) ((%jacobi_dn) u m) -2)
      ((mplus)
       ((mtimes) ((rat) -1 2)
	((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	((%jacobi_dn) u m)
	((mexpt) ((%jacobi_sn) u m) 2))
       ((mtimes) ((rat) -1 2) ((%jacobi_cn) u m)
	((%jacobi_sn) u m)
	((mplus) u
	 ((mtimes) -1
	  ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	  (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
	   m)))))))
  grad)

(defmfun simp-%jacobi_nd (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z))
	coef)
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   (/ (dn (float u) (float m))))
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (let ((u-r ($realpart u))
		 (u-i ($imagpart u))
		 (m-r ($realpart m))
		 (m-i ($imagpart m)))
	     (complexify (/ (dn (complex u-r u-i) (complex m-r m-i))))))
	  ((zerop1 u)
	   1)
	  ((zerop1 m)
	   ;; A&S 16.6.6
	   1)
	  ((onep1 m)
	   ;; A&S 16.6.6
	   `((%cosh) ,u))
	  ((and $trigsign (mminusp* u))
	   ;; nd is even
	   (cons-exp '%jacobi_nd (neg u) m))
	  ;; A&S 16.20 (Jacobi's Imaginary transformation)
	  ((and $%iargs (multiplep u '$%i))
	   ;; nd(i*u) = 1/dn(i*u) = 1/dc(u,1-m) = cd(u,1-m)
	   (cons-exp '%jacobi_cd (coeff u '$%i 1) (add 1 (neg m))))
	  ((setq coef (kc-arg2 u m))
	   ;; A&S 16.8.6
	   ;;
	   (destructuring-bind (lin const)
	       coef
	     (cond ((integerp lin)
		    ;; nd has period 2K
		    (ecase (mod lin 2)
		      (0
		       ;; nd(2*m*K+u) = nd(u)
		       ;; nd(0) = 1
		       (if (zerop1 const)
			   1
			   `((%jacobi_nd) ,const ,m)))
		      (1
		       ;; nd(2*m*K+K+u) = nd(K+u) = dn(u)/sqrt(1-m)
		       ;; nd(K) = 1/sqrt(1-m)
		       (if (zerop1 const)
			   `((mexpt simp)
			     ((mplus simp) 1 ((mtimes simp) -1 ,m))
			     ((rat simp) -1 2))
			   `((mtimes simp)
			     ((%jacobi_nd simp) ,const ,m)
			     ((mexpt simp)
			      ((mplus simp) 1 ((mtimes simp) -1 ,m))
			      ((rat simp) -1 2)))))))
		   (t
		    (eqtest (list '(%jacobi_nd) u m) form)))))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%jacobi_nd) u m) form)))))

;; jacobi_sc(u,m) = jacobi_sn/jacobi_cn
(defmfun  $jacobi_sc (u m)
  (simplify (list '(%jacobi_sc) (resimplify u) (resimplify m))))

(defprop %jacobi_sc simp-%jacobi_sc operators)

(defprop %jacobi_sc
    ((u m)
     ;; wrt u
     ((mtimes) ((mexpt) ((%jacobi_cn) u m) -2)
      ((%jacobi_dn) u m))
     ;; wrt m
     ((mplus)
      ((mtimes) ((mexpt) ((%jacobi_cn) u m) -1)
       ((mplus)
	((mtimes) ((rat) 1 2)
	 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	 ((mexpt) ((%jacobi_cn) u m) 2)
	 ((%jacobi_sn) u m))
	((mtimes) ((rat) 1 2) ((mexpt) m -1)
	 ((%jacobi_cn) u m) ((%jacobi_dn) u m)
	 ((mplus) u
	  ((mtimes) -1
	   ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	   (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
	    m))))))
      ((mtimes) -1 ((mexpt) ((%jacobi_cn) u m) -2)
       ((%jacobi_sn) u m)
       ((mplus)
	((mtimes) ((rat) -1 2)
	 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	 ((%jacobi_cn) u m)
	 ((mexpt) ((%jacobi_sn) u m) 2))
	((mtimes) ((rat) -1 2) ((mexpt) m -1)
	 ((%jacobi_dn) u m) ((%jacobi_sn) u m)
	 ((mplus) u
	  ((mtimes) -1
	   ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	   (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
	    m))))))))
  grad)

(defmfun simp-%jacobi_sc (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z))
	coef)
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   (let ((fu (float u))
		 (fm (float m)))
	     (/ (sn fu fm) (cn fu fm))))
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (let ((u-r ($realpart u))
		 (u-i ($imagpart u))
		 (m-r ($realpart m))
		 (m-i ($imagpart m)))
	     (complexify (/ (sn (complex u-r u-i) (complex m-r m-i))
			    (cn (complex u-r u-i) (complex m-r m-i))))))
	  ((zerop1 u)
	   0)
	  ((zerop1 m)
	   ;; A&S 16.6.9
	   `((%tan) ,u))
	  ((onep1 m)
	   ;; A&S 16.6.9
	   `((%sinh) ,u))
	  ((and $trigsign (mminusp* u))
	   ;; sc is odd
	   (neg (cons-exp '%jacobi_sc (neg u) m)))
	  ;; A&S 16.20 (Jacobi's Imaginary transformation)
	  ((and $%iargs (multiplep u '$%i))
	   ;; sc(i*u) = sn(i*u)/cn(i*u) = i*sc(u,m1)/nc(u,m1) = i*sn(u,m1)
	   (mul '$%i
		(cons-exp '%jacobi_sn (coeff u '$%i 1) (add 1 (neg m)))))
	  ((setq coef (kc-arg2 u m))
	   ;; A&S 16.8.9
	   ;; sc(2*m*K+u) = sc(u)
	   (destructuring-bind (lin const)
	       coef
	     (cond ((integerp lin)
		    (ecase (mod lin 2)
		      (0
		       ;; sc(2*m*K+ u) = sc(u)
		       ;; sc(0) = 0
		       (if (zerop1 const)
			   1
			   `((%jacobi_sc simp) ,const ,m)))
		      (1
		       ;; sc(2*m*K + K + u) = sc(K+u)= - cs(u)/sqrt(1-m)
		       ;; sc(K) = infinity
		       (if (zerop1 const)
			   (dbz-err1 'jacobi_sc)
			   `((mtimes simp) -1
			     ((mexpt simp)
			      ((mplus simp) 1 ((mtimes simp) -1 ,m))
			      ((rat simp) -1 2))
			     ((%jacobi_cs simp) ,const ,m))))))
		   ((and (alike1 lin 1//2)
			 (zerop1 const))
		    ;; (1-m)^(1/4)
		    `((mexpt simp)
		      ((mplus simp) 1 ((mtimes simp) -1 ,m))
		      ((rat simp) 1 4)))
		   (t
		    (eqtest (list '(%jacobi_sc) u m) form)))))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%jacobi_sc) u m) form)))))

;; jacobi_sd(u,m) = jacobi_sn/jacobi_dn
(defmfun  $jacobi_sd (u m)
  (simplify (list '(%jacobi_sd) (resimplify u) (resimplify m))))

(defprop %jacobi_sd simp-%jacobi_sd operators)

(defprop %jacobi_sd
    ((u m)
     ;; wrt u
     ((mtimes) ((%jacobi_cn) u m)
      ((mexpt) ((%jacobi_dn) u m) -2))
     ;; wrt m
     ((mplus)
      ((mtimes) ((mexpt) ((%jacobi_dn) u m) -1)
       ((mplus)
	((mtimes) ((rat) 1 2)
	 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	 ((mexpt) ((%jacobi_cn) u m) 2)
	 ((%jacobi_sn) u m))
	((mtimes) ((rat) 1 2) ((mexpt) m -1)
	 ((%jacobi_cn) u m) ((%jacobi_dn) u m)
	 ((mplus) u
	  ((mtimes) -1
	   ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	   (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
	    m))))))
      ((mtimes) -1 ((mexpt) ((%jacobi_dn) u m) -2)
       ((%jacobi_sn) u m)
       ((mplus)
	((mtimes) ((rat) -1 2)
	 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	 ((%jacobi_dn) u m)
	 ((mexpt) ((%jacobi_sn) u m) 2))
	((mtimes) ((rat) -1 2) ((%jacobi_cn) u m)
	 ((%jacobi_sn) u m)
	 ((mplus) u
	  ((mtimes) -1
	   ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	   (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
	    m))))))))
  grad)

(defmfun simp-%jacobi_sd (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z))
	coef)
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   (let ((fu (float u))
		 (fm (float m)))
	     (/ (sn fu fm) (dn fu fm))))
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (let ((u-r ($realpart u))
		 (u-i ($imagpart u))
		 (m-r ($realpart m))
		 (m-i ($imagpart m)))
	     (complexify (/ (sn (complex u-r u-i) (complex m-r m-i))
			    (dn (complex u-r u-i) (complex m-r m-i))))))
	  ((zerop1 u)
	   0)
	  ((zerop1 m)
	   ;; A&S 16.6.5
	   `((%sin) ,u))
	  ((onep1 m)
	   ;; A&S 16.6.5
	   `((%sinh) ,u))
	  ((and $trigsign (mminusp* u))
	   ;; sd is odd
	   (neg (cons-exp '%jacobi_sd (neg u) m)))
	  ;; A&S 16.20 (Jacobi's Imaginary transformation)
	  ((and $%iargs (multiplep u '$%i))
	   ;; sd(i*u) = sn(i*u)/dn(i*u) = i*sc(u,m1)/dc(u,m1) = i*sd(u,m1)
	   (mul '$%i
		(cons-exp '%jacobi_sd (coeff u '$%i 1) (add 1 (neg m)))))
	  ((setq coef (kc-arg2 u m))
	   ;; A&S 16.8.5
	   ;; sd(4*m*K+u) = sd(u)
	   (destructuring-bind (lin const)
	       coef
	     (cond ((integerp lin)
		    (ecase (mod lin 4)
		      (0
		       ;; sd(4*m*K+u) = sd(u)
		       ;; sd(0) = 0
		       (if (zerop1 const)
			   0
			   `((%jacobi_sd simp) ,const ,m)))
		      (1
		       ;; sd(4*m*K+K+u) = sd(K+u) = cn(u)/sqrt(1-m)
		       ;; sd(K) = 1/sqrt(m1)
		       (if (zerop1 const)
			   `((mexpt) ((mplus) 1 ((mtimes) -1 ,m))
			     ((rat) -1 2))
			   `((mtimes simp)
			     ((mexpt simp)
			      ((mplus simp) 1 ((mtimes simp) -1 ,m))
			      ((rat simp) -1 2))
			     ((%jacobi_cn simp) ,const ,m))))
		      (2
		       ;; sd(4*m*K+2*K+u) = sd(2*K+u) = -sd(u)
		       ;; sd(2*K) = 0
		       (if (zerop1 const)
			   0
			   (neg `((%jacobi_sd) ,const ,m))))
		      (3
		       ;; sd(4*m*K+3*K+u) = sd(3*K+u) = sd(2*K+K+u) =
		       ;; -sd(K+u) = -cn(u)/sqrt(1-m)
		       ;; sd(3*K) = -1/sqrt(m1)
		       (if (zerop1 const)
			   (neg `((mexpt)
				  ((mplus simp) 1 ((mtimes simp) -1 ,m))
				  ((rat) -1 2)))
			   (neg `((mtimes simp)
				  ((mexpt simp)
				   ((mplus simp) 1 ((mtimes simp) -1 ,m))
				   ((rat simp) -1 2))
				  ((%jacobi_cn simp) ,const ,m)))))))
		   ((and (alike1 lin 1//2)
			 (zerop1 const))
		    ;; jacobi_sn/jacobi_dn
		    `((mtimes)
		      ((%jacobi_sn) ((mtimes) ((rat) 1 2)
				     ((%elliptic_kc) ,m))
		       ,m)
		      ((mexpt)
		       ((%jacobi_dn) ((mtimes) ((rat) 1 2)
				      ((%elliptic_kc) ,m))
			,m)
		       -1)))
		   (t
		    (eqtest (list '(%jacobi_sd) u m) form)))))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%jacobi_sd) u m) form)))))

;; jacobi_cs(u,m) = jacobi_cn/jacobi_sn
(defmfun  $jacobi_cs (u m)
  (simplify (list '(%jacobi_cs) (resimplify u) (resimplify m))))

(defprop %jacobi_cs simp-%jacobi_cs operators)

(defprop %jacobi_cs
    ((u m)
     ;; wrt u
     ((mtimes) -1 ((%jacobi_dn) u m)
      ((mexpt) ((%jacobi_sn) u m) -2))
     ;; wrt m
     ((mplus)
      ((mtimes) -1 ((%jacobi_cn) u m)
       ((mexpt) ((%jacobi_sn) u m) -2)
       ((mplus)
	((mtimes) ((rat) 1 2)
	 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	 ((mexpt) ((%jacobi_cn) u m) 2)
	 ((%jacobi_sn) u m))
	((mtimes) ((rat) 1 2) ((mexpt) m -1)
	 ((%jacobi_cn) u m) ((%jacobi_dn) u m)
	 ((mplus) u
	  ((mtimes) -1
	   ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	   (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
	    m))))))
      ((mtimes) ((mexpt) ((%jacobi_sn) u m) -1)
       ((mplus)
	((mtimes) ((rat) -1 2)
	 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	 ((%jacobi_cn) u m)
	 ((mexpt) ((%jacobi_sn) u m) 2))
	((mtimes) ((rat) -1 2) ((mexpt) m -1)
	 ((%jacobi_dn) u m) ((%jacobi_sn) u m)
	 ((mplus) u
	  ((mtimes) -1
	   ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	   (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
	    m))))))))
  grad)

(defmfun simp-%jacobi_cs (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z))
	coef)
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   (let ((fu (float u))
		 (fm (float m)))
	     (/ (cn fu fm) (sn fu fm))))
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (let ((u-r ($realpart u))
		 (u-i ($imagpart u))
		 (m-r ($realpart m))
		 (m-i ($imagpart m)))
	     (complexify (/ (cn (complex u-r u-i) (complex m-r m-i))
			    (sn (complex u-r u-i) (complex m-r m-i))))))
	  ((zerop1 m)
	   ;; A&S 16.6.12
	   `(($cot) ,u))
	  ((onep1 m)
	   ;; A&S 16.6.12
	   `(($csch) ,u))
	  ((zerop1 u)
	   (dbz-err1 'jacobi_cs))
	  ((and $trigsign (mminusp* u))
	   ;; cs is odd
	   (neg (cons-exp '%jacobi_cs (neg u) m)))
	  ;; A&S 16.20 (Jacobi's Imaginary transformation)
	  ((and $%iargs (multiplep u '$%i))
	   ;; cs(i*u) = cn(i*u)/sn(i*u) = -i*nc(u,m1)/sc(u,m1) = -i*ns(u,m1)
	   (neg (mul '$%i
		     (cons-exp '%jacobi_ns (coeff u '$%i 1) (add 1 (neg m))))))
	  ((setq coef (kc-arg2 u m))
	   ;; A&S 16.8.12
	   ;; 
	   ;; cs(2*m*K + u) = cs(u)
	   (destructuring-bind (lin const)
	       coef
	     (cond ((integerp lin)
		    (ecase (mod lin 2)
		      (0
		       ;; cs(2*m*K + u) = cs(u)
		       ;; cs(0) = infinity
		       (if (zerop1 const)
			   (dbz-err1 'jacobi_cs)
			   `((%jacobi_cs simp) ,const ,m)))
		      (1
		       ;; cs(K+u) = -sqrt(1-m)*sc(u)
		       ;; cs(K) = 0
		       (if (zerop1 const)
			   0
			   `((mtimes simp) -1
			     ((mexpt simp)
			      ((mplus simp) 1 ((mtimes simp) -1 ,m))
			      ((rat simp) 1 2))
			     ((%jacobi_sc simp) ,const ,m))))))
		   ((and (alike1 lin 1//2)
			 (zerop1 const))
		    ;; 1/jacobi_sc
		    `((mexpt)
		      ((%jacobi_sc) ((mtimes) ((rat) 1 2)
				     ((%elliptic_kc) ,m)) ,m)
		      -1))
		   (t
		    (eqtest (list '(%jacobi_cs simp) u m) form)))))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%jacobi_cs simp) u m) form)))))

;; jacobi_cd(u,m) = jacobi_cn/jacobi_dn
(defmfun  $jacobi_cd (u m)
  (simplify (list '(%jacobi_cd) (resimplify u) (resimplify m))))

(defprop %jacobi_cd simp-%jacobi_cd operators)

(defprop %jacobi_cd
    ((u m)
     ;; wrt u
     ((mtimes) ((mplus) -1 m)
      ((mexpt) ((%jacobi_dn) u m) -2)
      ((%jacobi_sn) u m))
     ;; wrt m
     ((mplus)
      ((mtimes) -1 ((%jacobi_cn) u m)
       ((mexpt) ((%jacobi_dn) u m) -2)
       ((mplus)
	((mtimes) ((rat) -1 2)
	 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	 ((%jacobi_dn) u m)
	 ((mexpt) ((%jacobi_sn) u m) 2))
	((mtimes) ((rat) -1 2) ((%jacobi_cn) u m)
	 ((%jacobi_sn) u m)
	 ((mplus) u
	  ((mtimes) -1
	   ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	   (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
	    m))))))
      ((mtimes) ((mexpt) ((%jacobi_dn) u m) -1)
       ((mplus)
	((mtimes) ((rat) -1 2)
	 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	 ((%jacobi_cn) u m)
	 ((mexpt) ((%jacobi_sn) u m) 2))
	((mtimes) ((rat) -1 2) ((mexpt) m -1)
	 ((%jacobi_dn) u m) ((%jacobi_sn) u m)
	 ((mplus) u
	  ((mtimes) -1
	   ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	   (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
	    m))))))))
  grad)

(defmfun simp-%jacobi_cd (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z))
	coef)
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   (let ((fu (float u))
		 (fm (float m)))
	     (/ (cn fu fm) (dn fu fm))))
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (let ((u-r ($realpart u))
		 (u-i ($imagpart u))
		 (m-r ($realpart m))
		 (m-i ($imagpart m)))
	     (complexify (/ (cn (complex u-r u-i) (complex m-r m-i))
			    (dn (complex u-r u-i) (complex m-r m-i))))))
	  ((zerop1 u)
	   1)
	  ((zerop1 m)
	   ;; A&S 16.6.4
	   `((%cos) ,u))
	  ((onep1 m)
	   ;; A&S 16.6.4
	   1)
	  ((and $trigsign (mminusp* u))
	   ;; cd is even
	   (cons-exp '%jacobi_cd (neg u) m))
	  ;; A&S 16.20 (Jacobi's Imaginary transformation)
	  ((and $%iargs (multiplep u '$%i))
	   ;; cd(i*u) = cn(i*u)/dn(i*u) = nc(u,m1)/dc(u,m1) = nd(u,m1)
	   (cons-exp '%jacobi_nd (coeff u '$%i 1) (add 1 (neg m))))
	  ((setf coef (kc-arg2 u m))
	   ;; A&S 16.8.4
	   ;;
	   (destructuring-bind (lin const)
	       coef
	     (cond ((integerp lin)
		    (ecase (mod lin 4)
		      (0
		       ;; cd(4*m*K + u) = cd(u)
		       ;; cd(0) = 1
		       (if (zerop1 const)
			   1
			   `((%jacobi_cd) ,const ,m)))
		      (1
		       ;; cd(4*m*K + K + u) = cd(K+u) = -sn(u)
		       ;; cd(K) = 0
		       (if (zerop1 const)
			   0
			   (neg `((%jacobi_sn) ,const ,m))))
		      (2
		       ;; cd(4*m*K + 2*K + u) = cd(2*K+u) = -cd(u)
		       ;; cd(2*K) = -1
		       (if (zerop1 const)
			   -1
			   (neg `((%jacobi_cd) ,const ,m))))
		      (3
		       ;; cd(4*m*K + 3*K + u) = cd(2*K + K + u) =
		       ;; -cd(K+u) = sn(u)
		       ;; cd(3*K) = 0
		       (if (zerop1 const)
			   0
			   `((%jacobi_sn) ,const ,m)))))
		   ((and (alike1 lin 1//2)
			 (zerop1 const))
		    ;; jacobi_cn/jacobi_dn
		    `((mtimes)
		      ((%jacobi_cn) ((mtimes) ((rat) 1 2)
				     ((%elliptic_kc) ,m))
		       ,m)
		      ((mexpt)
		       ((%jacobi_dn) ((mtimes) ((rat) 1 2)
				      ((%elliptic_kc) ,m))
			,m)
		       -1)))
		   (t
		    ;; Nothing to do
		    (eqtest (list '(%jacobi_cd) u m) form)))))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%jacobi_cd) u m) form)))))

;; jacobi_ds(u,m) = jacobi_dn/jacobi_sn
(defmfun  $jacobi_ds (u m)
  (simplify (list '(%jacobi_ds) (resimplify u) (resimplify m))))

(defprop %jacobi_ds simp-%jacobi_ds operators)

(defprop %jacobi_ds
    ((u m)
     ;; wrt u
     ((mtimes) -1 ((%jacobi_cn) u m)
      ((mexpt) ((%jacobi_sn) u m) -2))
     ;; wrt m
     ((mplus)
      ((mtimes) -1 ((%jacobi_dn) u m)
       ((mexpt) ((%jacobi_sn) u m) -2)
       ((mplus)
	((mtimes) ((rat) 1 2)
	 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	 ((mexpt) ((%jacobi_cn) u m) 2)
	 ((%jacobi_sn) u m))
	((mtimes) ((rat) 1 2) ((mexpt) m -1)
	 ((%jacobi_cn) u m) ((%jacobi_dn) u m)
	 ((mplus) u
	  ((mtimes) -1
	   ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	   (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
	    m))))))
      ((mtimes) ((mexpt) ((%jacobi_sn) u m) -1)
       ((mplus)
	((mtimes) ((rat) -1 2)
	 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	 ((%jacobi_dn) u m)
	 ((mexpt) ((%jacobi_sn) u m) 2))
	((mtimes) ((rat) -1 2) ((%jacobi_cn) u m)
	 ((%jacobi_sn) u m)
	 ((mplus) u
	  ((mtimes) -1
	   ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	   (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
	    m))))))))
  grad)

(defmfun simp-%jacobi_ds (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z))
	coef)
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   (let ((fu (float u))
		 (fm (float m)))
	     (/ (dn fu fm) (sn fu fm))))
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (let ((u-r ($realpart u))
		 (u-i ($imagpart u))
		 (m-r ($realpart m))
		 (m-i ($imagpart m)))
	     (complexify (/ (dn (complex u-r u-i) (complex m-r m-i))
			    (sn (complex u-r u-i) (complex m-r m-i))))))
	  ((zerop1 m)
	   ;; A&S 16.6.11
	   `(($csc) ,u))
	  ((onep1 m)
	   ;; A&S 16.6.11
	   `(($csch) ,u))
	  ((zerop1 u)
	   (dbz-err1 'jacobi_ds))
	  ((and $trigsign (mminusp* u))
	   (neg (cons-exp '%jacobi_ds (neg u) m)))
	  ;; A&S 16.20 (Jacobi's Imaginary transformation)
	  ((and $%iargs (multiplep u '$%i))
	   ;; ds(i*u) = dn(i*u)/sn(i*u) = -i*dc(u,m1)/sc(u,m1) = -i*ds(u,m1)
	   (neg (mul '$%i
		     (cons-exp '%jacobi_ds (coeff u '$%i 1) (add 1 (neg m))))))
	  ((setf coef (kc-arg2 u m))
	   ;; A&S 16.8.11
	   (destructuring-bind (lin const)
	       coef
	     (cond ((integerp lin)
		    (ecase (mod lin 4)
		      (0
		       ;; ds(4*m*K + u) = ds(u)
		       ;; ds(0) = infinity
		       (if (zerop1 const)
			   (dbz-err1 'jacobi_ds)
			   `((%jacobi_ds) ,const ,m)))
		      (1
		       ;; ds(4*m*K + K + u) = ds(K+u) = sqrt(1-m)*nc(u)
		       ;; ds(K) = sqrt(1-m)
		       (if (zerop1 const)
			   `((mexpt simp)
			     ((mplus simp) 1 ((mtimes simp) -1 ,m))
			     ((rat simp) 1 2))
			   `((mtimes simp)
			     ((mexpt simp)
			      ((mplus simp) 1 ((mtimes simp) -1 ,m))
			      ((rat simp) 1 2))
			     ((%jacobi_nc simp) ,const ,m))))
		      (2
		       ;; ds(4*m*K + 2*K + u) = ds(2*K+u) = -ds(u)
		       ;; ds(0) = pole
		       (if (zerop1 const)
			   (dbz-err1 'jacobi_ds)
			   (neg `((%jacobi_ds) ,const ,m))))
		      (3
		       ;; ds(4*m*K + 3*K + u) = ds(2*K + K + u) =
		       ;; -ds(K+u) = -sqrt(1-m)*nc(u)
		       ;; ds(3*K) = -sqrt(1-m)
		       (if (zerop1 const)
			   (neg `((mexpt simp)
				  ((mplus simp) 1 ((mtimes simp) -1 ,m))
				  ((rat simp) 1 2)))
			   (neg `((mtimes simp)
				  ((mexpt simp)
				   ((mplus simp) 1 ((mtimes simp) -1 ,m))
				   ((rat simp) 1 2))
				  ((%jacobi_nc simp) ,const ,m)))))))
		   ((and (alike1 lin 1//2)
			 (zerop1 const))
		    ;; jacobi_dn/jacobi_sn
		    `((mtimes)
		      ((%jacobi_dn) ((mtimes) ((rat) 1 2)
				     ((%elliptic_kc) ,m))
		       ,m)
		      ((mexpt)
		       ((%jacobi_sn) ((mtimes) ((rat) 1 2)
				      ((%elliptic_kc) ,m))
			,m)
		       -1)))
		   (t
		    ;; Nothing to do
		    (eqtest (list '(%jacobi_ds) u m) form)))))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%jacobi_ds) u m) form)))))

;; jacobi_dc(u,m) = jacobi_dn/jacobi_cn
(defmfun  $jacobi_dc (u m)
  (simplify (list '(%jacobi_dc) (resimplify u) (resimplify m))))

(defprop %jacobi_dc simp-%jacobi_dc operators)

(defprop %jacobi_dc
    ((u m)
     ;; wrt u
     ((mtimes) ((mplus) 1 ((mtimes) -1 m))
      ((mexpt) ((%jacobi_cn) u m) -2)
      ((%jacobi_sn) u m))
     ;; wrt m
     ((mplus)
      ((mtimes) ((mexpt) ((%jacobi_cn) u m) -1)
       ((mplus)
	((mtimes) ((rat) -1 2)
	 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	 ((%jacobi_dn) u m)
	 ((mexpt) ((%jacobi_sn) u m) 2))
	((mtimes) ((rat) -1 2) ((%jacobi_cn) u m)
	 ((%jacobi_sn) u m)
	 ((mplus) u
	  ((mtimes) -1
	   ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	   (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
	    m))))))
      ((mtimes) -1 ((mexpt) ((%jacobi_cn) u m) -2)
       ((%jacobi_dn) u m)
       ((mplus)
	((mtimes) ((rat) -1 2)
	 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	 ((%jacobi_cn) u m)
	 ((mexpt) ((%jacobi_sn) u m) 2))
	((mtimes) ((rat) -1 2) ((mexpt) m -1)
	 ((%jacobi_dn) u m) ((%jacobi_sn) u m)
	 ((mplus) u
	  ((mtimes) -1
	   ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
	   (($elliptic_e) ((%asin) ((%jacobi_sn) u m))
	    m))))))))
  grad)

(defmfun simp-%jacobi_dc (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z))
	coef)
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   (let ((fu (float u))
		 (fm (float m)))
	     (/ (dn fu fm) (cn fu fm))))
	  ((zerop1 u)
	   1)
	  ((zerop1 m)
	   ;; A&S 16.6.7
	   `(($sec) ,u))
	  ((onep1 m)
	   ;; A&S 16.6.7
	   1)
	  ((and $trigsign (mminusp* u))
	   (cons-exp '%jacobi_dc (neg u) m))
	  ;; A&S 16.20 (Jacobi's Imaginary transformation)
	  ((and $%iargs (multiplep u '$%i))
	   ;; dc(i*u) = dn(i*u)/cn(i*u) = dc(u,m1)/nc(u,m1) = dn(u,m1)
	   (cons-exp '%jacobi_dn (coeff u '$%i 1) (add 1 (neg m))))
	  ((setf coef (kc-arg2 u m))
	   ;; See A&S 16.8.7
	   (destructuring-bind (lin const)
	       coef
	     (cond ((integerp lin)
		    (ecase (mod lin 4)
		      (0
		       ;; dc(4*m*K + u) = dc(u)
		       ;; dc(0) = 1
		       (if (zerop1 const)
			   1
			   `((%jacobi_dc) ,const ,m)))
		      (1
		       ;; dc(4*m*K + K + u) = dc(K+u) = -ns(u)
		       ;; dc(K) = pole
		       (if (zerop1 const)
			   (dbz-err1 'jacobi_dc)
			   (neg `((%jacobi_ns simp) ,const ,m))))
		      (2
		       ;; dc(4*m*K + 2*K + u) = dc(2*K+u) = -dc(u)
		       ;; dc(2K) = -1
		       (if (zerop1 const)
			   -1
			   (neg `((%jacobi_dc) ,const ,m))))
		      (3
		       ;; dc(4*m*K + 3*K + u) = dc(2*K + K + u) =
		       ;; -dc(K+u) = ns(u)
		       ;; dc(3*K) = ns(0) = inf
		       (if (zerop1 const)
			   (dbz-err1 'jacobi_dc)
			   `((%jacobi_dc simp) ,const ,m)))))
		   ((and (alike1 lin 1//2)
			 (zerop1 const))
		    ;; jacobi_dn/jacobi_cn
		    `((mtimes)
		      ((%jacobi_dn) ((mtimes) ((rat) 1 2)
				     ((%elliptic_kc) ,m))
		       ,m)
		      ((mexpt)
		       ((%jacobi_cn) ((mtimes) ((rat) 1 2)
				      ((%elliptic_kc) ,m))
			,m)
		       -1)))
		   (t
		    ;; Nothing to do
		    (eqtest (list '(%jacobi_dc) u m) form)))))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%jacobi_dc) u m) form)))))

;;; Other inverse Jacobian functions

;; inverse_jacobi_ns(x)
;;
;; Let u = inverse_jacobi_ns(x).  Then jacobi_ns(u) = x or
;; 1/jacobi_sn(u) = x or
;;
;; jacobi_sn(u) = 1/x
;;
;; so u = inverse_jacobi_sn(1/x)

(defmfun $inverse_jacobi_ns (u m)
  (simplify (list '(%inverse_jacobi_ns) (resimplify u) (resimplify m))))

(defprop %inverse_jacobi_ns
    ((x m)
     ;; -1/sqrt(1-x^2)/sqrt(x^2+m-1)
     ((mtimes) -1
      ((mexpt) ((mplus) -1 ((mexpt) x 2)) ((rat) -1 2))
      ((mexpt)
       ((mplus) ((mtimes simp ratsimp) -1 m) ((mexpt) x 2))
       ((rat) -1 2)))
     ;; wrt m
     ((%derivative) ((%inverse_jacobi_ns) x m) m 1))
  grad)

(defprop %inverse_jacobi_ns simp-%inverse_jacobi_ns operators)

(defmfun simp-%inverse_jacobi_ns (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ;; Numerically evaluate asn
	   ;;
	   ;; ans(x,m) = asn(1/x,m) = F(asin(1/x),m)
	   (elliptic-f (cl:asin (/ (float u))) (float m)))
	  ((and $numer (complex-number-p u)
		(complex-number-p m))
	   (complexify (elliptic-f (cl:asin (/ (complex ($realpart u) ($imagpart u))))
				   (complex ($realpart m) ($imagpart m)))))
	  ((zerop1 m)
	   ;; ans(x,0) = F(asin(1/x),0) = asin(1/x)
	   `((%elliptic_f) ((%asin) ((mexpt) ,u -1)) 0))
	  ((onep1 m)
	   ;; ans(x,1) = F(asin(1/x),1) = log(tan(pi/2+asin(1/x)/2))
	   `((%elliptic_f) ((%asin) ((mexpt) ,u -1)) 1))
	  ((onep1 u)
	   `((%elliptic_kc) ,m))
	  ((alike1 u -1)
	   (neg `((%elliptic_kc) ,m)))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%inverse_jacobi_ns) u m) form)))))

;; inverse_jacobi_nc(x)
;;
;; Let u = inverse_jacobi_nc(x).  Then jacobi_nc(u) = x or
;; 1/jacobi_cn(u) = x or
;;
;; jacobi_cn(u) = 1/x
;;
;; so u = inverse_jacobi_cn(1/x)

(defmfun $inverse_jacobi_nc (u m)
  (simplify (list '(%inverse_jacobi_nc) (resimplify u) (resimplify m))))

(defprop %inverse_jacobi_nc
    ((x m)
     ;; 1/sqrt((x^2-1)*(m+(1-m)*x^2)
     ((mtimes)
      ((mexpt) ((mplus) -1 ((mexpt) x 2)) ((rat) -1 2))
      ((mexpt)
       ((mplus) m
	((mtimes) -1 ((mplus) -1 m) ((mexpt) x 2)))
       ((rat) -1 2)))
     ;; wrt m
     ((%derivative) ((%inverse_jacobi_nc) x m) m 1))
  grad)

(defprop %inverse_jacobi_nc simp-%inverse_jacobi_nc operators)

(defmfun simp-%inverse_jacobi_nc (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ;;
	   ($inverse_jacobi_cn (/ u) m))
	  ((onep1 u)
	   0)
	  ((alike1 u -1)
	   `((mtimes) 2 ((%elliptic_kc) ,m)))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%inverse_jacobi_nc) u m) form)))))

;; inverse_jacobi_nd(x)
;;
;; Let u = inverse_jacobi_nd(x).  Then jacobi_nd(u) = x or
;; 1/jacobi_dn(u) = x or
;;
;; jacobi_dn(u) = 1/x
;;
;; so u = inverse_jacobi_dn(1/x)

(defmfun $inverse_jacobi_nd (u m)
  (simplify (list '(%inverse_jacobi_nd) (resimplify u) (resimplify m))))

(defprop %inverse_jacobi_nd
    ((x m)
     ;; -1/sqrt(1-x^2)/sqrt(x^2+m-1)
     ((mtimes) -1
      ((mexpt) ((mplus) -1 ((mexpt simp ratsimp) x 2))
       ((rat) -1 2))
      ((mexpt)
       ((mplus) 1
	((mtimes) ((mplus) -1 m) ((mexpt simp ratsimp) x 2)))
       ((rat) -1 2)))
     ;; wrt m
     ((%derivative) ((%inverse_jacobi_nd) x m) m 1))
  grad)

(defprop %inverse_jacobi_nd simp-%inverse_jacobi_nd operators)

(defmfun simp-%inverse_jacobi_nd (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ($inverse_jacobi_dn (/ u) m))
	  ((onep1 u)
	   0)
	  ((onep1 ($ratsimp (mul (power (sub 1 m) 1//2) u)))
	   ;; jacobi_nd(1/sqrt(1-m),m) = K(m).  This follows from
	   ;; jacobi_dn(sqrt(1-m),m) = K(m).
	   ($elliptic_kc m))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%inverse_jacobi_nd) u m) form)))))

;; inverse_jacobi_sc(x)
;;
;; Let u = inverse_jacobi_sc(x).  Then jacobi_sc(u) = x or
;; x = jacobi_sn(u)/jacobi_cn(u)
;;
;; x^2 = sn^2/cn^2
;;     = sn^2/(1-sn^2)
;;
;; so
;;
;; sn^2 = x^2/(1+x^2)
;;
;; sn(u) = x/sqrt(1+x^2)
;;
;; u = inverse_sn(x/sqrt(1+x^2))
;;

(defmfun $inverse_jacobi_sc (u m)
  (simplify (list '(%inverse_jacobi_sc) (resimplify u) (resimplify m))))

(defprop %inverse_jacobi_sc
    ((x m)
     ((mtimes)
      ((mexpt) ((mplus) 1 ((mexpt) x 2))
       ((rat) -1 2))
      ((mexpt)
       ((mplus) 1
	((mtimes) -1 ((mplus) -1 m) ((mexpt) x 2)))
       ((rat) -1 2)))
     ;; wrt m
     ((%derivative) ((%inverse_jacobi_sc) x m) m 1))
  grad)

(defprop %inverse_jacobi_sc simp-%inverse_jacobi_sc operators)

(defmfun simp-%inverse_jacobi_sc (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ($inverse_jacobi_sn (/ u (sqrt (+ 1 (* u u)))) m))
	  ((zerop1 u)
	   ;; jacobi_sc(0,m) = 0
	   0)
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%inverse_jacobi_sc) u m) form)))))

;; inverse_jacobi_sd(x)
;;
;; Let u = inverse_jacobi_sd(x).  Then jacobi_sd(u) = x or
;; x = jacobi_sn(u)/jacobi_dn(u)
;;
;; x^2 = sn^2/dn^2
;;     = sn^2/(1-m*sn^2)
;;
;; so
;;
;; sn^2 = x^2/(1+m*x^2)
;;
;; sn(u) = x/sqrt(1+m*x^2)
;;
;; u = inverse_sn(x/sqrt(1+m*x^2))
;;

(defmfun $inverse_jacobi_sd (u m)
  (simplify (list '(%inverse_jacobi_sd) (resimplify u) (resimplify m))))

(defprop %inverse_jacobi_sd
    ((x m)
     ((mtimes)
      ((mexpt)
       ((mplus) 1 ((mtimes) ((mplus) -1 m) ((mexpt) x 2)))
       ((rat) -1 2))
      ((mexpt) ((mplus) 1 ((mtimes) m ((mexpt) x 2)))
       ((rat) -1 2)))
     ;; wrt m
     ((%derivative) ((%inverse_jacobi_sd) x m) m 1))
  grad)

(defprop %inverse_jacobi_sd simp-%inverse_jacobi_sd operators)

(defmfun simp-%inverse_jacobi_sd (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ($inverse_jacobi_sn (/ u (sqrt (+ 1 (* m u u)))) m))
	  ((zerop1 u)
	   0)
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%inverse_jacobi_sd) u m) form)))))

;; inverse_jacobi_cs(x)
;;
;; Let u = inverse_jacobi_cs(x).  Then jacobi_cs(u) = x or
;; 1/x = 1/jacobi_cs(u) = jacobi_sc(u)
;;
;; u = inverse_sc(1/x)
;;

(defmfun $inverse_jacobi_cs (u m)
  (simplify (list '(%inverse_jacobi_cs) (resimplify u) (resimplify m))))

(defprop %inverse_jacobi_cs
    ((x m)
     ((mtimes) -1
      ((mexpt) ((mplus) 1 ((mexpt simp ratsimp) x 2))
       ((rat) -1 2))
      ((mexpt) ((mplus) 1
		((mtimes simp ratsimp) -1 m)
		((mexpt simp ratsimp) x 2))
       ((rat) -1 2)))
     ;; wrt m
     ((%derivative) ((%inverse_jacobi_cs) x m) m 1))
  grad)

(defprop %inverse_jacobi_cs simp-%inverse_jacobi_cs operators)

(defmfun simp-%inverse_jacobi_cs (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ($inverse_jacobi_sc (/ u) m))
	  ((zerop1 u)
	   `((%elliptic_kc) ,m))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%inverse_jacobi_cs) u m) form)))))

;; inverse_jacobi_cd(x)
;;
;; Let u = inverse_jacobi_cd(x).  Then jacobi_cd(u) = x or
;; x = jacobi_cn(u)/jacobi_dn(u)
;;
;; x^2 = cn^2/dn^2
;;     = (1-sn^2)/(1-m*sn^2)
;;
;; or
;;
;; sn^2 = (1-x^2)/(1-m*x^2)
;;
;; sn(u) = sqrt(1-x^2)/sqrt(1-m*x^2)
;;
;; u = inverse_sn(sqrt(1-x^2)/sqrt(1-m*x^2))
;;

(defmfun $inverse_jacobi_cd (u m)
  (simplify (list '(%inverse_jacobi_cd) (resimplify u) (resimplify m))))

(defprop %inverse_jacobi_cd
    ((x m)
     ((mtimes)
      ((mexpt)
       ((mplus) 1 ((mtimes) -1 ((mexpt) x 2)))
       ((rat) -1 2))
      ((mexpt)
       ((mplus) 1 ((mtimes) -1 m ((mexpt) x 2)))
       ((rat) -1 2)))
     ;; wrt m
     ((%derivative) ((%inverse_jacobi_cd) x m) m 1))
  grad)

(defprop %inverse_jacobi_cd simp-%inverse_jacobi_cd operators)

(defmfun simp-%inverse_jacobi_cd (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ($inverse_jacobi_sn (/ (sqrt (* (- 1 u) (+ 1 u)))
				  (sqrt (- 1 (* m u u)))) m))
	  ((onep1 u)
	   0)
	  ((zerop1 u)
	   `((%elliptic_kc) ,m))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%inverse_jacobi_cd) u m) form)))))

;; inverse_jacobi_ds(x)
;;
;; Let u = inverse_jacobi_ds(x).  Then jacobi_ds(u) = x or
;; 1/x = 1/jacobi_ds(u) = jacobi_sd(u)
;;
;; u = inverse_sd(1/x)
;;

(defmfun $inverse_jacobi_ds (u m)
  (simplify (list '(%inverse_jacobi_ds) (resimplify u) (resimplify m))))

(defprop %inverse_jacobi_ds
    ((x m)
     ((mtimes) -1
      ((mexpt)
       ((mplus) -1 m ((mexpt simp ratsimp) x 2))
       ((rat) -1 2))
      ((mexpt)
       ((mplus) m ((mexpt simp ratsimp) x 2))
       ((rat) -1 2)))
     ;; wrt m
     ((%derivative) ((%inverse_jacobi_ds) x m) m 1))
  grad)

(defprop %inverse_jacobi_ds simp-%inverse_jacobi_ds operators)

(defmfun simp-%inverse_jacobi_ds (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ($inverse_jacobi_sd (/ u) m))
	  ((and $trigsign (mminusp* u))
	   (neg (cons-exp '%inverse_jacobi_ds (neg u) m)))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%inverse_jacobi_ds) u m) form)))))


;; inverse_jacobi_dc(x)
;;
;; Let u = inverse_jacobi_dc(x).  Then jacobi_dc(u) = x or
;; 1/x = 1/jacobi_dc(u) = jacobi_cd(u)
;;
;; u = inverse_cd(1/x)
;;

(defmfun $inverse_jacobi_dc (u m)
  (simplify (list '(%inverse_jacobi_dc) (resimplify u) (resimplify m))))

(defprop %inverse_jacobi_dc
    ((x m)
     ((mtimes) -1
      ((mexpt)
       ((mplus) -1 ((mexpt simp ratsimp) x 2))
       ((rat) -1 2))
      ((mexpt)
       ((mplus)
	((mtimes simp ratsimp) -1 m)
	((mexpt simp ratsimp) x 2))
       ((rat) -1 2)))
     ;; wrt m
     ((%derivative) ((%inverse_jacobi_dc) x m) m 1))
  grad)

(defprop %inverse_jacobi_dc simp-%inverse_jacobi_dc operators)

(defmfun simp-%inverse_jacobi_dc (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ($inverse_jacobi_cd (/ u) m))
	  ((onep1 u)
	   0)
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%inverse_jacobi_dc) u m) form)))))

;; Convert an inverse Jacobian function into the equivalent elliptic
;; integral F.
;;
;; See A&S 17.4.41-17.4.52.
(defun make-elliptic-f (e)
  (cond ((atom e)
	 e)
	((member (caar e) '(%inverse_jacobi_sc %inverse_jacobi_cs
			    %inverse_jacobi_nd %inverse_jacobi_dn
			    %inverse_jacobi_sn %inverse_jacobi_cd
			    %inverse_jacobi_dc %inverse_jacobi_ns
			    %inverse_jacobi_nc %inverse_jacobi_ds
			    %inverse_jacobi_sd %inverse_jacobi_cn))
	 ;; We have some inverse Jacobi function.  Convert it to the F form.
	 (destructuring-bind ((fn &rest ops) u m)
	     e
	   (declare (ignore ops))
	   (ecase fn
	     (%inverse_jacobi_sc
	      ;; A&S 17.4.41
	      `(($elliptic_f) ((%atan) ,u) ,m))
	     (%inverse_jacobi_cs
	      ;; A&S 17.4.42
	      `(($elliptic_f) ((%atan) ((mexpt) ,u -1)) ,m))
	     (%inverse_jacobi_nd
	      ;; A&S 17.4.43
	      `(($elliptic_f)
		((%asin) ((mtimes)
			  ((mexpt) ,m ((rat) -1 2))
			  ((mexpt) ,u -1)
			  ((mexpt) ((mplus) -1 ((mexpt) ,u 2))
			   ((rat) 1 2))))
		,m))
	     (%inverse_jacobi_dn
	      ;; A&S 17.4.44
	      `(($elliptic_f)
		((%asin)
		 ((mtimes) ((mexpt) ,m ((rat) -1 2))
		  ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) ,u 2)))
		   ((rat) 1 2))))
		,m))
	     (%inverse_jacobi_sn
	      ;; A&S 17.4.45
	      `(($elliptic_f) ((%asin) ,u) ,m))
	     (%inverse_jacobi_cd
	      ;; A&S 17.4.46
	      `(($elliptic_f)
		((%asin)
		 ((mexpt) ((mtimes) ((mplus) 1
				     ((mtimes) -1 ((mexpt) ,u 2)))
			   ((mexpt) ((mplus) 1
				     ((mtimes) -1 ,m ((mexpt) ,u 2)))
			    -1))
		  ((rat) 1 2)))
		,m))
	     (%inverse_jacobi_dc
	      ;; A&S 17.4.47
	      `(($elliptic_f)
		((%asin)
		 ((mexpt)
		  ((mtimes) ((mplus) -1 ((mexpt) ,u 2))
		   ((mexpt) ((mplus) ((mtimes) -1 ,m) ((mexpt) ,u 2)) -1))
		  ((rat) 1 2)))
		,m))
	     (%inverse_jacobi_ns
	      ;; A&S 17.4.48
	      `(($elliptic_f) ((asin) ((mexpt) ,u -1)) ,m))
	     (%inverse_jacobi_nc
	      ;; A&S 17.4.49
	      `(($elliptic_f) ((acos) ((mexpt) ,u -1)) ,m))
	     (%inverse_jacobi_ds
	      ;; A&S 17.4.50
	      `(($elliptic_f)
		((%asin) ((mexpt) ((mplus) ,m ((mexpt) ,u 2))
			  ((rat) -1 2)))
		,m))
	     (%inverse_jacobi_sd
	      ;; A&S 17.4.51
	      `(($elliptic_f)
		((%asin)
		 ((mtimes) ,u
		  ((mexpt) ((mplus) 1 ((mtimes) ,m ((mexpt) ,u 2)))
		   ((rat) -1 2))))
		,m))
	     (%inverse_jacobi_cn
	      ;; A&S 17.4.52
	      `(($elliptic_f) ((%acos) ,u) ,m)))))
	(t
	 (recur-apply #'make-elliptic-f e))))

(defmfun $make_elliptic_f (e)
  (if (atom e)
      e
      (simplify (make-elliptic-f e))))

(defun make-elliptic-e (e)
  (cond ((atom e) e)
	((eq (caar e) '$elliptic_eu)
	 (destructuring-bind ((ffun &rest ops) u m) e
	   (declare (ignore ffun ops))
	   `(($elliptic_e) ((%asin) ((%jacobi_sn) ,u ,m)) ,m)))
	(t
	 (recur-apply #'make-elliptic-e e))))

(defmfun $make_elliptic_e (e)
  (if (atom e)
      e
      (simplify (make-elliptic-e e))))
  
	 
;; Eu(u,m) = integrate(jacobi_dn(v,m)^2,v,0,u)
;;         = integrate(sqrt((1-m*t^2)/(1-t^2)),t,0,jacobi_sn(u,m))
;;
;; Eu(u,m) = E(am(u),m)
;;
;; where E(u,m) is elliptic-e above.
;;
;; Checks.
;; Lawden gives the following relationships
;;
;; E(u+v) = E(u) + E(v) - m*sn(u)*sn(v)*sn(u+v)
;; E(u,0) = u, E(u,1) = tanh u
;;
;; E(i*u,k) = i*sc(u,k')*dn(u,k') - i*E(u,k') + i*u
;;
;; E(2*i*K') = 2*i*(K'-E')
;;
;; E(u + 2*i*K') = E(u) + 2*i*(K' - E')
;;
;; E(u+K) = E(u) + E - k^2*sn(u)*cd(u)
(defun elliptic-eu (u m)
  (cond ((realp u)
	 ;; E(u + 2*n*K) = E(u) + 2*n*E
	 (let ((ell-k (elliptic-k m))
	       (ell-e (elliptic-ec m)))
	   (multiple-value-bind (n u-rem)
	       (floor u (* 2 ell-k))
	     ;; 0 <= u-rem < 2*K
	     (+ (* 2 n ell-e)
		(cond ((>= u-rem ell-k)
		       ;; 0 <= u-rem < K so
		       ;; E(u + K) = E(u) + E - m*sn(u)*cd(u)
		       (let ((u-k (- u ell-k)))
			 (- (+ (elliptic-e (cl:asin (sn u-k m)) m)
			       ell-e)
			    (/ (* m (sn u-k m) (cn u-k m))
			       (dn u-k m)))))
		      (t
		       (elliptic-e (cl:asin (sn u m)) m)))))))
	((complexp u)
	 ;; From Lawden:
	 ;;
	 ;; E(u+i*v, m) = E(u,m) -i*E(v,m') + i*v + i*sc(v,m')*dn(v,m')
	 ;;                -i*m*sn(u,m)*sc(v,m')*sn(u+i*v,m)
	 ;;
	 (let ((u-r (realpart u))
	       (u-i (imagpart u))
	       (m1 (- 1 m)))
	   (+ (elliptic-eu u-r m)
	      (* #c(0 1)
		 (- (+ u-i
		       (/ (* (sn u-i m1) (dn u-i m1))
			  (cn u-i m1)))
		    (+ (elliptic-eu u-i m1)
		       (/ (* m (sn u-r m) (sn u-i m1) (sn u m))
			  (cn u-i m1))))))))))

(defprop $elliptic_eu simp-$elliptic_eu operators)
(defprop $elliptic_eu
    ((u m)
     ((mexpt) ((%jacobi_dn) u m) 2)
     ;; wrt m
     )
  grad)

(defmfun simp-$elliptic_eu (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   (let ((u-r ($realpart u))
		 (u-i ($imagpart u))
		 (m (float m)))
	     (complexify (elliptic-eu (complex u-r u-i) m))))
	  (t
	   (eqtest `(($elliptic_eu) ,u ,m) form)))))

(defmfun $elliptic_eu (u m)
  (simplify `(($elliptic_eu) ,(resimplify u) ,(resimplify m))))

(defun agm (a0 b0 phi)
  (let (c0)
    (dotimes (k 16)
      (psetq a0 (/ (+ a0 b0) 2)
	     b0 (sqrt (* a0 b0)))
      (setf c0 (/ (- a0 b0) 2))
      (setf phi (+ phi (cl:atan (* (/ b0 a0) (cl:tan phi)))))
      (format t "~A ~A ~A ~A~%" a0 b0 c0 phi))))

(defprop %jacobi_am simp-%jacobi_am operators)

(defmfun $jacobi_am (u m)
  (simplify `((%jacobi_am) ,(resimplify u) ,(resimplify m))))

(defmfun simp-%jacobi_am (form unused z)
  (declare (ignore unused))
  (twoargcheck form)
  (let ((u (simpcheck (cadr form) z))
	(m (simpcheck (caddr form) z)))
    (cond ((or (and (floatp u) (floatp m))
	       (and $numer (numberp u) (numberp m)))
	   ;; Numerically evaluate am
	   (cl:asin (sn (float u) (float m))))
	  (t
	   ;; Nothing to do
	   (eqtest (list '(%jacobi_am) u m) form)))))
  
