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)))) 743