src/pois3.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancments.                    ;;;;;
   4 ;;;                                                                    ;;;;;
   5 ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
   6 ;;;     All rights reserved                                            ;;;;;
   7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   8 ;;;     (c) Copyright 1981 Massachusetts Institute of Technology         ;;;
   9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10
  11 (in-package :maxima)
  12
  13 (macsyma-module pois3)
  14
  15 ;; GENERAL POISSON SERIES
  16
  17 (declare-top (special *argc *coef poisvals b* a* *a ss cc h* poishift
  18                       poistsm poists $poisz $pois1))
  19
  20 (defvar trim nil)
  21
  22 ;;; THESE ARE THE ONLY COEFFICIENT DEPENDENT ROUTINES.
  23
  24 ;;; POISCDECODE DECODES A COEFFICIENT
  25 (defun poiscdecode (x) x)
  26
  27 ;;; INTOPOISCO PUTS AN EXPRESSION INTO POISSON COEFFICIENT FORM
  28 (defun intopoisco (x) (simplifya x nil))
  29
  30 ;;; POISCO+ ADDS 2 COEFFICIENTS
  31 (defun poisco+ (r s) (simplifya (list '(mplus) r s) nil))
  32
  33 ;;; POISCO* MULTIPLIES 2 COEFFICIENTS
  34 (defun poisco* (r s) (simplifya (list '(mtimes) r s) nil))
  35
  36 ;;; HALVE DIVIDES A COEFFICIENT BY 2
  37 (defun halve (r)
  38   (simplifya (list '(mtimes) '((rat) 1 2) r) nil))
  39
  40 ;;; POISSUBSTCO SUBSTITUTES AN EXPRESSION FOR A VARIABLE IN A COEFFICIENT.
  41 (defun poissubstco (a b c)
  42   (maxima-substitute a b c))
  43
  44 ;;; THIS DIFFERENTIATES A COEFFICIENT
  45 (defun poiscodif (h var)
  46   ($diff h var))
  47
  48 ;;; THIS INTEGRATES A COEFFICIENT
  49 (defun poiscointeg (h var)
  50   (intopoisco($integrate (poiscdecode h) var)))
  51
  52 ;;; TEST FOR ZERO
  53 (defun poispzero (x) (zerop1 x))
  54
  55 (defun fumcheck (x)
  56   (not (and (atom x) (integerp x) (< (abs x) poistsm))))
  57
  58 (defun checkencode(r)
  59   (prog(q)
  60      (setq q ($coeff r '$u))
  61      (cond ((fumcheck q) (return nil))
  62            (t (setq r (simplifya (list '(mplus) r (list '(mtimes) -1  '$u q)) nil))))
  63      (setq q ($coeff r '$v))
  64      (cond ((fumcheck q)(return nil))
  65            (t (setq r (simplifya (list '(mplus) r (list '(mtimes) -1  '$v q)) nil))))
  66      (setq q ($coeff r '$w))
  67      (cond ((fumcheck q)(return nil))
  68            (t (setq r (simplifya (list '(mplus) r (list '(mtimes) -1  '$w q)) nil))))
  69      (setq q ($coeff r '$x))
  70      (cond ((fumcheck q)(return nil))
  71            (t (setq r (simplifya (list '(mplus) r (list '(mtimes) -1  '$x q)) nil))))
  72      (setq q ($coeff r '$y))
  73      (cond ((fumcheck q)(return nil))
  74            (t (setq r (simplifya (list '(mplus) r (list '(mtimes) -1  '$y q)) nil))))
  75      (setq q ($coeff r '$z))
  76      (cond ((fumcheck q)(return nil))
  77            (t (setq r (simplifya (list '(mplus) r (list '(mtimes) -1  '$z q)) nil))))
  78      (cond ((equal r 0)(return t))
  79            (t (return nil)))))
  80
  81 (defmfun $poissimp (x)
  82   (if (mbagp x)
  83       (cons (car x) (mapcar #'$poissimp (cdr x)))
  84       ($outofpois x)))
  85
  86 ;;;********
  87
  88 ;; ABOVE ASSUMES POISLIM(5) OR LESS ALSO REDEFINE ORDER< AND ORDER= TO BE < AND =
  89
  90 ;;; THIS TELLS THE EVALUATOR TO KEEP OUT OF POISSON $SERIES.
  91
  92 (defprop mpois (lambda (x) x) mfexpr*)
  93
  94 (defmfun $poisplus (a b)
  95   (setq a (intopois a) b (intopois b))
  96   (list '(mpois simp) (poismerge22 (cadr a) (cadr b)) (poismerge22 (caddr a) (caddr b))))
  97
  98 (declare-top (special *b *fn))
  99
 100 (defmfun $poismap (p sinfn cosfn)
 101   (prog (*b *fn)
 102      (setq p (intopois p))
 103      (setq *fn (list sinfn))
 104      (return (list (car p)
 105                    (poismap (cadr p))
 106                    (prog2 (setq *fn (list cosfn)) (poismap (caddr p)))))))
 107
 108 (defun poismap (y)
 109   (cond ((null y) nil)
 110         (t (setq *b (meval (list *fn (poiscdecode (cadr y)) (poisdecodec (car y)))))
 111            (tcons3(car y) (intopoisco  *b) (poismap (cddr y))))))
 112
 113 (defun poismerge22 (r s)
 114   (cond ((null r) s)
 115         ((null s) r)
 116         ((equal (car r) (car s))
 117          (prog (tt)
 118             (setq tt (poisco+ (cadr r) (cadr s)))
 119             (return (cond ((poispzero tt) (poismerge22 (cddr r) (cddr s)))
 120                           (t (cons (car s) (cons tt (poismerge22 (cddr r) (cddr s)))))))))
 121         ((< (car r) (car s)) (cons (car r) (cons (cadr r) (poismerge22 (cddr r) s))))
 122         (t (cons (car s) (cons (cadr s) (poismerge22 (cddr s) r))))))
 123
 124 (defun poiscosine (m)
 125   (setq m (poisencode m))
 126   (cond ((poisnegpred m) (setq m (poischangesign m))))
 127   (list '(mpois simp) nil (list m 1)))
 128
 129 (defun poissine (m)
 130   (setq m (poisencode m))
 131   (cond ((poisnegpred m) (list '(mpois simp) (list (poischangesign m) -1) nil))
 132         (t (list '(mpois simp) (list m 1) nil))))
 133
 134 (defmfun $intopois (x)
 135   (let (*a)
 136     (intopois x)))
 137
 138 (defun intopois (a)
 139   (cond ((atom a)
 140          (cond ((equal a 0) $poisz) (t (list '(mpois simp) nil (list poishift (intopoisco a))))))
 141         ((eq (caar a) 'mpois) a)
 142         ((eq (caar a) '%sin) (poissine (cadr a)))
 143         ((eq (caar a) '%cos) (poiscosine (cadr a)))
 144         ((and (eq (caar a) 'mexpt) (numberp (caddr a)) (> (caddr a) 0.))
 145          ($poisexpt (intopois (cadr a)) (caddr a)))
 146         ((eq (caar a) 'mplus)
 147          (setq *a (intopois (cadr a)))
 148          (mapc (function (lambda (z) (setq *a ($poisplus *a (intopois z))))) (cddr a))
 149          *a)
 150         ((eq (caar a) 'mtimes)
 151          (setq *a (intopois (cadr a)))
 152          (mapc (function (lambda (z) (setq *a ($poistimes *a (intopois z))))) (cddr a))
 153          *a)
 154         ((eq (caar a) 'mrat) (intopois (ratdisrep a)))
 155         (t (list '(mpois simp) nil (list poishift (intopoisco a))))))
 156
 157 (defun tcons (r s)
 158   (if (poispzero (car s))
 159       (cdr s)
 160       (cons r s)))
 161
 162 (defun poisnegpred ($n)
 163   (prog ($r)
 164    $loop (cond ((equal $n 0) (return nil))
 165                (t nil))
 166    (setq $r (- (rem $n poists) poistsm))
 167    (cond ((> $r 0) (return nil))
 168          ((> 0 $r) (return t))
 169          (t (setq $n (quotient $n poists))))
 170    (go $loop)))
 171
 172 (defun poischangesign ($n)
 173   (- (* poishift 2) $n))
 174
 175 (declare-top (special $u $v $w $x $y $z))
 176
 177 (defun poisencode (h*)
 178   (unless (checkencode h*)
 179     ;; NOT CLEAR WHAT IS ILLEGAL HERE
 180     (merror (intl:gettext "poissimp: illegal argument: ~M") h*))
 181   (apply #'(lambda ($z $y $x $w $v $u)
 182              (declare (special $u $v $w $x $y $z))
 183              (setq h* (meval h*))
 184              ;; NOT CLEAR WHAT IS ILLEGAL HERE EITHER
 185              (unless (integerp h*) (merror  (intl:gettext "poisson: illegal trigonometric argument.")))
 186              (+ poishift  h*))
 187          poisvals))
 188
 189 (let ((n 5))
 190    (setq poists (expt 2 n)
 191          poisvals (loop for i from 5 downto 0 collect (expt poists i))
 192          poistsm (expt 2 (1- n))
 193          poishift (loop for i from 0 to 5 sum (* poistsm (expt poists i)))
 194          $poisz '((mpois simp) nil nil)
 195          $pois1 (list '(mpois simp) nil (list poishift 1)))
 196    n)
 197
 198 (defun poisdecodec (m)
 199   (prog (arg h)
 200      (setq h m)
 201      (setq arg (list '(mtimes) (- (rem h poists) poistsm) '$u))
 202      (setq h (quotient h poists))
 203      (setq arg
 204            (list '(mplus)
 205                  arg
 206                  (list '(mtimes) (- (rem h poists) poistsm) '$v)))
 207      (setq h (quotient h poists))
 208      (setq arg
 209            (list '(mplus)
 210                  arg
 211                  (list '(mtimes) (- (rem h poists) poistsm) '$w)))
 212      (setq h (quotient h poists))
 213      (setq arg
 214            (list '(mplus)
 215                  arg
 216                  (list '(mtimes) (- (rem h poists) poistsm) '$x)))
 217      (setq h (quotient h poists))
 218      (setq arg
 219            (list '(mplus)
 220                  arg
 221                  (list '(mtimes) (- (rem h poists) poistsm) '$y)))
 222      (setq h (quotient h poists))
 223      (setq arg
 224            (list '(mplus)
 225                  arg
 226                  (list '(mtimes) (- (rem h poists) poistsm) '$z)))
 227      (return (simplifya arg nil))))
 228
 229
 230 ;;; THIS PROGRAM MULTIPLIES A POISSON SERIES P BY A NON-SERIES, C,
 231 ;;; WHICH IS FREE OF SINES AND COSINES .
 232
 233 (defmfun $poisctimes (c p)
 234   (list '(mpois simp) (poisctimes1 (setq c (intopoisco c)) (cadr p)) (poisctimes1 c (caddr p))))
 235
 236 (defmfun $outofpois (p)
 237   (prog (ans)
 238      (cond ((or (atom p) (not (eq (caar p) 'mpois))) (setq p (intopois p))))
 239
 240      ;; DO SINES
 241      (do ((m
 242            (cadr p)
 243            (cddr m)))(
 244                       (null m))
 245        (setq ans (cons (list '(mtimes)
 246                              (poiscdecode (cadr m))
 247                              (list '(%sin) (poisdecodec (car m))))
 248                        ans)))
 249
 250      ;; DO COSINES
 251      (do ((m
 252            (caddr p)
 253            (cddr m)))(
 254                       (null m))
 255        (setq ans (cons (list '(mtimes)
 256                              (poiscdecode (cadr m))
 257                              (cond ((equal (car m) poishift) 1)
 258                                    (t (list '(%cos) (poisdecodec (car m))))))
 259                        ans)))
 260      (return (cond ((null ans) 0.) (t (simplifya (cons '(mplus) ans) nil))))))
 261
 262 (defmfun $printpois (p)
 263   (prog nil
 264      (setq p (intopois p))
 265
 266      ;; DO SINES
 267      (do ((m
 268            (cadr p)
 269            (cddr m)))(
 270                       (null m))
 271        (displa (simplifya (list '(mtimes)
 272                                 (poiscdecode (cadr m))
 273                                 (list '(%sin) (poisdecodec (car m))))
 274                           t))
 275        (terpri)
 276        (finish-output))
 277
 278      ;; DO COSINES
 279      (do ((m
 280            (caddr p)
 281            (cddr m)))(
 282                       (null m))
 283        (displa (simplifya (list '(mtimes)
 284                                 (poiscdecode (cadr m))
 285                                 (cond ((equal (car m) poishift) 1.)
 286                                       (t (list '(%cos) (poisdecodec (car m))))))
 287                           t))
 288        (terpri)
 289        (finish-output))
 290      (return '$done)))
 291
 292
 293 ;;; $POISDIFF DIFFERENTIATES A POISSON SERIES WRT X, Y, Z, U, V, W, OR A COEFF VAR.
 294
 295
 296 (defmfun $poisdiff (p m)
 297   (declare (special m))
 298   (cond ((member m '($u $v $w $x $y $z) :test #'eq)
 299          (list (car p) (cosdif (caddr p) m) (sindif (cadr p) m)))
 300         (t (list (car p) (poisdif4(cadr p)) (poisdif4 (caddr p))))))
 301
 302
 303 (defun poisdif4 (y)
 304   (declare (special m))
 305   (cond ((null y) nil)
 306         (t (tcons3 (car y)(poiscodif (cadr y) m) (poisdif4 (cddr y))))))
 307
 308
 309 ;;; COSDIF DIFFERENTIATES COSINES TO GET SINES
 310
 311 (defun cosdif (h m)
 312   (cond ((null h) nil)
 313         (t (tcons (car h)
 314                   (cons (poisco* (intopoisco (- (poisxcoef (car h) m))) (cadr h))
 315                         (cosdif (cddr h) m))))))
 316
 317 (defun sindif (h m)
 318   (cond ((null h) nil)
 319         (t (tcons (car h)
 320                   (cons (poisco* (intopoisco (poisxcoef (car h) m)) (cadr h)) (sindif (cddr h) m))))))
 321
 322 (defun poisxcoef (h m)
 323   (- (rem (quotient h
 324                     (expt poists
 325                           (cadr (member m '($u 0 $v 1 $w 2 $x 3 $y 4 $z 5)))))
 326           poists)
 327      poistsm))
 328
 329
 330 ;;; AVL BALANCED TREE SEARCH AND INSERTION.
 331 ;;; NODE LOOKS LIKE (KEY (LLINK .  RLKINK) BALANCEFACTOR .  RECORD)
 332 ;;; PROGRAM FOLLOWS ALGORITHM GIVEN IN KNUTH VOL. 3 455-57
 333
 334 (declare-top (special ans))
 335
 336
 337 ;; MACROS TO EXTRACT FIELDS FROM NODE
 338
 339 (defmacro key  (&rest l) (cons 'car l))
 340
 341 (defmacro llink  (&rest l) (cons 'caadr l))
 342
 343 (defmacro rlink  (&rest l) (cons 'cdadr l))
 344
 345 (defmacro bp  (&rest l) (cons 'caddr l))
 346
 347 (defmacro rec  (&rest l) (cons 'cdddr l))
 348
 349
 350 ;;  FOR ORDERING KEYS
 351
 352 (defmacro order<  (&rest l) (cons '<  l))
 353 (defmacro order=  (&rest l) (cons '=  l))
 354
 355 ;; MACROS TO SET FIELDS IN NODE
 356
 357 (defmacro setrlink  (&rest l) (setq l (cons nil l))
 358           (list 'rplacd (list 'cadr (cadr l)) (caddr l)))
 359
 360 (defmacro setllink  (&rest l) (setq l (cons nil l))
 361           (list 'rplaca (list 'cadr (cadr l)) (caddr l)))
 362
 363 (defmacro setbp  (&rest l) (setq l (cons nil l))
 364           (list 'rplaca (list 'cddr (cadr l)) (caddr l)))
 365
 366 (defmacro setrec  (&rest l)(setq l (cons nil l))
 367           (list 'rplacd (list 'cddr (cadr l)) (caddr l)))
 368
 369
 370 (defun insert-it (pp newrec) (setrec pp (poisco+ (rec pp) newrec)))
 371
 372 (defun avlinsert (k newrec head)
 373   (prog (qq tt ss pp rr)
 374      (setq tt head)
 375      (setq ss (setq pp (rlink head)))
 376      a2   (cond ((order< k (key pp)) (go a3))
 377                 ((order< (key pp) k) (go a4))
 378                 (t (insert-it pp newrec) (return head)))
 379      a3   (setq qq (llink pp))
 380      (cond ((null qq) (setllink pp (cons k (cons (cons nil nil) (cons 0. newrec)))) (go a6))
 381            ((order= 0. (bp qq)) nil)
 382            (t (setq tt pp ss qq)))
 383      (setq pp qq)
 384      (go a2)
 385      a4   (setq qq (rlink pp))
 386      (cond ((null qq) (setrlink pp (cons k (cons (cons nil nil) (cons 0. newrec)))) (go a6))
 387            ((order= 0 (bp qq)) nil)
 388            (t (setq tt pp ss qq)))
 389      (setq pp qq)
 390      (go a2)
 391      a6   (cond ((order< k (key ss)) (setq rr (setq pp (llink ss)))) (t (setq rr (setq pp (rlink ss)))))
 392      a6loop
 393      (cond ((order< k (key pp)) (setbp pp -1) (setq pp (llink pp)))
 394            ((order< (key pp) k) (setbp pp 1) (setq pp (rlink pp)))
 395            ((order= k (key pp)) (go a7)))
 396      (go a6loop)
 397      a7   (cond ((order< k (key ss)) (go a7l)) (t (go a7r)))
 398      a7l  (cond ((order= 0. (bp ss)) (setbp ss -1) (setllink head (1+ (llink head))) (return head))
 399                 ((order= (bp ss) 1) (setbp ss 0) (return head)))
 400      (cond ((order= (bp rr) -1) nil) (t (go a9l)))
 401      (setq pp rr)
 402      (setllink ss (rlink rr))
 403      (setrlink rr ss)
 404      (setbp ss 0)
 405      (setbp rr 0)
 406      (go a10)
 407      a9l  (setq pp (rlink rr))
 408      (setrlink rr (llink pp))
 409      (setllink pp rr)
 410      (setllink ss (rlink pp))
 411      (setrlink pp ss)
 412      (cond ((order= (bp pp) -1.) (setbp ss 1.) (setbp rr 0.))
 413            ((order= (bp pp) 0.) (setbp ss 0.) (setbp rr 0.))
 414            ((order= (bp pp) 1.) (setbp ss 0.) (setbp rr -1.)))
 415      (setbp pp 0.)
 416      (go a10)
 417      a7r  (cond ((order= 0. (bp ss)) (setbp ss 1.) (setllink head (1+ (llink head))) (return head))
 418                 ((order= (bp ss) -1.) (setbp ss 0.) (return head)))
 419      (cond ((order= (bp rr) 1.) nil) (t (go a9r)))
 420      (setq pp rr)
 421      (setrlink ss (llink rr))
 422      (setllink rr ss)
 423      (setbp ss 0.)
 424      (setbp rr 0.)
 425      (go a10)
 426      a9r  (setq pp (llink rr))
 427      (setllink rr (rlink pp))
 428      (setrlink pp rr)
 429      (setrlink ss (llink pp))
 430      (setllink pp ss)
 431      (cond ((order= (bp pp) 1.) (setbp ss -1.) (setbp rr 0.))
 432            ((order= (bp pp) 0.) (setbp ss 0.) (setbp rr 0.))
 433            ((order= (bp pp) -1.) (setbp ss 0.) (setbp rr 1.)))
 434      (setbp pp 0.)
 435      a10  (cond ((eq ss (rlink tt)) (setrlink tt pp)) (t (setllink tt pp)))
 436      (return head)))
 437
 438 (defun avlinit (key rec)
 439   (cons 'top (cons (cons 0. (cons key (cons (cons nil nil) (cons 0. rec)))) (cons 0. nil))))
 440
 441
 442 ;; UNTREE CONVERTS THE TREE TO A LIST WHICH LOOKS LIKE ( SmALLEST-KEY RECORD NEXT-SMALLEST-KEY RECORD ....  LARGEST-KEY
 443 ;;RECORD)
 444
 445 (defun untree (h) (prog (ans) (untree1 (rlink h)) (return ans)))
 446
 447 (defun untree1 (h)
 448   (cond ((null h) ans)
 449         ((null (rlink h)) (setq ans (tcons3 (key h) (rec h) ans)) (untree1 (llink h)))
 450         (t (setq ans (tcons3 (key h) (rec h) (untree1 (rlink h)))) (untree1 (llink h)))))
 451
 452 (defun tcons3 (r s tt) (cond ((poispzero s) tt) (t (cons r (cons s tt)))))
 453
 454
 455 (defun poismerges (a ae l)
 456   (cond ((equal poishift ae) l)         ; SINE(0) IS 0
 457         ((poisnegpred ae) (poismerge (poisco* -1 a) (poischangesign ae) l))
 458         (t (poismerge a ae l))))
 459
 460 (defun poismergec (a ae l)
 461   (cond ((poisnegpred ae) (poismerge a (poischangesign ae) l)) (t (poismerge a ae l))))
 462
 463 (defun poismerge (a ae l) (cond ((poispzero a) nil) (t (merge11 a ae l))))
 464
 465 (defun poismerge2 (r s)
 466   (cond ((null r) s)
 467         ((null s) r)
 468         (t (prog (m n tt)
 469               (setq m (setq n (cons 0. r)))
 470               a    (cond ((null r) (rplacd m s) (return (cdr n)))
 471                          ((null s) (return (cdr n)))
 472                          ((equal (car r) (car s))
 473                           (setq tt (poisco+ (cadr r) (cadr s)))
 474                           (cond ((poispzero tt) (rplacd m (cddr r)) (setq r (cddr r) s (cddr s)))
 475                                 (t (rplaca (cdr r) tt) (setq s (cddr s) r (cddr r) m (cddr m)))))
 476                          ((> (car r) (car s))
 477                           (rplacd m s)
 478                           (setq s (cddr s))
 479                           (rplacd (cddr m) r)
 480                           (setq m (cddr m)))
 481                          (t (setq r (cddr r)) (setq m (cddr m))))
 482               (go a)))))
 483
 484 (defun merge11 (a ae l)
 485   (poismerge2 (list ae a) l))
 486
 487 (defun poismergesx (a ae l)
 488   (cond ((equal poishift ae) l)         ; SINE(0) IS 0
 489         ((poisnegpred ae) (avlinsert (poischangesign ae) (poisco* -1 a) l))
 490         (t (avlinsert ae a l))))
 491
 492 (defun poismergecx (a ae l)
 493   (cond ((poisnegpred ae) (avlinsert (poischangesign ae) a l)) (t (avlinsert ae a l))))
 494
 495 (defun poisctimes1 (c h)
 496   (cond ((null h) nil)
 497         ((and trim (trimf (car h))) (poisctimes1 c (cddr h)))
 498         (t (tcons (car h) (cons (poisco* c (cadr h)) (poisctimes1 c (cddr h)))))))
 499
 500 (defun trimf (m)
 501   (meval (list '($poistrim)
 502                (poisxcoef m '$u)
 503                (poisxcoef m '$v)
 504                (poisxcoef m '$w)
 505                (poisxcoef m '$x)
 506                (poisxcoef m '$y)
 507                (poisxcoef m '$z))))
 508
 509 (defmfun $poistimes (a b)
 510   (prog (slc clc temp ae aa zero trim t1 t2 f1 f2)
 511      (setq a (intopois a) b (intopois b))
 512      (cond ((or (getl-lm-fcn-prop '$poistrim '(expr subr))
 513                 (mget '$poistrim 'mexpr))
 514             (setq trim t)))
 515      (cond ((nonperiod a) (return ($poisctimes (cadr (caddr a)) b)))
 516            ((nonperiod b) (return ($poisctimes (cadr (caddr b)) a))))
 517      (setq slc (avlinit poishift (setq zero (intopoisco 0.))) clc (avlinit poishift zero))
 518      ;; PROCEED THROUGH ALL THE SINES IN ARGUMENT A
 519      (do ((sla
 520            (cadr a)
 521            (cddr sla)))(
 522                         (null sla))
 523        (setq aa (halve (cadr sla)) ae (car sla))
 524        ;; SINE(U)*SINE(V) ==> (-COSINE(U+V) + COSINE(U-V))/2
 525        (do ((slb
 526              (cadr b)
 527              (cddr slb)))(
 528                           (null slb))
 529          (setq t1 (+ ae poishift (- (car slb))) t2 (+ ae (- poishift) (car slb)))
 530          (cond(trim(setq f1(trimf t1) f2 (trimf t2)))
 531               (t (setq f1 nil f2 nil)))
 532          (setq temp (poisco* aa (cadr slb)))
 533          (cond ((poispzero temp) nil)
 534                (t (or f1 (poismergecx temp t1 clc))
 535                   (or f2 (poismergecx (poisco* -1 temp) t2 clc)))))
 536        ;; SINE*COSINE ==> SINE + SINE
 537        (do ((clb
 538              (caddr b)
 539              (cddr clb)))(
 540                           (null clb))
 541          (setq t1 (+ ae poishift (- (car clb))) t2 (+ ae (- poishift) (car clb)))
 542          (cond(trim(setq f1(trimf t1) f2 (trimf t2)))
 543               (t (setq f1 nil f2 nil)))
 544          (setq temp (poisco* aa (cadr clb)))
 545          (cond ((poispzero temp) nil)
 546                (t (or f1 (poismergesx temp t1 slc)) (or f2 (poismergesx temp t2 slc))))))
 547      ;; PROCEED THROUGH ALL THE COSINES IN ARGUMENT A
 548      (do ((cla
 549            (caddr a)
 550            (cddr cla)))(
 551                         (null cla))
 552        (setq aa (halve (cadr cla)) ae (car cla))
 553        ;; COSINE*SINE ==> SINE - SINE
 554        (do ((slb
 555              (cadr b)
 556              (cddr slb)))(
 557                           (null slb))
 558          (setq t1 (+ ae poishift (- (car slb)))
 559                t2 (+ ae (- poishift) (car slb)))
 560          (cond (trim (setq f1 (trimf t1) f2 (trimf t2)))
 561                (t (setq f1 nil f2 nil)))
 562          (cond (t (setq temp (poisco* aa (cadr slb)))
 563                   (cond ((poispzero temp) nil)
 564                         (t (or f1 (poismergesx (poisco* -1 temp) t1 slc))
 565                            (or f2 (poismergesx temp t2 slc)))))))
 566        ;; COSINE*COSINE ==> COSINE + COSINE
 567        (do ((clb (caddr b) (cddr clb)))
 568            ((null clb))
 569          (setq t1 (+ ae poishift (- (car clb)))
 570                t2 (+ ae (- poishift) (car clb)))
 571          (cond (trim (setq f1 (trimf t1) f2 (trimf t2)))
 572                (t (setq f1 nil f2 nil)))
 573          (cond
 574            (t (setq temp (poisco* aa (cadr clb)))
 575               (cond ((poispzero temp) nil)
 576                     (t (or f1 (poismergecx temp t1 clc))
 577                        (or f2 (poismergecx temp t2 clc))))))))
 578      (return (list '(mpois simp) (untree slc) (untree clc)))))
 579
 580 (defmfun $poisexpt (p n)
 581   (prog (u h)
 582      (cond ((oddp n) (setq u p)) (t (setq u (setq h (intopois 1.)))))
 583      a    (setq n (ash n -1))
 584      (cond ((zerop n) (return u)))
 585      (setq p ($poistimes p p))
 586      (cond ((oddp n) (setq u (cond ((equal u h) p) (t ($poistimes u p))))))
 587      (go a)))
 588
 589 (defmfun $poissquare (a) ($poisexpt a 2))
 590
 591 ;;; $POISINT INTEGRATES A POISSON SERIES WRT X,Y, Z, U, V, W.  THE VARIABLE OF
 592 ;;; INTEGRATION MUST OCCUR ONLY IN THE ARGUMENTS OF SIN OR COS,
 593 ;;; OR ONLY IN THE COEFFICIENTS.  POISCOINTEG IS CALLED TO INTEGRATE COEFFS.
 594
 595 ;;; NON-PERIODIC TERMS ARE REMOVED.
 596
 597 (defmfun $poisint (p m)
 598   (declare (special m))
 599   (prog (b*)
 600      (setq p (intopois p))
 601      (cond ((member m '($u $v $w $x $y $z) :test #'eq)
 602             (return (list (car p)
 603                           (cosint* (caddr p) m)
 604                           (sinint* (cadr p) m))))
 605            (t (return (list (car p)
 606                             (poisint4 (cadr p))
 607                             (poisint4 (caddr p))))))))
 608
 609 (defun poisint4 (y)
 610   (declare (special m))
 611   (cond ((null y) nil)
 612         (t (tcons3 (car y)(poiscointeg (cadr y) m) (poisint4 (cddr y))))))
 613
 614 ;;;COSINT* INTEGRATES COSINES TO GET SINES
 615
 616 (defun cosint* (h m)
 617   (cond ((null h) nil)
 618         ((equal 0 (setq b* (poisxcoef (car h) m))) (cosint* (cddr h) m))
 619         (t (tcons (car h)
 620                   (cons (poisco* (intopoisco (list '(mexpt) b* -1)) (cadr h))
 621                         (cosint* (cddr h) m))))))
 622
 623 (defun sinint* (h m)
 624   (cond ((null h) nil)
 625         ((equal 0 (setq b* (poisxcoef (car h) m))) (sinint* (cddr h) m))
 626         (t (tcons (car h)
 627                   (cons (poisco* (intopoisco (list '(mexpt) (- (poisxcoef (car h) m)) -1))
 628                                  (cadr h))
 629                         (sinint* (cddr h) m))))))
 630
 631
 632 ;;; $POISSUBST SUBSTITUTES AN EXPRESSION FOR A VARIABLE IN ARGUMENT OF TRIG FUNCTIONS OR
 633 ;;; COEFFICIENTS.
 634
 635 (defun poissubsta (a b* c)
 636   (prog (ss cc)
 637      (setq h* (- (poisencode (list '(mplus) a (list '(mtimes) -1 b*))) poishift))
 638      (poissubst1s (cadr c))
 639      (poissubst1c (caddr c))
 640      (return (list (car c) ss cc))))
 641
 642 (defun poissubst1s (c)
 643   (cond ((null c) nil)
 644         (t (setq ss (poismerges (cadr c) (argsubst (car c)) ss))
 645            (poissubst1s (cddr c)))))
 646
 647 (defun poissubst1c (c)
 648   (cond ((null c) nil)
 649         (t (setq cc (poismergec (cadr c) (argsubst (car c)) cc))
 650            (poissubst1c (cddr c)))))
 651
 652 (defun argsubst (c)
 653   (+ c (* h* (poisxcoef c b*))))
 654
 655 (defmfun $poissubst (aa bb cc &optional dd nn)
 656   (if (and dd nn)
 657       (fancypoissubst aa bb (intopois cc) (intopois dd) nn)
 658       (let ((a* aa) (b* bb) (c (intopois cc)))
 659         (if (member b* '($u $v $w $x $y $z) :test #'eq)
 660             (poissubsta a* b* c)
 661             (list (car c) (poissubstco1 (cadr c)) (poissubstco1 (caddr c)))))))
 662
 663 (declare-top (unspecial $u $v $w $x $y $z))
 664
 665 (defun poissubstco1 (c)
 666   (if (null c)
 667       nil
 668       (tcons (car c) (cons (poissubstco a* b* (cadr c)) (poissubstco1 (cddr c))))))
 669
 670 (declare-top (special dc ds *ans))
 671
 672 (defun fancypoissubst (a b* c d n)
 673   ;;SUBSTITUTES A+D FOR B IN C, WHERE D IS EXPANDED IN POWERSERIES TO ORDER N
 674   (prog (h* dc ds *ans)
 675      (setq *ans (list '(mpois simp) nil nil) d (intopois d) dc (intopois 1) ds (intopois 0))
 676      (when (equal n 0) (return ($poissubst a b* c)))
 677      (fancypois1s d 1 1 n)
 678      (setq h* (- (poisencode (list '(mplus) a (list '(mtimes) -1 b*))) poishift))
 679      (fancypas (cadr c))
 680      (fancypac (caddr c))
 681      (return *ans)))
 682
 683 (defun fancypois1s (d dp n lim) ; DP IS LAST POWER: D^(N-1), LIM IS HIGHEST TO
 684   (cond ((> n lim) nil)         ;GO
 685         (t (setq ds ($poisplus ds
 686                                ($poisctimes (list '(rat)
 687                                                   (expt -1 (ash (1- n) -1))
 688                                                   (factorial n))
 689                                             (setq dp ($poistimes dp d)))))
 690            (fancypois1c d dp (1+ n) lim))))
 691
 692 (defun fancypois1c (d dp n lim) ; DP IS LAST POWER: D^(N-1), LIM IS HIGHEST TO
 693   (cond ((> n lim) nil)         ;GO
 694         (t (setq dc
 695                  ($poisplus dc
 696                             ($poisctimes (list '(rat) (expt -1 (ash n -1)) (factorial n))
 697                                          (setq dp ($poistimes dp d)))))
 698            (fancypois1s d dp (1+ n) lim))))
 699
 700 ;;; COS(R+K*B) ==> K*COS(R+K*A)*DC - K*SIN(R+K*A)*DS
 701 ;;; SIN(R+K*B) ==> K*COS(R+K*A)*DS + K*SIN(R+K*A)*DC
 702
 703 (defun fancypac (c)
 704   (prog nil
 705      (cond ((null c) (return nil)))
 706      (setq *coef (poisxcoef (car c) b*))
 707      (cond ((equal *coef 0)
 708             (setq *ans ($poisplus *ans (list '(mpois simp) nil (list (car c) (cadr c)))))
 709             (go end)))
 710      (cond ((poispzero (setq *coef (poisco* (cadr c) (intopoisco *coef)))) (go end)))
 711      (setq *argc (argsubst (car c)))
 712      (setq *ans
 713            ($poisplus *ans
 714                       ($poisplus ($poistimes (list '(mpois simp)
 715                                                    nil
 716                                                    (poismergec *coef *argc nil))
 717                                              dc)
 718                                  ($poistimes (list '(mpois simp)
 719                                                    (poismerges (poisco* -1 *coef) *argc nil)
 720                                                    nil)
 721                                              ds))))
 722      end  (fancypac (cddr c))))
 723
 724 (defun fancypas (c)
 725   (prog nil
 726      (cond ((null c) (return nil)))
 727      (setq *coef (poisxcoef (car c) b*))
 728      (cond ((equal *coef 0.)
 729             (setq *ans ($poisplus *ans (list '(mpois simp) (list (car c) (cadr c)) nil)))
 730             (go end)))
 731      (cond ((poispzero (setq *coef (poisco* (cadr c) (intopoisco *coef)))) (go end)))
 732      (setq *argc (argsubst (car c)))
 733      (setq *ans ($poisplus *ans
 734                            ($poisplus ($poistimes (list '(mpois simp)
 735                                                         nil
 736                                                         (poismergec *coef *argc nil))
 737                                                   ds)
 738                                       ($poistimes (list '(mpois simp)
 739                                                         (poismerges *coef *argc nil)
 740                                                         nil)
 741                                                   dc))))
 742      end  (fancypas (cddr c))))