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 1982 Massachusetts Institute of Technology ;;; 9;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 10 11(in-package :maxima) 12 13(macsyma-module ratout) 14 15;; THIS IS THE OUT-OF-CORE SEGMENT OF THE RATIONAL FUNCTION PACKAGE. 16 17(declare-top (special $algebraic varlist ss *y* f $factorflag modulus 18 genvar *alpha *x* *p *max *var *res *chk *l 19 $ratfac u* $ratwtlvl *ratweights $ratweights)) 20 21(declare-top (special xv bigf1 bigf2 22 gcd $factorflag)) 23 24;; NEWGCD (X,Y) RETURNS A LIST OF THREE ITEMS, 25;; (GCD, X/GCD, Y/GCD) 26 27(defun newgcd (x y modulus) 28 (set-modulus modulus) 29 (let ((a (cond ((pcoefp x) 30 (cond ((zerop x) y) 31 ((pcoefp y) (cgcd x y)) 32 (t (pcontent1 (cdr y) x)))) 33 ((pcoefp y) (cond ((zerop y) x) (t (pcontent1 (cdr x) y)))) 34 ((pointergp (p-var x) (p-var y)) (oldcontent1 (cdr x) y)) 35 ((pointergp (p-var y) (p-var x)) (oldcontent1 (cdr y) x)) 36 (t nil)))) 37 (cond (a (list a (pquotient x a) (pquotient y a))) 38 (modulus (pgcdp x y modulus)) 39 (t (pgcdm x y))))) 40 41;;;*** PMODCONTENT COMPUTES CONTENT OF 42;;; P IN 43;; Z [X ] [X , X , ..., X ] 44;; P V 1 2 V-1 45 46;; PMODCONTENT OF 3*A*X IS A, IF MAINVAR IS X (=X ) 47;; V 48 49(defun pmodcontent (p) 50 (prog (*var *chk *res *max gcd) 51 (setq *chk (car p)) 52 (setq *max 0) 53 (setq *var (pnext (cdr p) nil)) 54 (cond ((pointergp xv *chk) (go ret1)) 55 ((null *var) (return (list p 1)))) 56 (pgath1 (cdr p)) 57 a (setq *res 0) 58 (pgath3 (cdr p)) 59 a2 (cond ((pcoefp *res) (cond ((pzerop *res) nil)(t(go ret1)))) 60 ((not (eq (car *res) *chk)) (go ret1)) 61 ((not (univar (cdr *res))) 62 (setq *res (car (pmodcontent *res))) 63 (go a2)) 64 (gcd (setq gcd (pgcdu gcd *res))) 65 (t (setq gcd *res))) 66 (cond ((pcoefp gcd) (go ret1)) 67 ((minusp (setq *max (1- *max))) 68 (return (list gcd (pquotient p gcd))))) 69 (go a) 70 ret1 (return (list 1 p)))) 71 72(defun pgathercoef (p *chk *res) 73 (if (not (eq (car p) *chk)) 1 (pgath2 (cdr p) nil))) 74 75(defun pgath1 (p) 76 (prog nil 77 (cond ((null p) (return *max)) 78 ((pcoefp (cadr p)) nil) 79 ((eq (caadr p) *var) (setq *max (max *max (cadadr p))))) 80 (return (pgath1 (cddr p))))) 81 82(defun pgath2 (p vmax) 83 (prog (v2) 84 (cond ((null p) (return *res)) 85 ((pcoefp (cadr p)) nil) 86 ((vgreat (setq v2 (pdegreer (cadr p))) vmax) 87 (setq *res (psimp *chk 88 (list (car p) (leadcoefficient (cadr p))))) 89 (setq vmax v2)) 90 ((equal vmax v2) 91 (setq *res 92 (pplus *res 93 (psimp *chk 94 (list (car p) (leadcoefficient (cadr p)))))))) 95 (return (pgath2 (cddr p) vmax)))) 96 97(defun pgath3 (p) 98 (prog (zz) 99 (cond ((null p) (return *res)) 100 ((pcoefp (cadr p)) 101 (cond ((equal *max 0) (setq zz (cadr p)) (go add)) (t (go ret)))) 102 ((eq (caadr p) *var) (setq zz (ptterm (cdadr p) *max)) (go add))) 103 (cond ((equal *max 0) (setq zz (cadr p))) (t (go ret))) 104 add (cond ((equal zz 0) (go ret))) 105 (setq *res (pplus *res (psimp *chk (list (car p) zz)))) 106 ret (return (pgath3 (cddr p))))) 107 108(defun pnext (x *l) 109 (pnext1 x) 110 (cond ((null *l) nil) 111 (t (car (sort *l #'pointergp))))) 112 113(defun pnext1 (x) 114 (prog nil 115 (cond ((null x) (return *l)) 116 ((or (pcoefp (cadr x)) (member (caadr x) *l :test #'eq)) nil) 117 (t (setq *l (cons (caadr x) *l)))) 118 (return (pnext1 (cddr x))))) 119 120(defun vgreat (x y) 121 (cond ((null x) nil) 122 ((null y) t) 123 ((pointergp (car x)(car y))t) 124 ((not (eq (car x)(car y)))nil) 125 ((> (cadr x)(cadr y)) t) 126 ((equal (cadr x)(cadr y))(vgreat (cddr x)(cddr y))) 127 (t nil))) 128 129(defun pdegreer (x) 130 (if (pcoefp x) () (cons (car x) (cons (cadr x) (pdegreer (caddr x)))))) 131 132;;*** PGCDP CORRESPONDS TO BROWN'S ALGORITHM P 133 134(defun pgcdp (bigf1 bigf2 modulus) 135 (prog (c c1 c2 n q 136 h1tilde h2tilde gstar h1star 137 h2star xv e b 138 gbar nubar nu1bar nu2bar 139 gtilde f1tilde f2tilde biggtilde 140 degree f1 f1f2) 141 (set-modulus modulus) 142 (cond ((and (univar (cdr bigf1)) (univar (cdr bigf2))) 143 (setq q (pgcdu bigf1 bigf2)) 144 (return (list q (pquotient bigf1 q) (pquotient bigf2 q))))) 145 (setq xv (car bigf1)) 146 (setq bigf1 (pmodcontent bigf1)) 147 (setq bigf2 (pmodcontent bigf2)) 148 (setq c (pgcdu (setq c1 (car bigf1)) (setq c2 (car bigf2)))) 149 (setq bigf1 (cadr bigf1)) 150 (setq bigf2 (cadr bigf2)) 151 (setq n 0) 152 (setq e (pdegreer bigf2)) 153 (setq degree (pdegreer bigf1)) 154 (cond ((vgreat e degree) (setq e degree))) 155 (setq b (ash modulus -1)) 156 (setq gbar 157 (pgcdu (setq f1 (pgathercoef bigf1 xv 0)) 158 (setq f1f2 159 (pgathercoef bigf2 xv 0)))) 160 (cond ((equal 0 f1f2) (go step15a))) 161 (setq nubar (pdegree gbar xv)) 162 (setq nu1bar (+ nubar (pdegree bigf1 xv))) 163 (setq nu2bar (+ nubar (pdegree bigf2 xv))) 164 (setq f1f2 (ptimes f1 f1f2)) 165 (setq nubar (max nu1bar nu2bar)) 166 step6(setq b (cplus b 1)) 167 (cond ((equal (pcsubst f1f2 b xv) 0) (go step6))) 168 ;; Step 7 169 (setq gtilde (pcsubst gbar b xv)) 170 (setq f1tilde (pcsubst bigf1 b xv)) 171 (setq f2tilde (pcsubst bigf2 b xv)) 172 (setq biggtilde 173 (ptimeschk gtilde 174 (car (setq h2tilde (newgcd f1tilde f2tilde modulus))))) 175 (cond ((pcoefp biggtilde) (go step15a))) 176 (setq h1tilde (cadr h2tilde)) 177 (setq h2tilde (caddr h2tilde)) 178 (setq degree (pdegreer biggtilde)) 179 (cond ((vgreat degree e) (go step6)) 180 ((vgreat e degree) (setq n 0) (setq e degree))) 181 (setq n (1+ n)) 182 (cond ((equal n 1) (setq q (list xv 1 1 0 (cminus b))) 183 (setq gstar biggtilde) 184 (setq h1star h1tilde) 185 (setq h2star h2tilde)) 186 (t (setq gstar (lagrange33 gstar biggtilde q b)) 187 (setq h1star (lagrange33 h1star h1tilde q b)) 188 (setq h2star (lagrange33 h2star h2tilde q b)) 189 (setq q (ptimes q (list xv 1 1 0 (cminus b)))))) 190 ;; Step 12 191 (cond ((not (> n nubar)) (go step6))) 192 ;; Step 13 193 (cond ((or (not (= nu1bar (+ (setq degree (pdegree gstar xv)) 194 (pdegree h1star xv)))) 195 (not (= nu2bar (+ degree (pdegree h2star xv))))) 196 (setq n 0) 197 (go step6))) 198 (setq gstar (cadr (pmodcontent gstar))) 199 ;; Step 15 200 (setq q (pgathercoef gstar xv 0)) 201 (return (monicgcd (ptimeschk c gstar) 202 (ptimeschk (pquotient c1 c) (pquotientchk h1star q)) 203 (ptimeschk (pquotient c2 c) (pquotientchk h2star q)) 204 (leadcoefficient gstar))) 205 step15a 206 (return (list c 207 (ptimeschk (pquotient c1 c) bigf1) 208 (ptimeschk (pquotient c2 c) bigf2))) )) 209 210 211(defun monicgcd (gcd x y lcf) 212 (cond ((equal lcf 1) (list gcd x y)) 213 (t (list (ptimes (crecip lcf) gcd) 214 (ptimes lcf x) 215 (ptimes lcf y) )) )) 216 217;;*** PGCDM CORRESPONDS TO BROWN'S ALGORITHM M 218 219 220(defun pgcdm 221 (bigf1 bigf2) 222 (prog (c c1 c2 f1 f2 n 223 e degree mubar p 224 gtilde h1tilde h2tilde 225 modulus 226 biggtilde q h1star h2star 227 gstar xv gbar) 228 (setq p *alpha) 229 (setq xv (car bigf1)) 230 ;; Step 1 231 (setq f1 (pcontent bigf1)) 232 (setq f2 (pcontent bigf2)) 233 (setq c (cgcd (setq c1 (car f1)) (setq c2 (car f2)))) 234 (setq bigf1 (cadr f1)) 235 (setq bigf2 (cadr f2)) 236 ;; Step 3 237 (setq f1 (leadcoefficient bigf1)) 238 (setq f2 (leadcoefficient bigf2)) 239 (setq gbar (cgcd f1 f2)) 240 ;; Step 4 241 (setq n 0) 242 (setq degree (pdegreer bigf1)) 243 (setq e (pdegreer bigf2)) 244 (cond ((vgreat e degree) (setq e degree))) 245 ;; Step 5 246 (setq mubar 247 (* 2 gbar (max (maxcoefficient bigf1) 248 (maxcoefficient bigf2)))) 249 (go step6a) 250 step6(setq p (newprime p)) 251 step6a 252 (cond ((or (zerop (rem f1 p)) (zerop (rem f2 p))) 253 (go step6))) 254 (set-modulus p) 255 ;; Step 7 256 (setq gtilde (pmod gbar)) 257 ;; Step 8 258 (setq biggtilde 259 (ptimeschk gtilde 260 (car (setq h2tilde 261 (newgcd (pmod bigf1) (pmod bigf2) 262 modulus))))) 263 (cond ((pcoefp biggtilde) (setq modulus nil) 264 (setq gstar 1) 265 (setq h1star bigf1) 266 (setq h2star bigf2) 267 (go step15))) 268 (cond ((null (cdr h2tilde)) 269 (setq h1tilde (pquotient (pmod bigf1) (car h2tilde))) 270 (setq h2tilde (pquotient (pmod bigf2) (car h2tilde)))) 271 (t (setq h1tilde (cadr h2tilde)) 272 (setq h2tilde (caddr h2tilde)))) 273 (setq degree (pdegreer biggtilde)) 274 (cond ((vgreat degree e) (go step6)) 275 ((vgreat e degree) (setq n 0) (setq e degree))) 276 (setq n (1+ n)) 277 ;; Step 11 278 (set-modulus nil) 279 (cond ((equal n 1) (setq q p) 280 (setq gstar biggtilde) 281 (setq h1star h1tilde) 282 (setq h2star h2tilde)) 283 (t (setq gstar (lagrange3 gstar biggtilde p q)) 284 (setq h1star (lagrange3 h1star h1tilde p q)) 285 (setq h2star (lagrange3 h2star h2tilde p q)) 286 (setq q (* p q)))) 287 ;; Step 12 288 (cond ((> mubar q) (go step6))) 289 (cond ((> (* 2 (max (* (setq gtilde (norm gstar)) (maxcoefficient h1star)) 290 (* gtilde (maxcoefficient h2star)))) 291 q) 292 (go step6))) 293 (set-modulus nil) 294 (setq gstar (cadr (pcontent gstar))) 295 step15 296 (setq q (leadcoefficient gstar)) 297 (return (list (ptimeschk c gstar) 298 (ptimeschk (cquotient c1 c) (pquotientchk h1star q)) 299 (ptimeschk (cquotient c2 c) (pquotientchk h2star q)))))) 300 301;; THE FUNCTIONS ON THIS PAGE ARE USED BY KRONECKER FACTORING 302 303(defun pkroneck (p) 304 (prog (maxexp i l *p factors factor) 305 (setq maxexp (quotient (cadr p) 2)) 306 (setq i 1) 307 a (when (> i maxexp) (return (cons p factors))) 308 (setq l (p1 (reverse (let ((p p) (i i) ($factorflag t)) 309 (pfactor2 p i))))) 310 b (when (null l) (go d)) 311 (setq *l (car l)) 312 (setq *p (car p)) 313 (ignore-rat-err 314 (setq factor (errset (pinterpolate *l *p)))) 315 (setq l (cdr l)) 316 (if (atom factor) 317 (go b) 318 (setq factor (car factor))) 319 (when (or (pcoefp factor) 320 (not (equal (car p) (car factor))) 321 (not (pzerop (prem p factor)))) 322 (go b)) 323 (cond (modulus (pmonicize (cdr factor))) 324 ((pminusp factor) (setq factor (pminus factor)))) 325 (setq p (pquotient p factor)) 326 (setq maxexp (quotient (cadr p) 2)) 327 (setq factors (cons factor factors)) 328 (go a) 329 d (incf i) 330 (go a))) 331 332(defun pfactor2 (p i) 333 (cond ((< i 0) nil) 334 (t (cons (pfactor (pcsubst p i (car p))) 335 (pfactor2 p (1- i)))))) 336 337(defun rpowerset (x n) 338 (cond ((null x) (quote (1 nil))) 339 ((equal x 1) (quote (1))) 340 (t (cons 1 (ptts1 x n x))))) 341 342 343(defun allprods (x y) 344 (cond ((null x) nil) 345 ((null y) nil) 346 (t (nconc (ap1 (car x) y) (allprods (cdr x) y))))) 347 348(defun al1 (f r len) 349 (prog (ss) 350 (cond 351 ((equal len 1) 352 (return (mapcar #'(lambda (*y*) (cons *y* nil)) f))) 353 ((null r) (return nil)) 354 (t 355 (mapc #'(lambda (*y*) 356 (setq ss 357 (nconc ss 358 (mapcar #'(lambda (z) (cons z *y*)) 359 f)))) 360 (al1 (car r) (cdr r) (1- len))) 361 (return ss))))) 362 363 364(defun ap1 (x l) 365 (cond ((null l) nil) 366 (t (cons (ptimes x (car l)) (ap1 x (cdr l)))))) 367 368(defun ptts1 (x n y) 369 (cond ((equal n 1) (list y)) 370 (t (cons y (ptts1 x (1- n) (ptimes x y)))))) 371 372(defun p1 (l) 373 (prog (a) 374 (setq a (mapcar #'p11 l)) 375 (return (cond ((null l) nil) 376 (t (cdr (al1 (car a) 377 (cdr a) 378 (length a)))))))) 379 380(defun p11 (ele) 381 (cond ((null (cddr ele)) (rpowerset (car ele) (cadr ele))) 382 (t (allprods (rpowerset (car ele) (cadr ele)) 383 (p11 (cddr ele)))))) 384 385(defun pinterpolate (l var) 386 (psimp var (pinterpolate1 (pinterpolate2 l 1) 387 (- (length l) 2)))) 388 389(defun pinterpolate1 (x n) 390 (pinterpolate4 (pinterpolate5 (reverse x) 1 n n) (1+ n))) 391 392(defun pinterpolate2 (x n) 393 (cond ((null (cdr x)) x) 394 (t (cons (car x) 395 (pinterpolate2 (pinterpolate3 x n) (1+ n)))))) 396 397(defun pinterpolate3 (x n) 398 (cond ((null (cdr x)) nil) 399 (t (cons (pquotient (pdifference (cadr x) (car x)) n) 400 (pinterpolate3 (cdr x) n))))) 401 402(defun pinterpolate4 (x n) 403 (cond ((null x) nil) 404 ((pzerop (car x)) (pinterpolate4 (cdr x) (1- n))) 405 (t (cons n (cons (car x) 406 (pinterpolate4 (cdr x) (1- n))))))) 407 408(defun pinterpolate5 (x i j n) 409 (cond ((> i n) x) 410 (t (pinterpolate5 (cons (car x) (pinterpolate6 x i j)) 411 (1+ i) 412 (1- j) 413 n)))) 414 415(defun pinterpolate6 (x i j) 416 (cond ((zerop i) (cdr x)) 417 (t (cons (pdifference (cadr x) (pctimes j (car x))) 418 (pinterpolate6 (cdr x) (1- i) j))))) 419 420;; THE N**(1.585) MULTIPLICATION SCHEME 421;;FOLLOWS. IT SHOULD BE USED ONLY WHEN BOTH INPUTS ARE MULTIVARIATE, 422;;DENSE, AND OF NEARLY THE SAME SIZE. OR ABSOLUTELY TREMENDOUS. 423;;(THE CLASSICAL MULTIPLICATION SCHEME IS N**2 WHERE N IS SIZE OF 424;;POLYNOMIAL (OR N*M FOR DIFFERENT SIZES). FOR THIS 425;;CASE, N IS APPX. THE SIZE OF LARGER. 426 427(defmfun $fasttimes (x y) 428 (cond ((and (not (atom x)) (not (atom y)) 429 (equal (car x) (car y)) (equal (caar x) 'mrat) 430 (equal (cddr x) 1) (equal (cddr y) 1)) 431 (cons (car x)(cons (fptimes (cadr x)(cadr y))1))) 432 (t (merror (intl:gettext "fasttimes: arguments must be CRE polynomials with same variables."))))) 433 434(defun fptimes (x y) 435 (cond ((or (pzerop x) (pzerop y)) (pzero)) 436 ((pcoefp x) (pctimes x y)) 437 ((pcoefp y) (pctimes y x)) 438 ((eq (car x) (car y)) 439 (cond((or(univar(cdr x))(univar(cdr y))) 440 (cons (car x) (ptimes1 (cdr x) (cdr y)))) 441 (t(cons (car x) (fptimes1 (cdr x)(cdr y)))))) 442 ((pointergp (car x) (car y)) 443 (cons (car x) (pctimes1 y (cdr x)))) 444 (t (cons (car y) (pctimes1 x (cdr y)))))) 445 446(defun fptimes1 (f g) 447 (prog (a b c d) 448 (cond ((or (null f) (null g)) (return nil)) 449 ((null (cddr f)) 450 (return (lsft (pctimes1 (cadr f) g) (car f)))) 451 ((null (cddr g)) 452 (return (lsft (pctimes1 (cadr g) f) (car g))))) 453 (setq d (ash (1+ (max (car f) (car g))) -1)) 454 (setq f (halfsplit f d) g (halfsplit g d)) 455 (setq a (fptimes1 (car f) (car g))) 456 (setq b 457 (fptimes1 (ptptplus (car f) (cdr f)) (ptptplus (car g) (cdr g)))) 458 (setq c (fptimes1 (cdr f) (cdr g))) 459 (setq b (ptptdiffer (ptptdiffer b a) c)) 460 (return (ptptplus (lsft a (ash d 1)) (ptptplus (lsft b d) c))))) 461 462(defun halfsplit (p d) 463 (do ((a) (p p (cddr p))) 464 ((or (null p) (< (car p) d)) (cons (nreverse a) p)) 465 (setq a (cons (cadr p) (cons (- (car p) d) a))))) 466 467(defun lsft (p n) 468 (do ((q p (cddr (rplaca q (+ (car q) n))))) 469 ((null q))) 470 p) 471 472(declare-top (special wtsofar xweight $ratwtlvl v *x*)) 473 474;;; TO TRUNCATE ON E, DO RATWEIGHT(E,1); 475;;;THEN DO RATWTLVL:N. ALL POWERS >N GO TO 0. 476 477(defmfun $ratweight (&rest args) 478 (when (oddp (length args)) 479 (merror (intl:gettext "ratweight: number of arguments must be a multiple of 2."))) 480 (do ((l args (cddr l))) 481 ((null l)) 482 (rplacd (or (assoc (first l) *ratweights :test #'equal) 483 (car (push (list (first l)) *ratweights))) 484 (second l))) 485 (setq $ratweights (cons '(mlist simp) (dot2l *ratweights))) 486 (if (null args) 487 $ratweights 488 (cons '(mlist) args))) 489 490(defun pweight (x) 491 (or (get x '$ratweight) 0)) 492 493(defun wtptimes (x y wtsofar) 494 (cond ((or (pzerop x) (pzerop y) (> wtsofar $ratwtlvl)) 495 (pzero)) 496 ((pcoefp x) (wtpctimes x y)) 497 ((pcoefp y) (wtpctimes y x)) 498 ((eq (car x) (car y)) 499 (palgsimp (car x) 500 (wtptimes1 (cdr x) 501 (cdr y) 502 (pweight (car x))) 503 (alg x))) 504 ((pointergp (car x) (car y)) 505 (psimp (car x) 506 (wtpctimes1 y (cdr x) (pweight (car x))))) 507 (t (psimp (car y) 508 (wtpctimes1 x (cdr y) (pweight (car y))))))) 509 510(defun wtptimes1 (*x* y xweight) 511 (prog (u* v) 512 (declare (special v)) 513 (setq v (setq u* (wtptimes2 y))) 514 a (setq *x* (cddr *x*)) 515 (cond ((null *x*) (return u*))) 516 (wtptimes3 y) 517 (go a))) 518 519 520(defun wtptimes2 (y) 521 (if (null y) 522 nil 523 (let ((ii (+ (* xweight (+ (car *x*) (car y))) wtsofar))) 524 (if (> ii $ratwtlvl) 525 (wtptimes2 (cddr y)) 526 (pcoefadd (+ (car *x*) (car y)) 527 (wtptimes (cadr *x*) (cadr y) ii) 528 (wtptimes2 (cddr y))))))) 529 530(defun wtptimes3 (y) 531 (prog ((e 0) u c) 532 (declare (special v)) 533 a1 (cond ((null y) (return nil))) 534 (setq e (+ (car *x*) (car y))) 535 (setq c (wtptimes (cadr y) (cadr *x*) (+ wtsofar (* xweight e)))) 536 (cond ((pzerop c) (setq y (cddr y)) (go a1)) 537 ((or (null v) (> e (car v))) (setq u* (setq v (ptptplus u* (list e c)))) (setq y (cddr y)) (go a1)) 538 ((equal e (car v)) 539 (setq c (pplus c (cadr v))) 540 (cond ((pzerop c) (setq u* (setq v (ptptdiffer u* (list (car v) (cadr v)))))) (t (rplaca (cdr v) c))) 541 (setq y (cddr y)) 542 (go a1))) 543 a (cond ((and (cddr v) (> (caddr v) e)) (setq v (cddr v)) (go a))) 544 (setq u (cdr v)) 545 b (cond ((or (null (cdr u)) (< (cadr u) e)) (rplacd u (cons e (cons c (cdr u)))) (go e))) 546 (cond ((pzerop (setq c (pplus (caddr u) c))) (rplacd u (cdddr u)) (go d)) (t (rplaca (cddr u) c))) 547 e (setq u (cddr u)) 548 d (setq y (cddr y)) 549 (cond ((null y) (return nil)) 550 ((pzerop 551 (setq c (wtptimes (cadr *x*) (cadr y) 552 (+ wtsofar (* xweight 553 (setq e (+ (car *x*) (car y)))))))) 554 (go d))) 555 c (cond ((and (cdr u) (> (cadr u) e)) (setq u (cddr u)) (go c))) 556 (go b))) 557 558 559(defun wtpctimes (c p) 560 (cond ((pcoefp p) (ctimes c p)) 561 (t (psimp (car p) (wtpctimes1 c (cdr p) (pweight (car p))))))) 562 563(defun wtpctimes1 (c x xwt) 564 (prog (cc) 565 (return 566 (cond ((null x) nil) 567 (t (setq cc (wtptimes c 568 (cadr x) 569 (+ wtsofar (* xwt (car x))))) 570 (cond ((pzerop cc) (wtpctimes1 c (cddr x) xwt)) 571 (t (cons (car x) 572 (cons cc 573 (wtpctimes1 c 574 (cddr x) 575 xwt)))))))))) 576 577(defun wtpexpt (x n) 578 (cond ((= n 0) 1) 579 ((= n 1) x) 580 ((evenp n) 581 (let ((xn2 (wtpexpt x (/ n 2)))) 582 (wtptimes xn2 xn2 0))) 583 (t (wtptimes x (wtpexpt x (1- n)) 0)))) 584 585(defmfun $horner (e &rest l) 586 (let (($ratfac nil) 587 (varlist (cdr $ratvars)) 588 genvar 589 (x nil) 590 (arg1 (taychk2rat e))) 591 (cond ((mbagp arg1) 592 (cons (car arg1) 593 (mapcar #'(lambda (u) (apply '$horner (cons u l))) (cdr arg1)))) 594 (t 595 (setq x (apply #'$rat (cons arg1 l))) 596 (mapc #'(lambda (y z) (putprop y z 'disrep)) (cadddr (car x)) (caddar x)) 597 (div* (hornrep (cadr x)) (hornrep (cddr x))))))) 598 599(defun hornrep (p) 600 (if (pcoefp p) 601 p 602 (horn+ (cdr p) (get (car p) 'disrep)))) 603 604(defun horn+ (l var) 605 (prog (ans last) 606 (setq ans (hornrep (cadr l))) 607 a (setq last (car l) l (cddr l)) 608 (cond ((null l) 609 (return (cond ((equal last 0) ans) 610 (t (list '(mtimes) 611 (list '(mexpt) var last) ans))))) 612 (t (setq ans (list '(mplus) 613 (hornrep (cadr l)) 614 (list '(mtimes) 615 (list '(mexpt) var (- last (car l))) 616 ans))))) 617 (go a))) 618 619(declare-top (special y genvar $savefactors checkfactors 620 exp var x $factorflag $ratfac 621 ratform 622 wholepart parnumer varlist n)) 623 624(defmfun $partfrac (exp var) 625 (cond ((mbagp exp) 626 (cons (car exp) (mapcar #'(lambda (u) ($partfrac u var)) (cdr exp)))) 627 ((and (atom var) (not (among var exp))) exp) 628 (t (let (($savefactors t) (checkfactors ()) (varlist (list var)) 629 $ratfac $algebraic $keepfloat ratform genvar) 630 (desetq (ratform . exp) (taychk2rat exp)) 631 (setq var (caadr (ratf var))) 632 (setq exp (partfrac exp var)) 633 (setq exp (cons (car exp) ;FULL DECOMP? 634 (mapcan #'partfraca (cdr exp)))) 635 (add2* (disrep (car exp)) 636 (cons '(mplus) 637 (mapcar #'(lambda (l) 638 (destructuring-let (((coef poly exp) l)) 639 (list '(mtimes) 640 (disrep coef) 641 (list '(mexpt) 642 (disrep poly) 643 (- exp))))) 644 (cdr exp)))))))) 645 646(defun partfraca (llist) 647 (destructuring-let (((coef poly exp) llist)) 648 (do ((nc (ratdivide coef poly) (ratdivide (car nc) poly)) 649 (n exp (1- n)) 650 (ans)) 651 ((rzerop (car nc)) (cons (list (cdr nc) poly n) ans)) 652 (push (list (cdr nc) poly n) ans)))) 653 654(defun partfrac (rat var) 655 (destructuring-let* (((wholepart frpart) (pdivide (car rat) (cdr rat))) 656 ((num . denom) (ratqu frpart (cdr rat)))) 657 (cond 658 ((pzerop num) (cons wholepart nil)) 659 ((or (pcoefp denom) (pointergp var (car denom))) (cons rat nil)) 660 (t (destructuring-let (((content bpart) (oldcontent denom))) 661 (let (apart y parnumer) 662 (loop 663 for (factor multiplicity) 664 on (pfactor bpart) by #'cddr 665 unless (zerop (pdegree factor var)) 666 do 667 (setq apart (pexpt factor multiplicity) 668 bpart (pquotient bpart apart) 669 y (bprog apart bpart) 670 frpart (cdr (ratdivide (ratti num (cdr y) t) 671 apart))) 672 (push (list (ratqu frpart content) factor multiplicity) 673 parnumer) 674 (desetq (num . content) 675 (cdr (ratdivide (ratqu (ratti num (car y) t) 676 content) 677 bpart)))) 678 (cons wholepart parnumer))))))) 679 680(declare-top (unspecial exp f n ss v var xv y *chk *l *max *p 681 *res u* *x* *y*)) 682 683;; $RATDIFF TAKES DERIVATIVES FAST. IT ASSUMES THAT THE 684;; ONLY ENTITY WHICH DEPENDS ON X IS X ITSELF. 685;; THAT IS, DEPENDENCIES DECLARED EXPLICITLY OR IMPLICITLY ARE 686;; TOTALLY IGNORED. RATDIFF(F(X),X) IS 0. RATDIFF(Y,X) IS 0. 687;; ANY OTHER USAGE MUST GO THROUGH $DIFF. 688;; FURTHERMORE, X IS ASSUMED TO BE AN ATOM OR A SINGLE ITEM ON 689;; VARLIST. E.G. X MIGHT BE SIN(U), BUT NOT 2*SIN(U). 690 691(declare-top (special varlist genvar x)) 692 693(defmfun $ratdiff (p x) 694 (if ($ratp p) 695 (setq p (minimize-varlist 696 (if (member 'trunc (cdar p) :test #'eq) ($taytorat p) p)))) 697 (let ((formflag ($ratp p)) (varlist) (genvar)) 698 (newvar x) (newvar p) 699 (or (every #'(lambda (exp) 700 (or (alike1 x exp) (free exp x))) 701 varlist) 702 (merror (intl:gettext "ratdiff: first argument must be a polynomial in ~M; found: ~M") x p)) 703 (setq p (ratf p)) 704 (setq x (caadr (ratf x))) 705 (setq p (cons (car p) (ratderivative (cdr p) x))) 706 (if formflag p ($ratdisrep p)))) 707 708(declare-top (unspecial x)) 709 710(declare-top (special $pfeformat varlist $factorflag m v dosimp)) 711 712(defmfun $pfet (m) 713 (prog (listov $pfeformat varlist $factorflag) 714 (setq $pfeformat t) 715 (newvar m) 716 (setq listov varlist) 717 (mapc #'(lambda (r) (setq m (pfet1 m r))) 718 listov) 719 (setq m (simplify m)) 720 (setq m (cond ((atom m) m) 721 ((eq (caar m) 'mplus) 722 (cons '(mplus) 723 (mapcar #'$ratexpand (cdr m)))) 724 (t ($ratexpand m)))) 725 (return (cond ((atom m) m) 726 ((eq (caar m) 'mplus) 727 (cons '(mplus) 728 (mapcar #'sssqfr (cdr m)))) 729 (t (sssqfr m)))))) 730 731(defun sssqfr (x) 732 (let ((dosimp t)) (simplify ($sqfr x)))) 733 734(defun pfet1 (m v) 735 (cond ((atom m) m) 736 ((eq (caar m) 'mplus) 737 (cons '(mplus) 738 (mapcar #'(lambda (s) ($partfrac s v)) 739 (cdr m)))) 740 (t ($partfrac m v)))) 741 742(declare-top (unspecial m v)) 743