;; -*- Mode: Lisp; Package: Climax; Base: 10; Syntax: Common-Lisp; -*-
;;;>
;;;> ** Copyright (c) 1982 - 1994 Macsyma Inc.  All rights reserved. 
;;;> ** Portions copyright (c) 1981 - 1982 Massachusetts Institute of Technology.
;;;> **                        All rights reserved.
;;;
;;; POIS2: MACSYMA Poisson Series package for coefficients in general representation

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; The programs in this file were written by Richard Fateman in 1972-74
;; while a faculty member at MIT.
;; Subsequent to this time, they were slightly modified by others.
;; Significant changes were made in 1982 by Richard Fateman, then
;; at Univ. Calif Berkeley, and the changes are used elsewhere
;; by permission of the author and the Univ. of Calif.
;; additional changes (c) 1986, Richard J. Fateman
;; The following copyright notice was inserted without the author's
;; consent.
;;;     (c) Copyright 1981 Massachusetts Institute of Technology         ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(IN-PACKAGE "CLIMAX")

;;;-------------------------------------------------------------------
;;; To run in Schelter's Maxima comment the above and uncomment these:
;(in-package "maxima")
;(defvar $poisvars '((mlist) $U $v $w $x $y $z))
;(defmacro wrong-number-of-args-error (op) `(wna-err ,op))
;(defmacro improper-arg-error (op arg) `(improper-arg-err ,arg ,op))
;(defmacro mexp-member (item list) `(member ,item ,list :test #'alike1))
;;;; This should replace the poissubstco at the end; 
;(defun poissubstco (uu vv ww) (maxima-substitute uu vv ww))
;;;-------------------------------------------------------------------

;;; NOTE: This file (originally in Franz Lisp, by Richard Fateman, converted by me)
;;; replaces wholesale the previous version of POIS2. -- Krausz 6/3/88

;;; Modifications by Bruce Miller @ NIST 7/94:
;;;   A] 1) minimize use of special variables, 2) replace recursions over the length of 
;;; poisson series by iterations, 3) inverted the sense of a negative angular multiple
;;; to simplify the predicate (now v-u is normalized, not u-v), 4) for poismult, reuse
;;; the AVL tree structure to build the resulting series, 5) use byte functions (ldb & dpb)
;;; to encode/decode angles, rather than arithmetic, 6) other random compulsive recoding...
;;;  B] 1) More liberal approach upon encoding poisson series: If $pois_encode_liberalize 
;;; is non-null, then trig. arguments that are NOT linear in poisvars (w/ integer multiples)
;;; are accepted, but taken into the coefficients.
;;; 2) Added an extra bit between fields in packed trig arguments to allow detection of
;;; over/underflow in multiplication & substitution. Added checks against *pois-guard* to
;;; detect these errors.
;;;  C] Variable width fields for encoded arguments: the number of bits allocated for each
;;; angle's multiple can be specified separately.
;;;  D] 1) Option for general encoded arguments: If $poislim is NIL (ie. "no limit"), then 
;;; the encoded arguments are simply lists of the multipliers.  This removes the restriction
;;; on integer multipliers:  any symbolic quantity is allowed.
;;;   2) Replaced fancypoissubst (the 5 arg version of poissubst) to use taylor's
;;; method (JPFitch suggestion); This gives DIFFERENT result than original code. I cant figure
;;; out what the original was meant to calculate.


;; General Poisson 
(proclaim '(special poisco1 poiscom1		;ss cc h* rwg 5/90
						; *argc *coef b* a* *a brm 7/94
		    poishift $wtlvl *pois1* $poisvars))

;;; The following were removed because they are no longer needed. Further, with the
;;; altered encoding scheme, anyone using them will likely get the `wrong thing'
;;    poisvals poistsm poissiz poists

(defvar *pois-initialized* nil) ;; tested by an EVAL-WHEN at the end of the file.
;; default set of poisson vars.
(defvar *poisz* '((mpois simp) nil nil))
(defvar trim nil) ;; don't trim terms.

;;; Alist of (var byteposition zero)
(defvar *pois-encoding* nil)
;;; Guard bit mask for encoded trig arguments.
(defvar *pois-guard* nil)
;;; Controls whether poisencode accepts trig arguments which are not strictly linear 
;;; in the poisvars.  NIL (default) signals errors. T allows the non `poisvar' part to be
;;; taken into the coefficients.
(defvar $pois_encode_liberalize nil)

(defconstant *min-$poislim* 30 			; 32 is ok too (only *pois-guard* is big)
  "Largest economical value for $poislim")

;;; this tells the evaluator to keep out of poisson series.
(defprop mpois (lambda (x) x) mfexpr*)

;;;********************************************************************************
;;; Iteration across (& construction of) a single sine or cosine list of a Poisson.
(defun rtcons3 (mlt cf pois)
  (if (poispzero cf) pois
      (cons cf (cons mlt pois))))

;;; executes BODY for each term in the sine/cosine list POIS.  MULT & COEF are bound to the
;;; packed argument & coefficient for each term.
;;; The internal macro (COLLECT MLT COEF) can (but needn't) be used to construct 
;;; a new sine/cosine list during execution, which is returned as the result.
(defmacro mapping-pois ((pois mult coef) &body body)
  (let ((result (make-symbol "RESULT"))
	(ptr (make-symbol "PTR")))
    `(let ((,result nil))
       (macrolet ((collect (m c) (list 'setq ',result (list 'rtcons3 m c ',result))))
	 (do ((,ptr ,pois (cddr ,ptr)))
	     ((null ,ptr) (nreverse ,result))
	   (let ((,mult (car ,ptr))(,coef (cadr ,ptr)))
	     ,@body))))))

;;;********************************************************************************
;;; Setup parameters for Poisson series argument encoding.
;;;  POISVARS is a list of variables
;;;  POISLIM is either
;;;    NIL for general poisson argument lists
;;;    a number: for the total number of bits to allocate to the set of variables
;;;    an mlist of numbers: for the number of bits for each variable.
;;; The number of bits for each variable includes the sign & guard bits:
;;;  thus 4 bits allows for a range of -3...3.

(defun pois-setup (poisvars poislim &optional verbose)
  (let (total limits (n (length (cdr poisvars)))
	      (pos 0)(guard 0)(encoding nil)(shift 0))
    (cond ((null poislim) (setq shift (make-list n :initial-element 0))
	   (when verbose (mtell "General Poisson multipliers are being used for ~M" poisvars)))
	  ((and (fixnump poislim)(plusp poislim)) 
	   (setq total (max poislim *min-$poislim*) 
		 limits (make-list n :initial-element (quotient total n))))
	  ((and ($listp poislim)
		(every #'(lambda (b) (and (fixnump b)(plusp b))) (cdr poislim))
		(= n (length (cdr poislim))))
	   (setq limits (cdr poislim)
		 total (addn limits t))
	   (let ((d (quotient (- *min-$poislim* total) n)))
	     (when (plusp d)
	       (setq limits (mapcar #'(lambda (x)(+ x d)) limits)
		     total (+ total (* n d))))))
	  (t (merror "POISSON: Improper encoding; POISLIM was ~M,
 must be False, a positive fixnum or a list of ~M positive fixnums" poislim n)))
    (when limits
      (when verbose
	(mtell "Packing Poisson arguments into ~:M bits;~%Maximum multipliers are" total))
      (do ((vars (cdr poisvars) (cdr vars))
	   (lims limits   (cdr lims)))
	  ((null vars))
	(let* ((v (car vars)) (m (car lims))
	       (b (byte (1- m) pos))
	       (z (expt 2 (- m 2))))
	  (push (list v b z) encoding)
	  (incf pos m)
	  (setq shift (dpb z b shift)
		guard (dpb 1 (byte 1 (1- pos)) guard))
	  (when verbose
	    (mtell " +/-~:M for ~:M~A" (1- z) v
		   (cond ((null (cdr vars)) ".")((null (cddr vars)) " and")(t ",")))))))
    (setq $poisvars poisvars
	  $poislim  poislim
	  *pois-encoding* (nreverse encoding)
	  *pois-guard* guard
	  poishift shift
	  *poisz*  '((mpois simp) nil nil)	; Poisson 0
	  *pois1*   (list '(mpois simp) nil (list poishift poisco1)))	; Poisson 1
    '$done))

;;; This is called when $poislim is set.
(defun poislim1 (u poislim)
  (declare (ignore u))
  (pois-setup $poisvars poislim t))

;;; This is (or rather should be) called when $poisvars is set
(defun poisvars1 (u vars)
  (declare (ignore u))
  (when (and ($listp $poislim) (not (= (length $poislim)(length vars))))
    (mtell "Poisvars: Setting poislim back to ~M"	; Get consistent $poislim.
	   (setq $poislim (max *min-$poislim* (addn (cdr $poislim) t)))))
  (pois-setup vars $poislim t))

;;; I assume this should be in some other file? (like for $poislim)
(setf (get '$poisvars 'assign) 'poisvars1)

;;; An alternative Macsyma-callable interface to variable width pois setup:
;;; Give an alternating list of variables & maximum multiplier (not #bits!)
(defun $pois_setup (&rest vars-and-multipliers)
  (do ((x vars-and-multipliers (cddr x)) (v)(m))
      ((null x)(pois-setup (cons '(mlist simp)(nreverse v))(cons '(mlist simp)(nreverse m)) t))
    (push (car x) v)
    (push (+ 2 (ceiling (log (cadr x) 2))) m)))

;;;********************************************************************************
;;; Poisson argument (ie. angle) operations.
;;; these routines deal with the internals of the encoding of arguments.

;;; Extract a given multiple from an encoded argument.
(defun poisarg-varpos (anglevar)
  (if $poislim
      (cdr (assoc anglevar *pois-encoding* :test #'alike1))
      (position anglevar (cdr $poisvars) :test #'alike1)))

(defun poisarg-extract (mlt encoding)
  (if $poislim
      (- (ldb (car encoding) mlt) (cadr encoding))
      (nth encoding mlt)))
  
;;; Return a list of the angle's multipliers from CODE
(defun poisarg-unpack (mlt)
  (if $poislim
      (mapcar #'(lambda (entry) (- (ldb (cadr entry) mlt) (caddr entry))) *pois-encoding*)
      mlt))

;;; Decode an argument
(defun poisarg-decode (mlt) 
  (addn (if $poislim
	    (mapcar #'(lambda (entry)(mul (car entry)(- (ldb (cadr entry) mlt)(caddr entry))))
		    *pois-encoding*)
	    (mapcar #'mul (cdr $poisvars) mlt)) t))

;;; Check for overflows in argument encoding.
(defun poisarg-overflow (mlt)
  (let ((entry (do ((es *pois-encoding* (cdr es)))
		   ((not (zerop (ldb (byte 1 (+ (byte-position (second (car es)))
						(byte-size   (second (car es))))) mlt)))
		    (car es)))))
    (merror "POISSON: Overflow in encoding: coefficient of ~M outside of range +/- ~M"
	    (first entry)(third entry))))

(defmacro poisarg-check (mlt)
  (let ((arg (make-symbol "ARG")))
    `(let ((,arg ,mlt))
       (when (logtest ,arg *pois-guard*)
	 (poisarg-overflow ,arg))
       ,arg)))

;;; Normalizing an encoded argument
;;; And return 2nd value: whether it was negated
;;; [Note this is the reverse ordering of the original]

(defun poisarg-normalize (mlt)
  (if $poislim
      (if (< mlt poishift)
	  (values (poisarg-check (+ poishift (- poishift mlt))) t)
	  (values mlt nil))
      (do ((m mlt (cdr m)))
	  ((or (null m) (not (zerop1 (car m))))
	   (if (and m (eq -1 (signum1 (car m))))
	       (values (mapcar #'neg mlt) t)
	       (values mlt nil))))))

(defun poisarg-zerop  (mlt)
  (if $poislim
      (= mlt poishift)
      (every #'zerop1 mlt)))

;;; Compare two encoded arguments. Return :equal, :less or :greater
(defun poisarg-cmp (mlt1 mlt2)
  (if $poislim
      (cond ((= mlt1 mlt2) :equal)
	    ((< mlt1 mlt2) :less)
	    (t :greater))
      (do ((m1 mlt1 (cdr m1))
	   (m2 mlt2 (cdr m2)))
	  ((or (null m1)(not (alike1 (car m1)(car m2))))
	   (cond ((null m1) :equal)
		 ((lessthan (car m1)(car m2)) :less)
		 (t :greater))))))

;;; Compute the sum & difference of two encoded arguments.
(defun poisarg+/- (mlt1 mlt2)
  (if $poislim
      (let ((d (- mlt2 poishift)))
	(values (poisarg-check (+ mlt1 d))(poisarg-check (- mlt1 d))))
      (values (mapcar #'add mlt1 mlt2) (mapcar #'sub mlt1 mlt2))))

;;; For substituting VV for UU.
;;; This precomputes something to pass as the subst argument to poisarg-subst.
(defun poisarg-substprep (uu vv)
  (let ((diff (poisencode `((mplus) ,uu ((mtimes) -1 ,vv)))))
    (if $poislim
	(- diff poishift)
	diff)))

(defun poisarg-subst (mlt subst loc)
  (let ((m (poisarg-extract mlt loc)))
    (if $poislim
	(poisarg-check (+ mlt (* subst m)))
	(mapcar #'(lambda (mm dd)(add mm (* m dd))) mlt subst))))

;;; Check the legality of the trig argument; return encoding or (err) if an error
;;; Liberal version.
(defun poisencode (argmnt)
  (let* ((rr ($expand argmnt))
	 (mlt (if $poislim
		  (do ((entries *pois-encoding* (cdr entries))
		       (result 0))
		      ((null entries) result)
		    (let* ((v (caar entries))(b (cadar entries))(z (caddar entries))
			   (cf ($coeff rr v)))
		      (cond ((integerp cf)
			     (unless (< (abs cf) z)
			       (merror "POISSON: multiple of ~M exceeds limit +/-~M in ~M"
				       v (1- z) rr))
			     (setq result (dpb (- cf z) b result)
				   rr (sub rr (mul cf v))))
			    ($pois_encode_liberalize)
			    (t (merror "POISSON: multiple of ~M is non-integral in ~M"
				       v argmnt)))))
		  (do ((vs (cdr $poisvars) (cdr vs))
		       (result nil))
		      ((null vs) (nreverse result))
		    (let* ((v (car vs))(cf ($coeff rr v)))
		      (push cf result)
		      (setq rr (sub rr (mul cf v))))))))
    (unless (or (zerop1 rr) $pois_encode_liberalize)
      (if (some #'(lambda (v)(not (freeof v rr))) (cdr $poisvars))
	  (merror "POISSON: argument is non-linear: ~M" argmnt)
	  (merror "POISSON: argument contains non-Poisvars: ~M" argmnt)))
    (values mlt rr)))

;;;********************************************************************************	
;;; Encoding & Decoding Poisson series

(defun $poissimp (x)
  (if (mbagp x) (cons (car x) (mapcar #'$poissimp (cdr x))) ($outofpois x)))

(defun $intopois (x) (intopois x)) 

(defun intopois (x) 
  (if (atom x)
      (if (equal x 0)
	  *poisz*
	  (list '(mpois simp) nil (list poishift (intopoisco x))))
      (case (caar x)
	(mpois x)
	(%sin (poissine (cadr x)))
	(%cos (poiscosine (cadr x)))
	(mexpt (if (and (numberp (caddr x)) (greaterp (caddr x) 0))
		   ($poisexpt (intopois (cadr x)) (caddr x))
		   (list '(mpois simp) nil (list poishift (intopoisco x)))))
	(mplus  (lisp:reduce #'$poisplus (cdr x) :from-end t))
	(mtimes (lisp:reduce #'$poistimes (cdr x)))
	(mrat   (intopois (ratdisrep x)))
	; shouldn't mbags be handled here?
	(t (list '(mpois simp) nil (list poishift (intopoisco x)))))))

;;; Liberal versions.
(defun poiscosine (angle) 
  (multiple-value-bind (mlt rest)(poisencode angle)
    (if (zerop1 rest)
	(list '(mpois simp) nil (list (poisarg-normalize mlt) poisco1))
	(multiple-value-bind (mlt negp)(poisarg-normalize mlt)
	  (list '(mpois simp)
		(list mlt (intopoisco (if negp ($sin rest) (neg ($sin rest)))))
		(list mlt (intopoisco ($cos rest))))))))

(defun poissine (angle) 
  (multiple-value-bind (mlt rest)(poisencode angle)
    (multiple-value-bind (mlt negp)(poisarg-normalize mlt)
      (if (zerop1 rest)
	  (list '(mpois simp) (list mlt (if negp poiscom1 poisco1)) nil)
	  (list '(mpois simp)
		(list mlt (intopoisco (if negp (neg ($cos rest))($cos rest))))
		(list mlt (intopoisco ($sin rest))))))))

;;;********************************************************************************	
;;; Mapping over Poisson series

(defun $poismap (pois sinfn cosfn) 
  (setq pois (intopois pois))
  (flet ((pm1 (list fcn)
	   (mapping-pois (list mmm ccc)
	     (collect mmm (mfuncall fcn (poiscdecode ccc) (poisarg-decode mmm))))))
    (list (car pois) (pm1 (cadr pois) sinfn) (pm1 (caddr pois) cosfn))))

;;;********************************************************************************	
;;; Decoding & Printing Poisson series

(defun $outofpois (pois) 
  (when (or (atom pois) (not (eq (caar pois) 'mpois)))
    (setq pois (intopois pois)))
  (let ((result nil))
    (mapping-pois ((cadr pois) mm cc)    ; do sines
      (push (mul (poiscdecode cc)(simplify (list '(%sin) (poisarg-decode mm)))) result))
    (mapping-pois ((caddr pois) mm cc)    ; do cosines
      (push (if (poisarg-zerop mm) (poiscdecode cc)
		(mul (poiscdecode cc)(simplify (list '(%cos) (poisarg-decode mm))))) result))
    (addn result t)))

(defun $printpois (pois)
  (let (($pfeformat t))				; typically there are lots of fractions;
    (declare (special $pfeformat))
    (setq pois (intopois pois))
    (mapping-pois ((cadr pois) mm cc)		; do sines
      (displa (mul (poiscdecode cc) (list '(%sin) (poisarg-decode mm)))))
    (mapping-pois ((caddr pois) mm cc)		; do cosines
      (displa (if (poisarg-zerop mm) (poiscdecode cc)
		  (mul (poiscdecode cc)(list '(%cos) (poisarg-decode mm))))))
    '$done))

;;;********************************************************************************	
;;; Addition of Poisson series

(defun $poisplus (x &rest ys) 
  (flet ((merge22 (x z)
	   (do ((result nil))
	       ((or (null x)(null z)) (nconc (nreverse result) (or z x)))
	     (case (poisarg-cmp (car x)(car z))
	       (:equal (let ((tt (poisco+ (cadr x) (cadr z))))
			 (unless (poispzero tt)
			   (setq result (cons tt (cons (car x) result))))
			 (setq x (cddr x) z (cddr z))))
	       (:less (setq result (cons (cadr x)(cons (car x) result)) x (cddr x)))
	       (t     (setq result (cons (cadr z)(cons (car z) result)) z (cddr z)))))))
    (setq x (intopois x))
    (let ((sins (cadr x))(coss (caddr x)))
      (dolist (yy ys)
	(setq yy (intopois yy)
	      sins (merge22 sins (cadr yy))
	      coss (merge22 coss (caddr yy))))
      (list '(mpois simp) sins coss))))

(defun $poisminus (x &rest ys)
  ($poisplus x ($poisctimes -1 (apply #'$poisplus ys))))

;;;********************************************************************************	
;;; Scaling:
;;; this program multiplies a poisson series p by a non-series, c,
;;; which is free of sines and cosines .


(defun poisctimes1 (x pois) 
  (flet ((pt1 (q)
	   (mapping-pois (q mlt cf)
	     (unless (and trim (trimf mlt))	; This gets set when called via $poismult
	       (collect mlt (poisco* x cf))))))
    (list '(mpois simp) (pt1 (cadr pois))(pt1 (caddr pois)))))

(defun $poisctimes (x pois) 
  (poisctimes1 (intopoisco x) (intopois pois)))

;;;********************************************************************************	
;;; $poisdiff differentiates a poisson series wrt 
;;; an element on  $poisvars (default x, y, z, u, v, w), or a coeff var.

(defun $poisdiff (pois var)
  (setq pois (intopois pois))
  (if (mexp-member var $poisvars)
      (let ((loc (poisarg-varpos var)))
	(flet ((d1 (q sgn)
		 (mapping-pois (q mm cc)
		   (collect mm (poisco* (intopoisco (* sgn (poisarg-extract mm loc))) cc)))))
	  (list (car pois) (d1 (caddr pois) -1) (d1 (cadr pois) +1))))
      (flet ((d1 (q)
	       (mapping-pois (q mm cc)
		 (collect mm (poiscodif cc var)))))
	(list (car pois) (d1 (cadr pois)) (d1 (caddr pois))))))

;;;********************************************************************************	
;;; $poisint integrates a poisson series wrt x,y, z, u, v, w.  the variable of
;;; integration must occur only in the arguments of sin or cos,
;;; or only in the coefficients.  poiscointeg is called to integrate coeffs.

;;; non-periodic terms are removed.

(defun $poisint (pois var) 
  (setq pois (intopois pois))
  (if (mexp-member var (cdr $poisvars))
      (let ((loc (poisarg-varpos var)))
	(flet ((i1 (q sgn)
		 (mapping-pois (q mm cc)
		   (let ((j (poisarg-extract mm loc)))
		     (unless (zerop1 j)
		       (collect mm (poisco* (intopoisco (list '(mexpt)(* j sgn) -1)) cc)))))))
	  (list (car pois) (i1 (caddr pois) +1) (i1 (cadr pois) -1))))
      (flet ((i1 (q)
	       (mapping-pois (q mm cc)
		 (collect mm (poiscointeg cc var)))))
	(list (car pois)(i1 (cadr pois))(i1 (caddr pois))))))
 
;;;********************************************************************************	
;;; Multiplication 

;;; avl balanced tree search and insertion.
;;; node looks like (key (llink .  rlink) balancefactor .  record)
;;; program follows algorithm given in knuth vol. 3 455-57

;; Try this: ( (key record . balancefactor) llink . rlink)
;; then, the poisson list can reuse the avl structure
;; (balance factor replaced by link to next term).
;; Note: balancefactor is always 1,0 or -1

(defmacro make-avlnode (key record &optional llink rlink)
  `(list* (list* ,key ,record 0) ,llink ,rlink))

;; macros to extract fields from node

(DEFMACRO KEY   (X) `(caar   ,X))
(DEFMACRO LLINK (X) `(cadr ,X)) 
(DEFMACRO RLINK (X) `(cddr ,X))
(DEFMACRO BP    (X) `(cddar ,X))
(DEFMACRO REC   (X) `(cadar ,X))

(defun avlinsert (key newrec head)
  (let* (qq (tt head) (pp (rlink head))(zz pp) kk bp)
						; This should be a loop with a case statement on poisarg-cmp.
						; But this is a tight, frequently called loop in multiplication, so at cost of modularity
						; I've open coded the branch here.
    (if $poislim
	(prog () loop
	      (setq kk (key pp))
	      (cond ((= key kk) (setf (rec pp) (poisco+ (rec pp) newrec))
		     (return-from avlinsert head))
		    ((< key kk) (unless (setq qq (llink pp))
				  (setf (llink pp) (make-avlnode key newrec))(return)))
		    (t (unless (setq qq (rlink pp))
			 (setf (rlink pp) (make-avlnode key newrec))(return))))
	      (unless (zerop (bp qq))
		(setq tt pp zz qq))
	      (setq pp qq)
	      (go loop))
	(prog () loop
	      (do ((m1 key (cdr m1))
		   (m2 (key pp) (cdr m2)))
		  ((or (null m1)(not (alike1 (car m1)(car m2))))
		   (cond ((null m1) (setf (rec pp) (poisco+ (rec pp) newrec))
			  (return-from avlinsert head))
			 ((lessthan (car m1)(car m2)) (unless (setq qq (llink pp))
							(setf (llink pp) (make-avlnode key newrec))))
			 (t (unless (setq qq (rlink pp))
			      (setf (rlink pp) (make-avlnode key newrec)))))))
	      (unless qq (return))
	      (unless (zerop (bp qq))
		(setq tt pp zz qq))
	      (setq pp qq)
	      (go loop)))
    (setq qq (setq pp (if (eq :less (poisarg-cmp key (key zz))) (llink zz) (rlink zz))))
    (prog () loop
	  (case (poisarg-cmp key (key pp))
	    (:equal (return))
	    (:less (setf (bp pp) -1 pp (llink pp)))
	    (t     (setf (bp pp)  1 pp (rlink pp))))
	  (go loop))
    (if (eql :less (poisarg-cmp key (key zz)))
	(case (bp zz)
	  (0 (setf (bp zz) -1) (incf (llink head)))
	  (1 (setf (bp zz) 0))
	  (t (if (= (bp qq) -1)
		 (setf pp         qq
		       (llink zz) (rlink qq)
		       (rlink qq) zz
		       (bp zz)    0
		       (bp qq)    0)
		 (setf pp         (rlink qq)
		       (rlink qq) (llink pp)
		       (llink pp) qq
		       (llink zz) (rlink pp)
		       (rlink pp) zz
		       bp         (bp pp)
		       (bp zz)    (if (eql bp -1) 1 0)
		       (bp qq)    (if (eql bp 1) -1 0)
		       (bp pp)    0))
	     (if (eq zz (rlink tt)) (setf (rlink tt) pp) (setf (llink tt) pp))))
	(case (bp zz)
	  (0  (setf (bp zz) 1) (incf (llink head)))
	  (-1 (setf (bp zz) 0))
	  (t (if (= (bp qq) 1)
		 (setf pp         qq
		       (rlink zz) (llink qq)
		       (llink qq) zz
		       (bp zz)    0
		       (bp qq)    0)
		 (setf pp         (llink qq)
		       (llink qq) (rlink pp)
		       (rlink pp) qq
		       (rlink zz) (llink pp)
		       (llink pp) zz
		       bp         (bp pp)
		       (bp zz)    (if (eql bp 1) -1 0)
		       (bp qq)    (if (eql bp -1) 1 0)
		       (bp pp)    0))
	     (if (eq zz (rlink tt)) (setf (rlink tt) pp) (setf (llink tt) pp)))))
    head))

(defun avlinit (key rec) 
  (make-avlnode 'top nil 0 (make-avlnode key rec))) 

;; untree converts the tree to a list which looks like
;;  ( smallest-key record next-smallest-key record ....  largest-key record)

(defun untree (tree)
  (let ((result nil))
    (labels ((un1 (tree)
	       (when tree
		 (un1 (rlink tree))
		 (unless (poispzero (rec tree))
		   (setf (cddar tree) result)
		   (setq result (car tree)))
		 (un1 (llink tree)))))
      (un1 (rlink tree))
      result)))

;;;********************************************************************************	
(defun nonperiod (p)
  (and (null (cadr p)) (poisarg-zerop (caaddr p)) (null (cddr (caddr p))))) 

(defun trimf (mlt)
  (mapply '$poistrim (poisarg-unpack mlt) '$poistrim))

(defun $poistimes (a &rest more) 
  (let ((trim (or (getl '$poistrim '(expr subr)) (mget '$poistrim 'mexpr))))
    (setq a (intopois a))
    (dolist (b more)
      (setq a (poistimes2 a (intopois b))))
    a))

(defun poistimes2 (poisa poisb) 
  (cond ((and (null (cadr poisa)) (null (caddr poisa))) (return-from poistimes2 poisa))
	((and (null (cadr poisb)) (null (caddr poisb))) (return-from poistimes2 poisb))
	((nonperiod poisa) (return-from poistimes2 (poisctimes1 (cadr (caddr poisa)) poisb)))
	((nonperiod poisb) (return-from poistimes2 (poisctimes1 (cadr (caddr poisb)) poisa))))
  (let* ((zero (intopoisco 0.))
	 (slc (avlinit poishift zero))
	 (clc (avlinit poishift zero)))
    (do ((las (cdr poisa)(cdr las))		; first sines of a, then cosines
	 (asinp t nil))				; Are elements of lista sines?
	((null las))
      (mapping-pois ((car las) mlta cfaa)	; through each sine/cosine in la
	(do ((lbs (cdr poisb)(cdr lbs))		; first sines of b, then cosines
	     (cfa (halve cfaa))
	     (bsinp t nil))			; is B list of sines?
	    ((null lbs))
	  (let* ((sinp (xor asinp bsinp))	;result is sine?
		 (lc   (if sinp slc clc))
		 (neg- (and (not asinp) bsinp))	; does u-v get - sign?
		 (neg+ (and asinp bsinp)))	; does u+v get - sign?
	    (mapping-pois ((car lbs) mltb cfb)	; each sine/cosine of b
	      (let ((cc (poisco* cfa cfb)))
		(unless (poispzero cc)
		  (multiple-value-bind (mlt+ mlt-)(poisarg+/- mlta mltb)
		    (unless (or (and trim (trimf mlt-))	; trim this term?
				(and sinp (poisarg-zerop mlt-)))	; sin(0)
		      (multiple-value-bind (mm negp)(poisarg-normalize mlt-)
			(avlinsert mm (if (xor (and negp sinp) neg-)
					  (poisco* poiscom1 cc) cc)  lc)))
		    (unless (and trim (trimf mlt+))	; trim this term?
		      (avlinsert mlt+ (if neg+ (poisco* poiscom1 cc) cc) lc))))))))))
    (list '(mpois simp) (untree slc) (untree clc))))

(defun $poisexpt (pois pow) 
  (prog (result one) 
	(if (or (not (integerp pow)) (minusp pow))
	    (improper-arg-error '$poisexpt pow))
	(setq result (if (oddp pow) pois (setq one (intopois 1))))
     a    (setq pow (lsh pow -1))
	(if (zerop pow) (return result))
	(setq pois ($poistimes pois pois))
	(if (oddp pow)
	    (setq result (if (equal result one) pois ($poistimes result pois))))
	(go a)))

(defun $poissquare (x) ($poistimes x x))

;;;********************************************************************************	
;;; $poissubst substitutes an expression for a variable in argument of trig functions or
;;; coefficients.

;;; From KMP's version...
(defun $poissubst (uu vv ww &rest xx-and-oo)
  (case (length xx-and-oo)
    (0 (let ((ww (intopois ww)))
	 (if (mexp-member vv (cdr $poisvars))
	     (poissubsta uu vv ww)
	     (flet ((s1 (q)
		      (mapping-pois (q mm cc)
			(collect mm (poissubstco uu vv cc)))))
	       (list (car ww) (s1 (cadr ww)) (s1 (caddr ww)))))))
    (2 (fancypoissubst uu vv (intopois ww) (intopois (car xx-and-oo)) (cadr xx-and-oo)))
    (otherwise (wrong-number-of-args-error '$poissubst))))

(defun poissubsta (uu vv ww)
  (let* ((loc (poisarg-varpos vv))
	 (subst (poisarg-substprep uu vv))
	 (zero (intopoisco 0.))
	 (sins (avlinit poishift zero))
	 (coss (avlinit poishift zero)))
    (mapping-pois ((cadr ww) mm cf)
      (let ((mm (poisarg-subst mm subst loc)))
	(unless (poisarg-zerop mm)		; sine(0) is 0
	  (multiple-value-bind (mm negp)(poisarg-normalize mm)
	    (avlinsert mm (if negp (poisco* poiscom1 cf) cf) sins)))))
    (mapping-pois ((caddr ww) mm cf)
      (avlinsert (poisarg-normalize (poisarg-subst mm subst loc)) cf coss))
    (list (car ww)(untree sins)(untree coss))))

; I think this is what poissubst was intended to do.
; I've followed J.P.Fitch's suggestion: Given that poisson series are easy to differentiate,
; computing the expansion by Taylor's method is quite reasonable.
(defun fancypoissubst (uu vv ww xx oo)
  (let ((xx^n 1)
	(dwwdvv ww)
	(result (poissubsta uu vv ww)))		; 0-th order
    (dotimes (i oo)
      (setq xx^n ($poisctimes (list '(rat) 1 (1+ i))($poistimes xx^n xx))
	    dwwdvv ($poisdiff dwwdvv vv)
	    result ($poisplus result ($poistimes xx^n (poissubsta uu vv dwwdvv)))))
    result))

;;;********************************************************************************	

; the following comment fragment predates my marauding --rwg 5/90
;argument  do not exceed some predefined bound in absolute value

;;;********************************************************************************	
;;; these are the only coefficient dependent routines. 

;;; poiscdecode decodes a coefficient

(defun poiscdecode (x) x) 

;;; intopoisco puts an expression into poisson coefficient form

(defun intopoisco (x) (simplifya x nil)) 

;;; poisco+ adds 2 coefficients

(defun poisco+ (x z) (add x z)) 

;;; poisco* multiplies 2 coefficients

(defun poisco* (x z) (mul x z)) 

;;; halve divides a coefficient by 2

(defun halve (x) (mul  '((rat) 1. 2.) x))

;;; poissubstco substitutes an expression for a variable in a coefficient.

(defun poissubstco (uu vv ww) (substitute uu vv ww)) 

;;; this differentiates a coefficient

(defun poiscodif (cf vv) (sdiff cf vv))

;;; this integrates a coefficient
(defun poiscointeg (cf vv) (intopoisco ($integrate (poiscdecode cf) vv)))

;;; test for zero

(defun poispzero (x) (zerop1 x)) 

;;; the number 1 in coefficient arithmetic, the number -1

(defvar poisco1 1)
(defvar poiscom1 -1)
(setq poisco1 1 poiscom1 -1)

;;; If the user has not set poislim, do it.
(eval-when (load eval)
  (when (null *pois-initialized*)
    (poislim1 nil 30.)
    (setq *pois-initialized* t)))

;;; -------------------------------------------------------------------------