src/ratout.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 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))