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 nrat4) 14 15(declare-top (special $ratsimpexpons *exp *exp2 *radsubst *loglist $radsubstflag 16 $logsimp *v *var radcanp)) 17 18(defmvar $radsubstflag nil 19 "`radsubstflag' `t' makes `ratsubs' call `radcan' when it appears useful") 20 21 22(defun pdis (x) ($ratdisrep (pdis* x))) 23 24(defun pdis* (x) `((mrat simp ,varlist ,genvar) ,x . 1)) 25 26(defun rdis (x) ($ratdisrep (rdis* x))) 27 28(defun rdis* (x) `((mrat simp ,varlist ,genvar) . ,x)) 29 30(defun rform (x) (cdr (ratf x))) 31 32(setq radcanp nil) 33 34(defmfun $ratcoef (e x &optional (n 1)) 35 (ratcoeff e x n)) ; The spelling "ratcoeff" is nicer. 36 37(defun ratcoeff (a b c) 38 (let* ((formflag ($ratp a)) 39 (taylorform (and formflag (member 'trunc (cdar a) :test #'eq)))) 40 (cond ((zerop1 b) (improper-arg-err b '$ratcoeff)) 41 ((mbagp a) (cons (car a) 42 (mapcar #'(lambda (a) (ratcoeff a b c)) 43 (cdr a)))) 44 ((and taylorform (mnump c) (assolike b (cadddr (cdar a)))) 45 (pscoeff1 a b c)) 46 ((and taylorform (mexptp b) (mnump c) (mnump (caddr b)) 47 (assolike (cadr b) (cadddr (cdar a)))) 48 (pscoeff1 a (cadr b) (mul2 c (caddr b)))) 49 ((and taylorform (equal c 0)) a) 50 (t (if taylorform (setq a (ratdisrep a))) 51 (setq a (let ($ratwtlvl) 52 (if (equal c 0) 53 (ratcoef (mul2* a b) b) 54 (ratcoef a (if (equal c 1) b (list '(mexpt) b c)))))) 55 (if (and formflag (not taylorform)) 56 (minimize-varlist a) 57 (ratdisrep a)))))) 58 59(defun minimize-varlist (ratfun) 60 (if (not ($ratp ratfun)) (setq ratfun (ratf ratfun))) 61 (minvarlist-mrat (caddr (car ratfun)) (cadddr (car ratfun)) 62 (cdr ratfun))) 63 64(defun minvarlist-mrat (vars gens ratform) 65 (let ((newgens (union* (listovars (car ratform)) 66 (listovars (cdr ratform))))) 67 (do ((lv vars (cdr lv)) 68 (lg gens (cdr lg)) 69 (nlv ()) 70 (nlg ())) 71 ((null lg) 72 (cons (list 'mrat 'simp (nreverse nlv) (nreverse nlg)) 73 ratform)) 74 (cond ((member (car lg) newgens :test #'eq) 75 (push (car lg) nlg) 76 (push (car lv) nlv)))))) 77 78(defun ratcoef (exp var) 79 (prog (varlist genvar $ratfac $algebraic $ratwtlvl bas minvar) 80 (setq var (ratdisrep var)) 81 (setq bas (if (and (mexptp var) (mnump (caddr var))) (cadr var) var)) 82 (newvar var) 83 (newvar bas) 84 (setq minvar (car varlist)) 85 (newvar exp) 86 (setq exp (cdr (ratrep* exp))) 87 (setq var (cdr (ratrep* var))) 88 (setq bas (cadr (ratrep* bas))) 89 (if (and (onep1 (cdr exp)) (onep1 (cdr var)) (pureprod (car var))) 90 (return (pdis* (prodcoef (car var) (car exp))))) 91 (setq exp (ratquotient exp var)) 92 (if (null minvar) (return (pdis* (prodcoef (cdr exp) (car exp))))) 93 (setq minvar (caadr (ratrep* minvar))) 94 loop (if (or (pcoefp (cdr exp)) (pointergp minvar (cadr exp))) 95 (return (rdis* (cdr (ratdivide exp bas))))) 96 (setq exp (ratcoef1 (car exp) (cdr exp))) 97 (go loop))) 98 99(defun ratcoef1 (num den) 100 (cond ((pcoefp num) (rzero)) 101 ((eq (car num) (car den)) (car (pdivide num den))) 102 ((pointergp (car den) (car num)) (rzero)) 103 (t (ratcoef1 (constcoef (cdr num)) den)))) 104 105(defun constcoef (p) 106 (cond ((null p) 0) 107 ((zerop (car p)) (cadr p)) 108 (t (constcoef (cddr p))))) 109 110(setq *radsubst nil) 111 112(defmfun $ratsubst (a b c) ; NEEDS CODE FOR FAC. FORM 113 (prog (varlist newvarlist dontdisrepit $ratfac genvar $keepfloat) 114 ;; hard to maintain user ordering info. 115 (if ($ratp c) (setq dontdisrepit t)) 116 (if (and $radsubstflag 117 (prog2 (newvar b) (some #'mexptp varlist))) 118 (let (($factorflag t) *exp *exp2 *radsubst) 119 (setq b (fullratsimp b)) 120 (setq c (fullratsimp c)) 121 (setq varlist nil) 122 (fnewvar b) 123 (fnewvar c) 124 (setq *exp (cdr (ratrep* b))) 125 (setq *exp2 (cdr (ratrep* c))) 126 ;; since *radsubst is t, both *exp and *exp2 will be radcan simplified 127 (setq *radsubst t) 128 (spc0) 129 (setq b (rdis *exp) c (rdis *exp2)) 130 (setq varlist nil)) 131 (setq varlist nil)) 132 (setq a ($ratdisrep a) b ($ratdisrep b) c ($ratdisrep c)) 133 (cond ((integerp b) (setq c (ratf (maxima-substitute a b c))) 134 (return (cond (dontdisrepit c) (t ($ratdisrep c)))))) 135 (newvar c) 136 (setq 137 newvarlist 138 (mapcar 139 #'(lambda (z) 140 (cond ((atom z) z) 141 (t (resimplify 142 (cons (car z) 143 (mapcar #'(lambda (zz) 144 (cond ((alike1 zz b) a) 145 ((atom zz) zz) 146 (t ($ratdisrep 147 ($ratsubst a b zz))))) 148 (cdr z))))))) 149 varlist)) 150 (newvar a) (newvar b) 151 (setq newvarlist (reverse (pairoff (reverse varlist) 152 (reverse newvarlist)))) 153 (setq a (cdr (ratrep* a))) 154 (setq b (cdr (ratrep* b))) 155 (setq c (cdr (ratrep* c))) 156 (when (pminusp (car b)) 157 (setq b (ratminus b)) 158 (setq a (ratminus a))) 159 (when (and (equal 1 (car b)) 160 (not (equal 1 (cdr b))) 161 (not (equal 0 (car a)))) 162 (setq a (ratinvert a)) 163 (setq b (ratinvert b))) 164 (cond ((not (equal 1 (cdr b))) 165 (setq a (rattimes a (cons (cdr b) 1) t)) 166 (setq b (cons (car b) 1)))) 167 (setq c 168 (cond ((member (car b) '(0 1) :test #'equal) 169 (ratf (maxima-substitute (rdis a) b (rdis c)))) 170 (t (cons (list 'mrat 'simp varlist genvar) 171 (if (equal (cdr a) 1) 172 (ratreduce (everysubst0 (car a) (car b) (car c)) 173 (everysubst0 (car a) (car b) (cdr c))) 174 (allsubst00 a b c)))))) 175 (unless (alike newvarlist varlist) 176 (setq varlist newvarlist 177 c (rdis (cdr c)) 178 varlist nil 179 c (ratf c))) 180 (return (cond (dontdisrepit c) (t ($ratdisrep c)))))) 181 182(defun xptimes (x y) (if $ratwtlvl (wtptimes x y 0) (ptimes x y))) 183 184(defun allsubst00 (a b c) 185 (cond ((equal a b) c) 186 ((not (equal (cdr b) 1)) c) 187 (t (ratquotient (everysubst00 a (car b) (car c)) 188 (everysubst00 a (car b) (cdr c)))))) 189 190(defun everysubst00 (x i z) 191 (loop with ans = (rzero) 192 for (exp coef) on (everysubst i z *alpha) by #'cddr 193 do (setq ans (ratplus ans (rattimes (cons coef 1) (ratexpt x exp) t))) 194 finally (return ans))) 195 196(defun everysubst0 (x i z) 197 (loop with ans = (pzero) 198 for (exp coef) on (everysubst i z *alpha) by #'cddr 199 do (setq ans (pplus ans (xptimes coef (pexpt x exp)))) 200 finally (return ans))) 201 202(defun everysubst1 (a b maxpow) 203 (loop for (exp coef) on (p-terms b) by #'cddr 204 for part = (everysubst a coef maxpow) 205 nconc (if (= 0 exp) part 206 (everysubst2 part (make-poly (p-var b) exp 1))))) 207 208(defun everysubst2 (l h) 209 (do ((ptr l (cddr ptr))) 210 ((null ptr) l) 211 (setf (cadr ptr) (ptimes h (cadr ptr))))) 212 213 214(defun pairoff (l m) 215 (cond ((null m) l) (t (cons (car m) (pairoff (cdr l) (cdr m)))))) 216 217;;(DEFUN PAIROFF (L M) 218;; ;(COND ((NULL M) L) (T (CONS (CAR M) (PAIROFF (CDR L) (CDR M))))) 219;; (let ((ans nil)) 220;; (dolist (x m (nreconc ans l)) 221;; (push x ans) (setq l (cdr l))))) 222 223(defun everysubst (a b maxpow) 224 (cond ((pcoefp a) 225 (cond ((equal a 1) (list maxpow b)) 226 ((pcoefp b) 227 (list (setq maxpow 228 (do ((b b (quotient b a)) 229 (ans 0 (1+ ans))) 230 ((or (> (abs a) (abs b)) 231 (equal maxpow ans)) 232 ans))) 233 (quotient b (setq maxpow (expt a maxpow))) 234 0 235 (rem b maxpow))) 236 (t (everysubst1 a b maxpow)))) 237 ((or (pcoefp b) (pointergp (car a) (car b))) (list 0 b)) 238 ((eq (car a) (car b)) 239 (cond ((null (cdddr a)) (everypterms b (caddr a) (cadr a) maxpow)) 240 (t (substforsum a b maxpow)))) 241 (t (everysubst1 a b maxpow)))) 242 243(defun everypterms (x p n maxpow) 244 (if (< (cadr x) n) 245 (list 0 x) 246 (prog (k ans q part) 247 (setq k (car x)) 248 (setq x (cdr x)) 249 l (setq q (min maxpow (quotient (car x) n))) 250 m (when (equal q 0) 251 (return (if (null x) 252 ans 253 (cons 0 (cons (psimp k x) ans))))) 254 (setq part (everysubst p (cadr x) q)) 255 (setq ans (nconc (everypterms1 part k n (car x)) ans)) 256 (setq x (cddr x)) 257 (when (null x) 258 (setq q 0) 259 (go m)) 260 (go l)))) 261 262(defun everypterms1 (l k n j) 263 (do ((ptr l (cddr ptr))) 264 ((null ptr) l) 265 (setf (cadr ptr) 266 (ptimes (psimp k (list (- j (* n (car ptr))) 1)) 267 (cadr ptr))))) 268 269(defun substforsum (a b maxpow) 270 (do ((pow 0 (1+ pow)) 271 (quot) (zl-rem) (ans)) 272 ((not (< pow maxpow)) (list* maxpow b ans)) 273 (desetq (quot zl-rem) (pdivide b a)) 274 (unless (and (equal (cdr quot) 1) 275 (not (pzerop (car quot))) 276 (equal (cdr zl-rem) 1)) 277 (return (cons pow (cons b ans)))) 278 (unless (pzerop (car zl-rem)) 279 (setq ans (cons pow (cons (car zl-rem) ans)))) 280 (setq b (car quot)))) 281 282(defun prodcoef (a b) 283 (cond ((pcoefp a) 284 (cond ((pcoefp b) (quotient b a)) (t (prodcoef1 a b)))) 285 ((pcoefp b) (pzero)) 286 ((pointergp (car a) (car b)) (pzero)) 287 ((eq (car a) (car b)) 288 (cond ((null (cdddr a)) 289 (prodcoef (caddr a) (ptterm (cdr b) (cadr a)))) 290 (t (sumcoef a b)))) 291 (t (prodcoef1 a b)))) 292 293(defun sumcoef (a b) 294 (desetq (a b) (pdivide b a)) 295 (if (and (equal (cdr a) 1) (equal (cdr b) 1)) 296 (car a) 297 (pzero))) 298 299(defun prodcoef1 (a b) 300 (loop with ans = (pzero) 301 for (bexp bcoef) on (p-terms b) by #'cddr 302 for part = (prodcoef a bcoef) 303 unless (pzerop part) 304 do (setq ans (pplus ans (psimp (p-var b) (list bexp part)))) 305 finally (return ans))) 306 307(defun pureprod (x) 308 (or (atom x) 309 (and (not (atom (cdr x))) 310 (null (cdddr x)) 311 (pureprod (caddr x))))) 312 313(defmfun $bothcoef (r var) 314 (prog (*var h varlist genvar $ratfac) 315 (unless ($ratp r) 316 (return `((mlist) 317 ,(setq h (coeff r var 1.)) 318 ((mplus) ,r ((mtimes) -1 ,h ,var))))) 319 (newvar var) 320 (setq h (and varlist (car varlist))) 321 (newvar r) 322 (setq var (cdr (ratrep* var))) 323 (setq r (cdr (ratrep* r))) 324 (and h (setq h (caadr (ratrep* h)))) 325 (cond ((and h (or (pcoefp (cdr r)) (pointergp h (cadr r))) 326 (equal 1 (cdr var))) 327 (setq var (bothprodcoef (car var) (car r))) 328 (return (list '(mlist) 329 (rdis* (ratreduce (car var) (cdr r))) 330 (rdis* (ratreduce (cdr var) (cdr r)))))) 331 (t 332 ;; CAN'T TELL WHAT BROUGHT US TO THIS POINT, SORRY 333 (merror (intl:gettext "bothcoef: invalid arguments.")))))) 334 335;;COEFF OF A IN B 336 337(defun bothprodcoef (a b) 338 (let ((c (prodcoef a b))) 339 (if (pzerop c) (cons (pzero) b) (cons c (pdifference b (ptimes c a)))))) 340 341(defvar argsfreeofp nil) 342 343(defun argsfreeof (var e) 344 (let ((argsfreeofp t)) (freeof var e))) 345 346;;; This is a version of freeof for a list first argument 347(defmfun $lfreeof (l e) "`freeof' for a list first argument" 348 (unless ($listp l) 349 (merror (intl:gettext "lfreeof: first argument must be a list; found: ~M") l)) 350 (let ((exp ($totaldisrep e))) 351 (dolist (var (margs l) t) 352 (unless (freeof ($totaldisrep var) exp) (return nil))))) 353 354(defmfun $freeof (&rest args) 355 (prog (l e) 356 (setq l (mapcar #'$totaldisrep (nreverse args)) 357 e (car l)) 358 loop (or (setq l (cdr l)) (return t)) 359 (if (freeof (getopr (car l)) e) (go loop)) 360 (return nil))) 361 362(defun freeof (var e) 363 (cond ((alike1 var e) nil) 364 ((atom e) t) 365 ((and (not argsfreeofp) 366 (or (alike1 var ($verbify (caar e))) 367 (alike1 var ($nounify (caar e))))) 368 nil) 369 ((and (or (member (caar e) '(%product %sum %laplace) :test #'eq) 370 (and (eq (caar e) '%integrate) (cdddr e)) 371 (and (eq (caar e) '%limit) (cddr e))) 372 (alike1 var (caddr e))) 373 (freeofl var (cdddr e))) 374 ((eq (caar e) '%at) 375 (cond ((not (freeofl var (hand-side (caddr e) 'r))) nil) 376 ((not (freeofl var (hand-side (caddr e) 'l))) t) 377 (t (freeof var (cadr e))))) 378 ((and (eq (caar e) 'lambda) 379 (not (member 'array (cdar e) :test #'eq)) 380 ($listp (cadr e)) 381 ; Check if var appears in the lambda list in any of the 382 ; following ways: var, 'var, [var] or ['var]. 383 (some (lambda (v) 384 (or (eq v var) 385 (alike1 v `((mquote) ,var)) 386 (alike1 v `((mlist) ,var)) 387 (alike1 v `((mlist) ((mquote) ,var))))) 388 (cdadr e))) 389 t) 390 ;; Check for a local variable in a block. 391 ((and (eq (caar e) 'mprog) 392 ($listp (cadr e)) 393 ; Check if var appears in the variable list alone or 394 ; in an assignment 395 (some (lambda (v) 396 (or (eq v var) 397 (and (msetqp v) 398 (eq (cadr v) var)))) 399 (cdadr e))) 400 t) 401 ;; Check for a loop variable. 402 ((and (member (caar e) '(mdo mdoin) :test #'eq) 403 (alike1 var (cadr e))) 404 t) 405 (argsfreeofp (freeofl var (margs e))) 406 (t (freeofl var (cdr e))))) 407 408(defun freeofl (var l) (loop for x in l always (freeof var x))) 409 410(defun hand-side (e flag) 411 (setq e (if (eq (caar e) 'mequal) (ncons e) (cdr e))) 412 (mapcar #'(lambda (u) (if (eq flag 'l) (cadr u) (caddr u))) e)) 413 414;; subtitle radcan 415 416(defmfun $radcan (exp) 417 (cond ((mbagp exp) (cons (car exp) (mapcar '$radcan (cdr exp)))) 418 (t (let (($ratsimpexpons t)) 419 (simplify (let (($expop 0) ($expon 0)) 420 (radcan1 (fr1 exp nil)))))))) 421 422(defun radcan1 (*exp) 423 (cond ((atom *exp) *exp) 424 (t (let (($factorflag t) varlist genvar $ratfac $norepeat 425 ($gcd (or $gcd (car *gcdl*))) 426 (radcanp t)) 427 (newvar *exp) 428 (setq *exp (cdr (ratrep* *exp))) 429 (setq varlist 430 (mapcar 431 #'(lambda (x) (cond 432 ((atom x) x) 433 (t (cons (car x) 434 (mapcar 'radcan1 (cdr x)))))) 435 varlist)) 436 (spc0) 437 (fr1 (rdis *exp) nil))))) 438 439(defun spc0 () 440 (prog (*v *loglist) 441 (if (allatoms varlist) (return nil)) 442 (setq varlist (mapcar #'spc1 varlist)) ;make list of logs 443 (setq *loglist (factorlogs *loglist)) 444 (mapc #'spc2 *loglist) ;subst log factorizations 445 (mapc #'spc3 varlist genvar) ;expand exponents 446 (mapc #'spc4 varlist) ;make exponent list 447 (desetq (varlist . genvar) (spc5 *v varlist genvar)) 448 ;find expon dependencies 449 (setq varlist (mapcar #'rjfsimp varlist)) ;restore radicals 450 (mapc #'spc7 varlist))) ;simplify radicals 451 452(defun allatoms (l) 453 (loop for x in l always (atom x))) 454 455(defun rjfsimp (x &aux expon) 456 (cond ((and *radsubst $radsubstflag) x) 457 ((not (m$exp? (setq x (let ($logsimp) (resimplify x))))) x) 458 ((mlogp (setq expon (caddr x))) (cadr expon)) 459 ((not (and (mtimesp expon) (or $logsimp *var))) x) 460 (t (do ((rischflag (and *var (not $logsimp) (not (freeof *var x)))) 461 (power (cdr expon) (cdr power))) ;POWER IS A PRODUCT 462 ((null power) x) 463 (cond ((numberp (car power))) 464 ((mlogp (car power)) 465 (and rischflag (cdr power) (return x)) 466 (return 467 `((mexpt) ,(cadar power) 468 ,(muln (remove (car power) (cdr expon) :count 1 :test #'equal) 469 nil)))) 470 (rischflag (return x))))))) 471 472(defun dsubsta (x y zl) 473 (cond ((null zl) zl) 474 (t (cond ((alike1 y (car zl)) (rplaca zl x)) 475 ((not (atom (car zl))) (dsubsta x y (cdar zl)))) 476 (dsubsta x y (cdr zl)) 477 zl))) 478 479(defun radsubst (a b) 480 (setq *exp (allsubst00 a b *exp)) 481 (if *radsubst (setq *exp2 (allsubst00 a b *exp2)))) 482 483(setq *var nil) 484 485(defun spc1 (x) 486 (cond ((mlogp x) (putonloglist x)) 487 ((and (mexptp x) (not (eq (cadr x) '$%e))) 488 ($exp-form (list '(mtimes) 489 (caddr x) 490 (putonloglist (list '(%log simp ratsimp) 491 (cadr x)))))) 492 (t x))) 493 494(defun putonloglist (l) 495 (unless (memalike l *loglist) (push l *loglist)) 496 l) 497 498(defun spc2 (p) 499 (radsubst (rform (cdr p)) (rform (car p))) 500 (dsubsta (cdr p) (car p) varlist)) 501 502(defun spc2a (x) ;CONVERTS FACTORED 503 (let ((sum (mapcar #'spc2b x))) ;RFORM LOGAND TO SUM 504 (if (cdr sum) ;OF LOGS 505 (cons '(mplus) sum) 506 (car sum)))) 507 508(defun spc2b (x) 509 (let ((log `((%log simp ratsimp irreducible) ,(pdis (car x))))) 510 (if (equal 1 (cdr x)) log 511 (list '(mtimes) (cdr x) log)))) 512 513(defun spc3 (x v &aux y) 514 (when (and (m$exp? x) 515 (not (atom (setq y (caddr x)))) 516 (mplusp (setq y (expand1 (if *var ($partfrac y *var) y) 10 10)))) 517 (setq y (cons '(mtimes) 518 (mapcar #'(lambda (z) ($ratsimp ($exp-form z))) (cdr y)))) 519 (radsubst (rform y) (rget v)) 520 (dsubsta y x varlist))) 521 522(defun spc4 (x) 523 (if (and (m$exp? x) 524 (not (memalike (caddr x) *v))) 525 (push (caddr x) *v))) 526 527(defun rzcontent (r) 528 (destructuring-let (((c1 p) (pcontent (car r))) 529 ((c2 q) (pcontent (cdr r)))) 530 (if (pminusp p) (setq p (pminus p) c1 (cminus c1))) 531 (cons (cons c1 c2) (cons p q)))) 532 533;;The GCDLIST looks like (( GCM1pair occurrencepair11 occurrencepair12 ...) ... 534;;(GCMnpair occurrencepairn1 occurrencepairn2 ...)) 535;;where GCMpairs are lists of ratforms and prefix forms for the greatest common 536;;multiple of the occurrencepairs. Each of these pairs is a list of a ratform 537;;and a prefix form. The prefix form is a pointer into the varlist. 538;;The occurrences are exponents of the base %E. 539 540(defun spc5 (vl oldvarlist oldgenvar &aux gcdlist varlist genvar) 541 (dolist (v vl) 542 (destructuring-let* ((((c1 . c) . r) (rzcontent (rform v))) 543 (g (assoc r gcdlist :test #'equal))) 544 (cond (g (setf (cadr g) (plcm c (cadr g))) 545 (push (list ($exp-form (div* v c1)) c) (cddr g))) 546 (t (push (list r c (list ($exp-form (div* v c1)) c)) gcdlist))))) 547 (dolist (g gcdlist) 548 (let ((rd (rdis (car g)))) 549 (when (and (mlogp rd) (memalike (cadr rd) oldvarlist)) 550 (push (list (cadr rd) 1) (cddr g))) 551 (rplaca g ($exp-form (div rd (cadr g)))))) 552 (spc5b gcdlist oldvarlist oldgenvar)) 553 554;;(DEFUN SPC5B (V VARLIST GENVAR) 555;; (DOLIST (L V) 556;; (DOLIST (X (CDDR L)) 557;; (UNLESS (EQUAL (CADR L) (CADR X)) 558;; (RADSUBST (RATEXPT (RFORM (CAR L)) 559;; (CAR (QUOTIENT (CADR X) (CADR L)))) 560;; (RFORM (CAR X)))))) 561;; (CONS VARLIST GENVAR)) 562 563 564(defun spc5b (v varlist genvar) 565 (dolist (l v) 566 (dolist (x (cddr l)) 567 (unless (equal (cadr l) (cadr x)) 568 (radsubst (ratexpt (rform (car l)) 569 (quotient (cadr l) (cadr x))) 570 (rform (car x)))))) 571 (cons varlist genvar)) 572 573(defun spc7 (x) 574 (if (eq x '$%i) (setq x '((mexpt) -1 ((rat) 1 2)))) 575 (when (and (mexptp x) 576 (ratnump (caddr x))) 577 (let ((rad (rform x)) 578 (rbase (rform (cadr x))) 579 (expon (caddr x))) 580 (radsubst (ratexpt rbase (cadr expon)) 581 (ratexpt rad (caddr expon)))))) 582 583 584(defun goodform (l) ;;bad -> good 585 (loop for (exp coef) on l by #'cddr 586 collect (cons exp coef))) 587 588(defun factorlogs (l) 589 (prog (negl posl maxpl maxnl maxn) 590 (dolist (log l) 591 (setq log 592 (cons log (goodform 593 (ratfact (rform (radcan1 (cadr log))) 594 #'pfactor)))) 595 (cond ((equal (caadr log) -1) (push log negl)) 596 (t (push log posl)))) 597 (setq negl (flsort negl) posl (flsort posl) l (append negl posl)) 598 (setq negl (mapcar #'cdr negl) 599 posl (mapcar #'cdr posl)) 600 a (setq negl (delete '((-1 . 1)) negl :test #'equal)) 601 (or negl 602 (return (mapc #'(lambda (x) (rplacd x (spc2a (cdr x)))) l))) 603 (setq maxnl (flmaxl negl) 604 maxn (caaar maxnl)) 605 b (setq maxpl (flmaxl posl)) 606 (cond ((and maxpl (flgreat (caaar maxpl) maxn)) 607 (setq posl (flred posl (caaar maxpl))) 608 (go b)) 609 ((and maxpl 610 (not (equal (caaar maxpl) maxn))) 611 (setq maxpl nil))) 612 (cond ((and (flevenp maxpl) (not (flevenp maxnl))) 613 (mapc #'(lambda (fp) (rplaca (car fp) (pminus (caar fp))) 614 (cond ((oddp (cdar fp)) 615 (setq fp (delete '(-1 . 1) fp :test #'equal)) 616 (setq negl (delete fp negl :test #'equal)) 617 (and (cdr fp) (push (cdr fp) posl))))) 618 maxnl) 619 (go a)) 620 (t (setq posl (flred posl maxn) 621 negl (flred negl maxn)) 622 (go a))))) 623 624(defun flevenp (pl) 625 (loop for l in pl never (oddp (cdar l)))) 626 627(defun flred (pl p) 628 (mapl #'(lambda (x) (if (equal p (caaar x)) 629 (rplaca x (cdar x)))) 630 pl) 631 (delete nil pl :test #'equal)) 632 633(defun flmaxl (fpl) ;lists of fac. polys 634 (cond ((null fpl) nil) 635 (t (do ((maxl (list (car fpl)) 636 (cond ((equal (caaar maxl) (caaar ll)) 637 (cons (car ll) maxl)) 638 ((flgreat (caaar maxl) (caaar ll)) maxl) 639 (t (list (car ll))))) 640 (ll (cdr fpl) (cdr ll))) 641 ((null ll) maxl))))) 642 643(defun flsort (fpl) 644 (mapc #'(lambda (x) (rplacd x (sort (cdr x) #'flgreat :key #'car))) 645 fpl)) 646 647(defun nmt (p any) 648 (cond ((pcoefp p) 649 (if (or any (cminusp p)) 1 0)) 650 (t (loop for lp on (p-terms p) by #'cddr 651 sum (nmt (cadr lp) any))))) 652 653(defun nmterms (p) 654 (cond ((equal p -1) (cons 0 0)) 655 (t (cons (nmt p nil) (nmt p t))))) 656 657(defun flgreat (p q) 658 (let ((pn (nmterms p)) (qn (nmterms q))) 659 (cond ((> (car pn) (car qn)) t) 660 ((< (car pn) (car qn)) nil) 661 ((> (cdr pn) (cdr qn)) t) 662 ((< (cdr pn) (cdr qn)) nil) 663 (t (flgreat1 p q))))) 664 665(defun flgreat1 (p q) 666 (cond ((numberp p) 667 (cond ((numberp q) (> p q)) 668 (t nil))) 669 ((numberp q) t) 670 ((pointergp (car p) (car q)) t) 671 ((pointergp (car q) (car p)) nil) 672 ((> (cadr p) (cadr q)) t) 673 ((< (cadr p) (cadr q)) nil) 674 (t (flgreat1 (caddr p) (caddr q))))) 675