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
 277      ;; DO COSINES
 278      (do ((m
 279            (caddr p)
 280            (cddr m)))(
 281                       (null m))
 282        (displa (simplifya (list '(mtimes)
 283                                 (poiscdecode (cadr m))
 284                                 (cond ((equal (car m) poishift) 1.)
 285                                       (t (list '(%cos) (poisdecodec (car m))))))
 286                           t))
 287        (terpri))
 288      (return '$done)))
 289
 290
 291 ;;; $POISDIFF DIFFERENTIATES A POISSON SERIES WRT X, Y, Z, U, V, W, OR A COEFF VAR.
 292
 293
 294 (defmfun $poisdiff (p m)
 295   (declare (special m))
 296   (cond ((member m '($u $v $w $x $y $z) :test #'eq)
 297          (list (car p) (cosdif (caddr p) m) (sindif (cadr p) m)))
 298         (t (list (car p) (poisdif4(cadr p)) (poisdif4 (caddr p))))))
 299
 300
 301 (defun poisdif4 (y)
 302   (declare (special m))
 303   (cond ((null y) nil)
 304         (t (tcons3 (car y)(poiscodif (cadr y) m) (poisdif4 (cddr y))))))
 305
 306
 307 ;;; COSDIF DIFFERENTIATES COSINES TO GET SINES
 308
 309 (defun cosdif (h m)
 310   (cond ((null h) nil)
 311         (t (tcons (car h)
 312                   (cons (poisco* (intopoisco (- (poisxcoef (car h) m))) (cadr h))
 313                         (cosdif (cddr h) m))))))
 314
 315 (defun sindif (h m)
 316   (cond ((null h) nil)
 317         (t (tcons (car h)
 318                   (cons (poisco* (intopoisco (poisxcoef (car h) m)) (cadr h)) (sindif (cddr h) m))))))
 319
 320 (defun poisxcoef (h m)
 321   (- (rem (quotient h
 322                     (expt poists
 323                           (cadr (member m '($u 0 $v 1 $w 2 $x 3 $y 4 $z 5) :test #'eq))))
 324           poists)
 325      poistsm))
 326
 327
 328 ;;; AVL BALANCED TREE SEARCH AND INSERTION.
 329 ;;; NODE LOOKS LIKE (KEY (LLINK .  RLKINK) BALANCEFACTOR .  RECORD)
 330 ;;; PROGRAM FOLLOWS ALGORITHM GIVEN IN KNUTH VOL. 3 455-57
 331
 332 (declare-top (special ans))
 333
 334
 335 ;; MACROS TO EXTRACT FIELDS FROM NODE
 336
 337 (defmacro key  (&rest l) (cons 'car l))
 338
 339 (defmacro llink  (&rest l) (cons 'caadr l))
 340
 341 (defmacro rlink  (&rest l) (cons 'cdadr l))
 342
 343 (defmacro bp  (&rest l) (cons 'caddr l))
 344
 345 (defmacro rec  (&rest l) (cons 'cdddr l))
 346
 347
 348 ;;  FOR ORDERING KEYS
 349
 350 (defmacro order<  (&rest l) (cons '<  l))
 351 (defmacro order=  (&rest l) (cons '=  l))
 352
 353 ;; MACROS TO SET FIELDS IN NODE
 354
 355 (defmacro setrlink  (&rest l) (setq l (cons nil l))
 356           (list 'rplacd (list 'cadr (cadr l)) (caddr l)))
 357
 358 (defmacro setllink  (&rest l) (setq l (cons nil l))
 359           (list 'rplaca (list 'cadr (cadr l)) (caddr l)))
 360
 361 (defmacro setbp  (&rest l) (setq l (cons nil l))
 362           (list 'rplaca (list 'cddr (cadr l)) (caddr l)))
 363
 364 (defmacro setrec  (&rest l)(setq l (cons nil l))
 365           (list 'rplacd (list 'cddr (cadr l)) (caddr l)))
 366
 367
 368 (defun insert-it (pp newrec) (setrec pp (poisco+ (rec pp) newrec)))
 369
 370 (defun avlinsert (k newrec head)
 371   (prog (qq tt ss pp rr)
 372      (setq tt head)
 373      (setq ss (setq pp (rlink head)))
 374      a2   (cond ((order< k (key pp)) (go a3))
 375                 ((order< (key pp) k) (go a4))
 376                 (t (insert-it pp newrec) (return head)))
 377      a3   (setq qq (llink pp))
 378      (cond ((null qq) (setllink pp (cons k (cons (cons nil nil) (cons 0. newrec)))) (go a6))
 379            ((order= 0. (bp qq)) nil)
 380            (t (setq tt pp ss qq)))
 381      (setq pp qq)
 382      (go a2)
 383      a4   (setq qq (rlink pp))
 384      (cond ((null qq) (setrlink pp (cons k (cons (cons nil nil) (cons 0. newrec)))) (go a6))
 385            ((order= 0 (bp qq)) nil)
 386            (t (setq tt pp ss qq)))
 387      (setq pp qq)
 388      (go a2)
 389      a6   (cond ((order< k (key ss)) (setq rr (setq pp (llink ss)))) (t (setq rr (setq pp (rlink ss)))))
 390      a6loop
 391      (cond ((order< k (key pp)) (setbp pp -1) (setq pp (llink pp)))
 392            ((order< (key pp) k) (setbp pp 1) (setq pp (rlink pp)))
 393            ((order= k (key pp)) (go a7)))
 394      (go a6loop)
 395      a7   (cond ((order< k (key ss)) (go a7l)) (t (go a7r)))
 396      a7l  (cond ((order= 0. (bp ss)) (setbp ss -1) (setllink head (1+ (llink head))) (return head))
 397                 ((order= (bp ss) 1) (setbp ss 0) (return head)))
 398      (cond ((order= (bp rr) -1) nil) (t (go a9l)))
 399      (setq pp rr)
 400      (setllink ss (rlink rr))
 401      (setrlink rr ss)
 402      (setbp ss 0)
 403      (setbp rr 0)
 404      (go a10)
 405      a9l  (setq pp (rlink rr))
 406      (setrlink rr (llink pp))
 407      (setllink pp rr)
 408      (setllink ss (rlink pp))
 409      (setrlink pp ss)
 410      (cond ((order= (bp pp) -1.) (setbp ss 1.) (setbp rr 0.))
 411            ((order= (bp pp) 0.) (setbp ss 0.) (setbp rr 0.))
 412            ((order= (bp pp) 1.) (setbp ss 0.) (setbp rr -1.)))
 413      (setbp pp 0.)
 414      (go a10)
 415      a7r  (cond ((order= 0. (bp ss)) (setbp ss 1.) (setllink head (1+ (llink head))) (return head))
 416                 ((order= (bp ss) -1.) (setbp ss 0.) (return head)))
 417      (cond ((order= (bp rr) 1.) nil) (t (go a9r)))
 418      (setq pp rr)
 419      (setrlink ss (llink rr))
 420      (setllink rr ss)
 421      (setbp ss 0.)
 422      (setbp rr 0.)
 423      (go a10)
 424      a9r  (setq pp (llink rr))
 425      (setllink rr (rlink pp))
 426      (setrlink pp rr)
 427      (setrlink ss (llink pp))
 428      (setllink pp ss)
 429      (cond ((order= (bp pp) 1.) (setbp ss -1.) (setbp rr 0.))
 430            ((order= (bp pp) 0.) (setbp ss 0.) (setbp rr 0.))
 431            ((order= (bp pp) -1.) (setbp ss 0.) (setbp rr 1.)))
 432      (setbp pp 0.)
 433      a10  (cond ((eq ss (rlink tt)) (setrlink tt pp)) (t (setllink tt pp)))
 434      (return head)))
 435
 436 (defun avlinit (key rec)
 437   (cons 'top (cons (cons 0. (cons key (cons (cons nil nil) (cons 0. rec)))) (cons 0. nil))))
 438
 439
 440 ;; UNTREE CONVERTS THE TREE TO A LIST WHICH LOOKS LIKE ( SmALLEST-KEY RECORD NEXT-SMALLEST-KEY RECORD ....  LARGEST-KEY
 441 ;;RECORD)
 442
 443 (defun untree (h) (prog (ans) (untree1 (rlink h)) (return ans)))
 444
 445 (defun untree1 (h)
 446   (cond ((null h) ans)
 447         ((null (rlink h)) (setq ans (tcons3 (key h) (rec h) ans)) (untree1 (llink h)))
 448         (t (setq ans (tcons3 (key h) (rec h) (untree1 (rlink h)))) (untree1 (llink h)))))
 449
 450 (defun tcons3 (r s tt) (cond ((poispzero s) tt) (t (cons r (cons s tt)))))
 451
 452
 453 (defun poismerges (a ae l)
 454   (cond ((equal poishift ae) l)         ; SINE(0) IS 0
 455         ((poisnegpred ae) (poismerge (poisco* -1 a) (poischangesign ae) l))
 456         (t (poismerge a ae l))))
 457
 458 (defun poismergec (a ae l)
 459   (cond ((poisnegpred ae) (poismerge a (poischangesign ae) l)) (t (poismerge a ae l))))
 460
 461 (defun poismerge (a ae l) (cond ((poispzero a) nil) (t (merge11 a ae l))))
 462
 463 (defun poismerge2 (r s)
 464   (cond ((null r) s)
 465         ((null s) r)
 466         (t (prog (m n tt)
 467               (setq m (setq n (cons 0. r)))
 468               a    (cond ((null r) (rplacd m s) (return (cdr n)))
 469                          ((null s) (return (cdr n)))
 470                          ((equal (car r) (car s))
 471                           (setq tt (poisco+ (cadr r) (cadr s)))
 472                           (cond ((poispzero tt) (rplacd m (cddr r)) (setq r (cddr r) s (cddr s)))
 473                                 (t (rplaca (cdr r) tt) (setq s (cddr s) r (cddr r) m (cddr m)))))
 474                          ((> (car r) (car s))
 475                           (rplacd m s)
 476                           (setq s (cddr s))
 477                           (rplacd (cddr m) r)
 478                           (setq m (cddr m)))
 479                          (t (setq r (cddr r)) (setq m (cddr m))))
 480               (go a)))))
 481
 482 (defun merge11 (a ae l)
 483   (poismerge2 (list ae a) l))
 484
 485 (defun poismergesx (a ae l)
 486   (cond ((equal poishift ae) l)         ; SINE(0) IS 0
 487         ((poisnegpred ae) (avlinsert (poischangesign ae) (poisco* -1 a) l))
 488         (t (avlinsert ae a l))))
 489
 490 (defun poismergecx (a ae l)
 491   (cond ((poisnegpred ae) (avlinsert (poischangesign ae) a l)) (t (avlinsert ae a l))))
 492
 493 (defun poisctimes1 (c h)
 494   (cond ((null h) nil)
 495         ((and trim (trimf (car h))) (poisctimes1 c (cddr h)))
 496         (t (tcons (car h) (cons (poisco* c (cadr h)) (poisctimes1 c (cddr h)))))))
 497
 498 (defun trimf (m)
 499   (meval (list '($poistrim)
 500                (poisxcoef m '$u)
 501                (poisxcoef m '$v)
 502                (poisxcoef m '$w)
 503                (poisxcoef m '$x)
 504                (poisxcoef m '$y)
 505                (poisxcoef m '$z))))
 506
 507 (defmfun $poistimes (a b)
 508   (prog (slc clc temp ae aa zero trim t1 t2 f1 f2)
 509      (setq a (intopois a) b (intopois b))
 510      (cond ((or (getl-lm-fcn-prop '$poistrim '(expr subr))
 511                 (mget '$poistrim 'mexpr))
 512             (setq trim t)))
 513      (cond ((nonperiod a) (return ($poisctimes (cadr (caddr a)) b)))
 514            ((nonperiod b) (return ($poisctimes (cadr (caddr b)) a))))
 515      (setq slc (avlinit poishift (setq zero (intopoisco 0.))) clc (avlinit poishift zero))
 516      ;; PROCEED THROUGH ALL THE SINES IN ARGUMENT A
 517      (do ((sla
 518            (cadr a)
 519            (cddr sla)))(
 520                         (null sla))
 521        (setq aa (halve (cadr sla)) ae (car sla))
 522        ;; SINE(U)*SINE(V) ==> (-COSINE(U+V) + COSINE(U-V))/2
 523        (do ((slb
 524              (cadr b)
 525              (cddr slb)))(
 526                           (null slb))
 527          (setq t1 (+ ae poishift (- (car slb))) t2 (+ ae (- poishift) (car slb)))
 528          (cond(trim(setq f1(trimf t1) f2 (trimf t2)))
 529               (t (setq f1 nil f2 nil)))
 530          (setq temp (poisco* aa (cadr slb)))
 531          (cond ((poispzero temp) nil)
 532                (t (or f1 (poismergecx temp t1 clc))
 533                   (or f2 (poismergecx (poisco* -1 temp) t2 clc)))))
 534        ;; SINE*COSINE ==> SINE + SINE
 535        (do ((clb
 536              (caddr b)
 537              (cddr clb)))(
 538                           (null clb))
 539          (setq t1 (+ ae poishift (- (car clb))) t2 (+ ae (- poishift) (car clb)))
 540          (cond(trim(setq f1(trimf t1) f2 (trimf t2)))
 541               (t (setq f1 nil f2 nil)))
 542          (setq temp (poisco* aa (cadr clb)))
 543          (cond ((poispzero temp) nil)
 544                (t (or f1 (poismergesx temp t1 slc)) (or f2 (poismergesx temp t2 slc))))))
 545      ;; PROCEED THROUGH ALL THE COSINES IN ARGUMENT A
 546      (do ((cla
 547            (caddr a)
 548            (cddr cla)))(
 549                         (null cla))
 550        (setq aa (halve (cadr cla)) ae (car cla))
 551        ;; COSINE*SINE ==> SINE - SINE
 552        (do ((slb
 553              (cadr b)
 554              (cddr slb)))(
 555                           (null slb))
 556          (setq t1 (+ ae poishift (- (car slb)))
 557                t2 (+ ae (- poishift) (car slb)))
 558          (cond (trim (setq f1 (trimf t1) f2 (trimf t2)))
 559                (t (setq f1 nil f2 nil)))
 560          (cond (t (setq temp (poisco* aa (cadr slb)))
 561                   (cond ((poispzero temp) nil)
 562                         (t (or f1 (poismergesx (poisco* -1 temp) t1 slc))
 563                            (or f2 (poismergesx temp t2 slc)))))))
 564        ;; COSINE*COSINE ==> COSINE + COSINE
 565        (do ((clb (caddr b) (cddr clb)))
 566            ((null clb))
 567          (setq t1 (+ ae poishift (- (car clb)))
 568                t2 (+ ae (- poishift) (car clb)))
 569          (cond (trim (setq f1 (trimf t1) f2 (trimf t2)))
 570                (t (setq f1 nil f2 nil)))
 571          (cond
 572            (t (setq temp (poisco* aa (cadr clb)))
 573               (cond ((poispzero temp) nil)
 574                     (t (or f1 (poismergecx temp t1 clc))
 575                        (or f2 (poismergecx temp t2 clc))))))))
 576      (return (list '(mpois simp) (untree slc) (untree clc)))))
 577
 578 (defmfun $poisexpt (p n)
 579   (prog (u h)
 580      (cond ((oddp n) (setq u p)) (t (setq u (setq h (intopois 1.)))))
 581      a    (setq n (ash n -1))
 582      (cond ((zerop n) (return u)))
 583      (setq p ($poistimes p p))
 584      (cond ((oddp n) (setq u (cond ((equal u h) p) (t ($poistimes u p))))))
 585      (go a)))
 586
 587 (defmfun $poissquare (a) ($poisexpt a 2))
 588
 589 ;;; $POISINT INTEGRATES A POISSON SERIES WRT X,Y, Z, U, V, W.  THE VARIABLE OF
 590 ;;; INTEGRATION MUST OCCUR ONLY IN THE ARGUMENTS OF SIN OR COS,
 591 ;;; OR ONLY IN THE COEFFICIENTS.  POISCOINTEG IS CALLED TO INTEGRATE COEFFS.
 592
 593 ;;; NON-PERIODIC TERMS ARE REMOVED.
 594
 595 (defmfun $poisint (p m)
 596   (declare (special m))
 597   (prog (b*)
 598      (setq p (intopois p))
 599      (cond ((member m '($u $v $w $x $y $z) :test #'eq)
 600             (return (list (car p)
 601                           (cosint* (caddr p) m)
 602                           (sinint* (cadr p) m))))
 603            (t (return (list (car p)
 604                             (poisint4 (cadr p))
 605                             (poisint4 (caddr p))))))))
 606
 607 (defun poisint4 (y)
 608   (declare (special m))
 609   (cond ((null y) nil)
 610         (t (tcons3 (car y)(poiscointeg (cadr y) m) (poisint4 (cddr y))))))
 611
 612 ;;;COSINT* INTEGRATES COSINES TO GET SINES
 613
 614 (defun cosint* (h m)
 615   (cond ((null h) nil)
 616         ((equal 0 (setq b* (poisxcoef (car h) m))) (cosint* (cddr h) m))
 617         (t (tcons (car h)
 618                   (cons (poisco* (intopoisco (list '(mexpt) b* -1)) (cadr h))
 619                         (cosint* (cddr h) m))))))
 620
 621 (defun sinint* (h m)
 622   (cond ((null h) nil)
 623         ((equal 0 (setq b* (poisxcoef (car h) m))) (sinint* (cddr h) m))
 624         (t (tcons (car h)
 625                   (cons (poisco* (intopoisco (list '(mexpt) (- (poisxcoef (car h) m)) -1))
 626                                  (cadr h))
 627                         (sinint* (cddr h) m))))))
 628
 629
 630 ;;; $POISSUBST SUBSTITUTES AN EXPRESSION FOR A VARIABLE IN ARGUMENT OF TRIG FUNCTIONS OR
 631 ;;; COEFFICIENTS.
 632
 633 (defun poissubsta (a b* c)
 634   (prog (ss cc)
 635      (setq h* (- (poisencode (list '(mplus) a (list '(mtimes) -1 b*))) poishift))
 636      (poissubst1s (cadr c))
 637      (poissubst1c (caddr c))
 638      (return (list (car c) ss cc))))
 639
 640 (defun poissubst1s (c)
 641   (cond ((null c) nil)
 642         (t (setq ss (poismerges (cadr c) (argsubst (car c)) ss))
 643            (poissubst1s (cddr c)))))
 644
 645 (defun poissubst1c (c)
 646   (cond ((null c) nil)
 647         (t (setq cc (poismergec (cadr c) (argsubst (car c)) cc))
 648            (poissubst1c (cddr c)))))
 649
 650 (defun argsubst (c)
 651   (+ c (* h* (poisxcoef c b*))))
 652
 653 (defmfun $poissubst (aa bb cc &optional dd nn)
 654   (if (and dd nn)
 655       (fancypoissubst aa bb (intopois cc) (intopois dd) nn)
 656       (let ((a* aa) (b* bb) (c (intopois cc)))
 657         (if (member b* '($u $v $w $x $y $z) :test #'eq)
 658             (poissubsta a* b* c)
 659             (list (car c) (poissubstco1 (cadr c)) (poissubstco1 (caddr c)))))))
 660
 661 (declare-top (unspecial $u $v $w $x $y $z))
 662
 663 (defun poissubstco1 (c)
 664   (if (null c)
 665       nil
 666       (tcons (car c) (cons (poissubstco a* b* (cadr c)) (poissubstco1 (cddr c))))))
 667
 668 (declare-top (special dc ds *ans))
 669
 670 (defun fancypoissubst (a b* c d n)
 671   ;;SUBSTITUTES A+D FOR B IN C, WHERE D IS EXPANDED IN POWERSERIES TO ORDER N
 672   (prog (h* dc ds *ans)
 673      (setq *ans (list '(mpois simp) nil nil) d (intopois d) dc (intopois 1) ds (intopois 0))
 674      (when (equal n 0) (return ($poissubst a b* c)))
 675      (fancypois1s d 1 1 n)
 676      (setq h* (- (poisencode (list '(mplus) a (list '(mtimes) -1 b*))) poishift))
 677      (fancypas (cadr c))
 678      (fancypac (caddr c))
 679      (return *ans)))
 680
 681 (defun fancypois1s (d dp n lim) ; DP IS LAST POWER: D^(N-1), LIM IS HIGHEST TO
 682   (cond ((> n lim) nil)         ;GO
 683         (t (setq ds ($poisplus ds
 684                                ($poisctimes (list '(rat)
 685                                                   (expt -1 (ash (1- n) -1))
 686                                                   (factorial n))
 687                                             (setq dp ($poistimes dp d)))))
 688            (fancypois1c d dp (1+ n) lim))))
 689
 690 (defun fancypois1c (d dp n lim) ; DP IS LAST POWER: D^(N-1), LIM IS HIGHEST TO
 691   (cond ((> n lim) nil)         ;GO
 692         (t (setq dc
 693                  ($poisplus dc
 694                             ($poisctimes (list '(rat) (expt -1 (ash n -1)) (factorial n))
 695                                          (setq dp ($poistimes dp d)))))
 696            (fancypois1s d dp (1+ n) lim))))
 697
 698 ;;; COS(R+K*B) ==> K*COS(R+K*A)*DC - K*SIN(R+K*A)*DS
 699 ;;; SIN(R+K*B) ==> K*COS(R+K*A)*DS + K*SIN(R+K*A)*DC
 700
 701 (defun fancypac (c)
 702   (prog nil
 703      (cond ((null c) (return nil)))
 704      (setq *coef (poisxcoef (car c) b*))
 705      (cond ((equal *coef 0)
 706             (setq *ans ($poisplus *ans (list '(mpois simp) nil (list (car c) (cadr c)))))
 707             (go end)))
 708      (cond ((poispzero (setq *coef (poisco* (cadr c) (intopoisco *coef)))) (go end)))
 709      (setq *argc (argsubst (car c)))
 710      (setq *ans
 711            ($poisplus *ans
 712                       ($poisplus ($poistimes (list '(mpois simp)
 713                                                    nil
 714                                                    (poismergec *coef *argc nil))
 715                                              dc)
 716                                  ($poistimes (list '(mpois simp)
 717                                                    (poismerges (poisco* -1 *coef) *argc nil)
 718                                                    nil)
 719                                              ds))))
 720      end  (fancypac (cddr c))))
 721
 722 (defun fancypas (c)
 723   (prog nil
 724      (cond ((null c) (return nil)))
 725      (setq *coef (poisxcoef (car c) b*))
 726      (cond ((equal *coef 0.)
 727             (setq *ans ($poisplus *ans (list '(mpois simp) (list (car c) (cadr c)) nil)))
 728             (go end)))
 729      (cond ((poispzero (setq *coef (poisco* (cadr c) (intopoisco *coef)))) (go end)))
 730      (setq *argc (argsubst (car c)))
 731      (setq *ans ($poisplus *ans
 732                            ($poisplus ($poistimes (list '(mpois simp)
 733                                                         nil
 734                                                         (poismergec *coef *argc nil))
 735                                                   ds)
 736                                       ($poistimes (list '(mpois simp)
 737                                                         (poismerges *coef *argc nil)
 738                                                         nil)
 739                                                   dc))))
 740      end  (fancypas (cddr c))))