src/ratout.lisp

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