home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / ratout.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  24.0 KB  |  750 lines
  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. (macsyma-module ratout)
  13.  
  14. ;; THIS IS THE OUT-OF-CORE SEGMENT OF THE RATIONAL FUNCTION PACKAGE. 
  15.  
  16. (DECLARE-TOP
  17.  (SPECIAL $ALGEBRAIC ERRRJFFLAG VARLIST SS *Y* F $FACTORFLAG MODULUS HMODULUS
  18.       GENVAR *A* *ALPHA *VAR* *X* *P *MAX *VAR *RES *CHK *L $INTFACLIM
  19.       $RATFAC U* $RATWTLVL *RATWEIGHTS $RATWEIGHTS $KEEPFLOAT)
  20.  (*LEXPR $RAT)
  21.  (GENPREFIX A_O))
  22.  
  23. (LOAD-MACSYMA-MACROS RATMAC)
  24.  
  25. ;; This splitfile contains Brown's Modular gcd algorithm
  26. (DECLARE-TOP(SPLITFILE MODGCD))
  27.  
  28. (DECLARE-TOP(SPECIAL $GCD XV BIGF1 BIGF2 NONLINDEG $LINHACK
  29.           $INTFACLIM BIGF1TILDE BIGF2TILDE
  30.           GCD $FACTORFLAG *GCDL* LAST-GOOD-PRIME))
  31.  
  32. ;    NEWGCD (X,Y) RETURNS A LIST OF THREE ITEMS,
  33. ;    (GCD, X/GCD, Y/GCD)
  34.  
  35. (DEFUN NEWGCD (X Y MODULUS &AUX HMODULUS)
  36.   (SETQMODULUS MODULUS)
  37.   (LET ((A (COND ((PCOEFP X)
  38.           (COND ((ZEROP X) Y)
  39.             ((PCOEFP Y) (CGCD X Y))
  40.             (T (PCONTENT1 (CDR Y) X))))
  41.          ((PCOEFP Y) (COND ((ZEROP Y) X) (T (PCONTENT1 (CDR X) Y))))
  42.          ((POINTERGP (P-VAR X) (P-VAR Y)) (OLDCONTENT1 (CDR X) Y))
  43.          ((POINTERGP (P-VAR Y) (P-VAR X)) (OLDCONTENT1 (CDR Y) X))
  44.          (T NIL))))
  45.     (COND (A (LIST A (PQUOTIENT X A) (PQUOTIENT Y A)))
  46.       (MODULUS (PGCDP X Y MODULUS))
  47.       (T (PGCDM X Y)))))
  48.  
  49. ;;;***    PMODCONTENT COMPUTES CONTENT OF
  50. ;;;    P IN 
  51. ;;    Z [X ] [X , X , ..., X   ] 
  52. ;;        P  V    1   2        V-1
  53.  
  54. ;;    PMODCONTENT OF 3*A*X IS A, IF MAINVAR IS X (=X )
  55. ;;                              V
  56.  
  57. (DEFUN PMODCONTENT (P)
  58.        (PROG (*VAR *CHK *RES *MAX GCD)
  59.          (SETQ *CHK (CAR P))
  60.          (SETQ *MAX 0)
  61.          (SETQ *VAR (PNEXT (CDR P) NIL))
  62.          (COND ((POINTERGP XV *CHK) (GO RET1))
  63.            ((NULL *VAR) (RETURN (LIST P 1))))
  64.          (PGATH1 (CDR P))
  65.         A    (SETQ *RES 0)
  66.          (PGATH3 (CDR P))
  67.     A2   (COND ((PCOEFP *RES) (COND ((PZEROP *RES) NIL)(T(GO RET1))))
  68.            ((NOT (EQ (CAR *RES) *CHK)) (GO RET1))
  69.            ((NOT (UNIVAR (CDR *RES)))
  70.             (SETQ *RES (CAR (PMODCONTENT *RES)))
  71.             (GO A2))
  72.            (GCD (SETQ GCD (PGCDU GCD *RES)))
  73.            (T (SETQ GCD *RES)))
  74.          (COND ((PCOEFP GCD) (GO RET1))
  75.            ((MINUSP (SETQ *MAX (SUB1 *MAX)))
  76.             (RETURN (LIST GCD (PQUOTIENT P GCD)))))
  77.          (GO A)
  78.     RET1 (RETURN (LIST 1 P))))
  79.  
  80. (DEFUN PGATHERCOEF (P *CHK *RES)
  81.   (IF (NOT (EQ (CAR P) *CHK)) 1 (PGATH2 (CDR P) NIL)))
  82.  
  83. (DEFUN PGATH1 (P)
  84.   (PROG NIL
  85.     (COND ((NULL P) (RETURN *MAX))
  86.           ((PCOEFP (CADR P)) NIL)
  87.           ((EQ (CAADR P) *VAR) (SETQ *MAX (MAX *MAX (CADADR P)))))
  88.     (RETURN (PGATH1 (CDDR P)))))
  89.  
  90. (DEFUN PGATH2 (P VMAX)
  91.        (PROG (V2)
  92.          (COND ((NULL P) (RETURN *RES))
  93.            ((PCOEFP (CADR P)) NIL)
  94.            ((VGREAT (SETQ V2 (PDEGREER (CADR P))) VMAX)
  95.             (SETQ *RES (PSIMP *CHK
  96.                 (LIST (CAR P) (LEADCOEFFICIENT (CADR P)))))
  97.             (SETQ VMAX V2))
  98.            ((EQUAL VMAX V2)
  99.             (SETQ *RES
  100.               (PPLUS *RES
  101.                  (PSIMP *CHK
  102.                     (LIST (CAR P) (LEADCOEFFICIENT (CADR P))))))))
  103.          (RETURN (PGATH2 (CDDR P) VMAX))))
  104.  
  105. (DEFUN PGATH3 (P)
  106.        (PROG (ZZ)
  107.          (COND ((NULL P) (RETURN *RES))
  108.            ((PCOEFP (CADR P))
  109.             (COND ((EQN *MAX 0) (SETQ ZZ (CADR P)) (GO ADD)) (T (GO RET))))
  110.            ((EQ (CAADR P) *VAR) (SETQ ZZ (PTERM (CDADR P) *MAX)) (GO ADD)))
  111.          (COND ((EQN *MAX 0) (SETQ ZZ (CADR P))) (T (GO RET)))
  112.     ADD  (COND ((EQN ZZ 0) (GO RET)))
  113.          (SETQ *RES (PPLUS *RES (PSIMP *CHK (LIST (CAR P) ZZ))))
  114.     RET  (RETURN (PGATH3 (CDDR P)))))
  115.  
  116.  
  117. (DEFUN PNEXT (X *L) (PNEXT1 X) (COND ((NULL *L) NIL)
  118.                      (T (CAR (SORT *L #'POINTERGP)))))
  119.  
  120. (DEFUN PNEXT1 (X) (PROG NIL
  121.             (COND ((NULL X) (RETURN *L))
  122.                   ((OR (PCOEFP (CADR X)) (MEMQ (CAADR X) *L)) NIL)
  123.                   (T (SETQ *L (CONS (CAADR X) *L))))
  124.             (RETURN (PNEXT1 (CDDR X)))))
  125.  
  126. (DEFUN VGREAT (X Y) (COND ((NULL X) NIL)
  127.               ((NULL Y) T)
  128.               ((POINTERGP (CAR X)(CAR Y))T)
  129.               ((NOT (EQ (CAR X)(CAR Y)))NIL)
  130.               ((GREATERP (CADR X)(CADR Y)) T)
  131.               ((EQN (CADR X)(CADR Y))(VGREAT (CDDR X)(CDDR Y)))
  132.               (T NIL)))
  133.  
  134. (DEFUN PDEGREER (X)
  135.   (IF (PCOEFP X) () (CONS (CAR X) (CONS (CADR X) (PDEGREER (CADDR X))))))
  136.  
  137.  
  138. ;;***    PGCDP CORRESPONDS TO BROWN'S ALGORITHM P
  139.  
  140. (DEFUN PGCDP (BIGF1 BIGF2 MODULUS)
  141.        (PROG (C C1        C2        N        Q
  142.         H1TILDE        H2TILDE        GSTAR        H1STAR
  143.         H2STAR        XV        E        B
  144.         GBAR        NUBAR        NU1BAR        NU2BAR
  145.         GTILDE        F1TILDE        F2TILDE        BIGGTILDE
  146.         DEGREE        F1        F1F2        HMODULUS)
  147.          (SETQMODULUS MODULUS)
  148.          (COND ((AND (UNIVAR (CDR BIGF1)) (UNIVAR (CDR BIGF2)))
  149.             (SETQ Q (PGCDU BIGF1 BIGF2))
  150.             (RETURN (LIST Q (PQUOTIENT BIGF1 Q) (PQUOTIENT BIGF2 Q)))))
  151.          (SETQ XV (CAR BIGF1))
  152.          (SETQ BIGF1 (PMODCONTENT BIGF1))
  153.          (SETQ BIGF2 (PMODCONTENT BIGF2))
  154.          (SETQ C (PGCDU (SETQ C1 (CAR BIGF1)) (SETQ C2 (CAR BIGF2))))
  155.          (SETQ BIGF1 (CADR BIGF1))
  156.          (SETQ BIGF2 (CADR BIGF2))
  157.          (SETQ N 0)
  158.          (SETQ E (PDEGREER BIGF2))
  159.          (SETQ DEGREE (PDEGREER BIGF1))
  160.          (COND ((VGREAT E DEGREE) (SETQ E DEGREE)))
  161.          (SETQ B (LSH MODULUS -1))
  162.          (SETQ GBAR
  163.            (PGCDU (SETQ F1 (PGATHERCOEF BIGF1 XV 0))
  164.               (SETQ F1F2
  165.                 (PGATHERCOEF BIGF2 XV 0))))
  166.          (COND ((EQUAL 0 F1F2) (GO STEP15A)))
  167.          (SETQ NUBAR (PDEGREE GBAR XV))
  168.          (SETQ NU1BAR (f+ NUBAR (PDEGREE BIGF1 XV)))
  169.          (SETQ NU2BAR (f+ NUBAR (PDEGREE BIGF2 XV)))
  170.          (SETQ F1F2 (PTIMES F1 F1F2))
  171.          (SETQ NUBAR (MAX NU1BAR NU2BAR))
  172.     STEP6(SETQ B (CPLUS B 1))
  173.          (COND ((EQUAL (pcsubst F1F2 B XV) 0) (GO STEP6)))
  174.          ;; Step 7
  175.          (SETQ GTILDE (pcsubst GBAR B XV))
  176.          (SETQ F1TILDE (pcsubst BIGF1 B XV))
  177.          (SETQ F2TILDE (pcsubst BIGF2 B XV))
  178.          (SETQ BIGGTILDE
  179.            (PTIMESCHK GTILDE
  180.                   (CAR (SETQ H2TILDE (NEWGCD F1TILDE F2TILDE MODULUS)))))
  181.          (COND ((PCOEFP BIGGTILDE) (GO STEP15A)))
  182.           (SETQ H1TILDE (CADR H2TILDE))
  183.              (SETQ H2TILDE (CADDR H2TILDE))
  184.          (SETQ DEGREE (PDEGREER BIGGTILDE))
  185.          (COND ((VGREAT DEGREE E) (GO STEP6))
  186.            ((VGREAT E DEGREE) (SETQ N 0) (SETQ E DEGREE)))
  187.          (SETQ N (ADD1 N))
  188.          (COND ((EQUAL N 1) (SETQ Q (LIST XV 1 1 0 (CMINUS B)))
  189.                 (SETQ GSTAR BIGGTILDE)
  190.                 (SETQ H1STAR H1TILDE)
  191.                 (SETQ H2STAR H2TILDE))
  192.            (T (SETQ GSTAR (LAGRANGE33 GSTAR BIGGTILDE Q B))
  193.               (SETQ H1STAR (LAGRANGE33 H1STAR H1TILDE Q B))
  194.               (SETQ H2STAR (LAGRANGE33 H2STAR H2TILDE Q B))
  195.               (SETQ Q (PTIMES Q (LIST XV 1 1 0 (CMINUS B))))))
  196.          ;; Step 12
  197.          (COND ((NOT (> N NUBAR)) (GO STEP6)))
  198.          ;; Step 13
  199.          (COND ((OR (NOT (= NU1BAR (f+ (SETQ DEGREE (PDEGREE GSTAR XV))
  200.                          (PDEGREE H1STAR XV))))
  201.             (NOT (= NU2BAR (f+ DEGREE (PDEGREE H2STAR XV)))))
  202.             (SETQ N 0)
  203.             (GO STEP6)))
  204.          (SETQ GSTAR (CADR (PMODCONTENT GSTAR)))
  205.          ;; Step 15
  206.          (SETQ Q (PGATHERCOEF GSTAR XV 0))
  207.          (RETURN (MONICGCD  (PTIMESCHK C GSTAR)
  208.                 (PTIMESCHK (PQUOTIENT C1 C) (PQUOTIENTCHK H1STAR Q))
  209.                 (PTIMESCHK (PQUOTIENT C2 C) (PQUOTIENTCHK H2STAR Q))
  210.                 (LEADCOEFFICIENT GSTAR)))
  211.     STEP15A
  212.          (RETURN (LIST C
  213.             (PTIMESCHK (PQUOTIENT C1 C) BIGF1)
  214.             (PTIMESCHK (PQUOTIENT C2 C) BIGF2))) ))
  215.  
  216.  
  217. (DEFUN MONICGCD (GCD X Y LCF)
  218.     (COND ((EQN LCF 1) (LIST GCD X Y))
  219.           (T (LIST    (PTIMES (CRECIP LCF) GCD)
  220.             (PTIMES LCF X)
  221.             (PTIMES LCF Y) )) ))
  222.  
  223.  
  224. ;***    PGCDM CORRESPONDS TO BROWN'S ALGORITHM M
  225.  
  226.  
  227. (DEFUN PGCDM
  228.        (BIGF1 BIGF2)
  229.        (PROG (C C1        C2        F1        F2    N
  230.         E        DEGREE        MUBAR        P
  231.         NONLINDEG    GTILDE        H1TILDE        H2TILDE
  232.         MODULUS        HMODULUS    BIGF1TILDE    BIGF2TILDE
  233.         BIGGTILDE    Q        H1STAR        H2STAR
  234.         GSTAR        XV              GBAR)
  235.          (SETQ P *ALPHA)
  236.          (SETQ XV (CAR BIGF1))
  237.          ;; Step 1
  238.          (SETQ F1 (PCONTENT BIGF1))
  239.          (SETQ F2 (PCONTENT BIGF2))
  240.          (SETQ C (CGCD (SETQ C1 (CAR F1)) (SETQ C2 (CAR F2))))
  241.          (SETQ BIGF1 (CADR F1))
  242.          (SETQ BIGF2 (CADR F2))
  243.          ;; Step 3
  244.          (SETQ F1 (LEADCOEFFICIENT BIGF1))
  245.          (SETQ F2 (LEADCOEFFICIENT BIGF2))
  246.          (SETQ GBAR (CGCD F1 F2))
  247.          ;; Step 4
  248.          (SETQ N 0)
  249.          (SETQ DEGREE (PDEGREER BIGF1))
  250.          (SETQ E (PDEGREER BIGF2))
  251.          (COND ((VGREAT E DEGREE) (SETQ E DEGREE)))
  252.          ;; Step 5
  253.          (SETQ MUBAR
  254.            (TIMES 2 GBAR (MAX (MAXCOEFFICIENT BIGF1)
  255.                       (MAXCOEFFICIENT BIGF2))))
  256.          (GO STEP6A)
  257.     STEP6(SETQ P (NEWPRIME P))
  258.     STEP6A
  259.          (COND ((OR (EQUAL 0 (REMAINDER F1 P)) (EQUAL 0 (REMAINDER F2 P)))
  260.             (GO STEP6)))
  261.          (SETQMODULUS P)
  262.          ;; Step 7
  263.          (SETQ GTILDE (PMOD GBAR))
  264.          ;; Step 8
  265.          (SETQ BIGGTILDE
  266.            (PTIMESCHK GTILDE
  267.                   (CAR (SETQ H2TILDE
  268.                      (NEWGCD (PMOD BIGF1) (PMOD BIGF2)
  269.                          MODULUS)))))
  270.          (COND ((PCOEFP BIGGTILDE) (SETQ MODULUS NIL)
  271.                        (SETQ GSTAR 1)
  272.                        (SETQ H1STAR BIGF1)
  273.                        (SETQ H2STAR BIGF2)
  274.                        (GO STEP15)))
  275.          (COND ((NULL (CDR H2TILDE))
  276.             (SETQ H1TILDE (PQUOTIENT (PMOD BIGF1) (CAR H2TILDE)))
  277.             (SETQ H2TILDE (PQUOTIENT (PMOD BIGF2) (CAR H2TILDE))))
  278.            (T (SETQ H1TILDE (CADR H2TILDE))
  279.               (SETQ H2TILDE (CADDR H2TILDE))))
  280.          (SETQ DEGREE (PDEGREER BIGGTILDE))
  281.          (COND ((VGREAT DEGREE E) (GO STEP6))
  282.            ((VGREAT E DEGREE) (SETQ N 0) (SETQ E DEGREE)))
  283.          (SETQ N (ADD1 N))
  284.          ;; Step 11
  285.          (SETQMODULUS NIL)
  286.          (COND ((EQUAL N 1) (SETQ Q P)
  287.                 (SETQ GSTAR BIGGTILDE)
  288.                 (SETQ H1STAR H1TILDE)
  289.                 (SETQ H2STAR H2TILDE))
  290.            (T (SETQ GSTAR (LAGRANGE3 GSTAR BIGGTILDE P Q))
  291.               (SETQ H1STAR (LAGRANGE3 H1STAR H1TILDE P Q))
  292.               (SETQ H2STAR (LAGRANGE3 H2STAR H2TILDE P Q))
  293.               (SETQ Q (TIMES P Q))))
  294.          ;; Step 12
  295.          (COND ((GREATERP MUBAR Q) (GO STEP6)))
  296.          (COND ((GREATERP (TIMES 2 (MAX
  297.     (TIMES (SETQ GTILDE (NORM GSTAR)) (MAXCOEFFICIENT H1STAR))
  298.     (TIMES GTILDE (MAXCOEFFICIENT H2STAR)) ))
  299.                   Q)
  300.             (GO STEP6)))
  301.          (SETQMODULUS NIL)
  302.          (SETQ GSTAR (CADR (PCONTENT GSTAR)))
  303.     STEP15
  304.          (SETQ LAST-GOOD-PRIME P)
  305.              (SETQ Q (LEADCOEFFICIENT GSTAR))
  306.          (RETURN (LIST (PTIMESCHK C GSTAR)
  307.                (PTIMESCHK (CQUOTIENT C1 C) (PQUOTIENTCHK H1STAR Q))
  308.                (PTIMESCHK (CQUOTIENT C2 C) (PQUOTIENTCHK H2STAR Q))))))
  309.  
  310. ;    THE FUNCTIONS ON THIS PAGE ARE USED BY KRONECKER FACTORING
  311.  
  312. (DECLARE-TOP(SPLITFILE KRONEC))
  313.  
  314. (DEFUN PKRONECK (P) (PROG (MAXEXP I L *P FACTORS FACTOR ERRRJFFLAG)
  315.                (SETQ MAXEXP (QUOTIENT (CADR P) 2))
  316.                (SETQ I 1)
  317.           A    (COND ((GREATERP I MAXEXP) (RETURN (CONS P FACTORS))))
  318.                (SETQ L (P1 (REVERSE ((LAMBDA (P I $FACTORFLAG)
  319.                              (PFACTOR2 P I))
  320.                          P
  321.                          I
  322.                          T))))
  323.           B    (COND ((NULL L) (GO D)))
  324.                (SETQ *L (CAR L))
  325.                (SETQ *P (CAR P))
  326.                (SETQ ERRRJFFLAG T)
  327.                (SETQ FACTOR (ERRSET (PINTERPOLATE *L *P) NIL))
  328.                (SETQ ERRRJFFLAG NIL)
  329.                (SETQ L (CDR L))
  330.                (COND ((ATOM FACTOR) (GO B))
  331.                  (T (SETQ FACTOR (CAR FACTOR))))
  332.                (COND ((OR (PCOEFP FACTOR)
  333.                   (NOT (EQN (CAR P) (CAR FACTOR)))
  334.                   (NOT (PZEROP (PREM P FACTOR))))
  335.                   (GO B)))
  336.                (COND (MODULUS (PMONICIZE (CDR FACTOR)))
  337.                  ((PMINUSP FACTOR) (SETQ FACTOR (PMINUS FACTOR))))
  338.                (SETQ P (PQUOTIENT P FACTOR))
  339.                (SETQ MAXEXP (QUOTIENT (CADR P) 2))
  340.                (SETQ FACTORS (CONS FACTOR FACTORS))
  341.                (COND ((OR (EQN P 1) (EQN P -1)) (RETURN FACTORS)))
  342.                (GO A)
  343.           D    (SETQ I (ADD1 I))
  344.                (GO A)
  345.      ))
  346.  
  347. (DEFUN PFACTOR2 (P I) (COND ((LESSP I 0.) NIL)
  348.                    (T (CONS (PFACTOR (PCSUBST P I (CAR P)))
  349.                         (PFACTOR2 P (SUB1 I))))))
  350.  
  351. (DEFUN POWERSET (X N) (COND ((NULL X) (QUOTE (1 NIL)))
  352.                  ((EQUAL X 1) (QUOTE (1)))
  353.                  (T (CONS 1 (PTTS1 X N X)))))
  354.      
  355.  
  356. (DEFUN ALLPRODS (X Y) (COND ((NULL X) NIL)
  357.                  ((NULL Y) NIL)
  358.                  (T (NCONC (AP1 (CAR X) Y) (ALLPRODS (CDR X) Y)))))
  359.      
  360. (DEFUN AL1 (F R LEN)
  361.   (PROG (SS)
  362.     (COND
  363.      ((EQUAL LEN 1)
  364.       (RETURN (MAPCAR #'(LAMBDA (*Y*) (CONS *Y* NIL)) F)))
  365.      ((NULL R) (RETURN NIL))
  366.      (T
  367.       (MAPC #'(LAMBDA (*Y*)
  368.               (SETQ SS
  369.                 (NCONC SS
  370.                        (MAPCAR #'(LAMBDA (Z) (CONS Z *Y*))
  371.                            F))))
  372.         (AL1 (CAR R) (CDR R) (SUB1 LEN)))
  373.       (RETURN SS)))))
  374.  
  375.  
  376. (DEFUN AP1 (X L) (COND ((NULL L) NIL)
  377.                (T (CONS (PTIMES X (CAR L)) (AP1 X (CDR L))))))
  378.  
  379. (DEFUN PTTS1 (X N Y) (COND ((EQN N 1) (LIST Y))
  380.                (T (CONS Y (PTTS1 X (SUB1 N) (PTIMES X Y))))))
  381.  
  382. (DEFUN P1 (L) (PROG (A)
  383.             (SETQ A (MAPCAR #'P11 L))
  384.             (RETURN (COND ((NULL L) NIL)
  385.                   (T (CDR (AL1 (CAR A)
  386.                            (CDR A)
  387.                            (LENGTH A))))))))
  388.      
  389. (DEFUN P11 (ELE) (COND ((NULL (CDDR ELE)) (POWERSET (CAR ELE) (CADR ELE)))
  390.                (T (ALLPRODS (POWERSET (CAR ELE) (CADR ELE))
  391.                     (P11 (CDDR ELE))))))
  392.  
  393. (DEFUN PINTERPOLATE (L VAR)
  394.        (PSIMP VAR (PINTERPOLATE1 (PINTERPOLATE2 L 1)
  395.                  (DIFFERENCE (LENGTH L) 2))))
  396.  
  397. (DEFUN PINTERPOLATE1 (X N)
  398.        (PINTERPOLATE4 (PINTERPOLATE5 (REVERSE X) 1 N N) (ADD1 N)))
  399.      
  400. (DEFUN PINTERPOLATE2 (X N)
  401.        (COND ((NULL (CDR X)) X)
  402.          (T (CONS (CAR X)
  403.               (PINTERPOLATE2 (PINTERPOLATE3 X N) (ADD1 N))))))
  404.  
  405. (DEFUN PINTERPOLATE3 (X N)
  406.        (COND ((NULL (CDR X)) NIL)
  407.          (T (CONS (PQUOTIENT (PDIFFERENCE (CADR X) (CAR X)) N)
  408.               (PINTERPOLATE3 (CDR X) N)))))
  409.      
  410. (DEFUN PINTERPOLATE4 (X N)
  411.        (COND ((NULL X) NIL)
  412.          ((PZEROP (CAR X)) (PINTERPOLATE4 (CDR X) (SUB1 N)))
  413.          (T (CONS N (CONS (CAR X)
  414.                   (PINTERPOLATE4 (CDR X) (SUB1 N)))))))
  415.  
  416. (DEFUN PINTERPOLATE5 (X I J N)
  417.        (COND ((GREATERP I N) X)
  418.          (T (PINTERPOLATE5 (CONS (CAR X) (PINTERPOLATE6 X I J))
  419.                    (ADD1 I)
  420.                    (SUB1 J)
  421.                    N))))
  422.      
  423. (DEFUN PINTERPOLATE6 (X I J)
  424.        (COND ((ZEROP I) (CDR X))
  425.          (T (CONS (PDIFFERENCE (CADR X) (PCTIMES J (CAR X)))
  426.               (PINTERPOLATE6 (CDR X) (SUB1 I) J)))))
  427.  
  428.  
  429. (declare-top (SPLITFILE FASTT))
  430.  
  431. ;; THE N**(1.585) MULTIPLICATION SCHEME
  432. ;;FOLLOWS.  IT SHOULD BE USED ONLY WHEN BOTH INPUTS ARE MULTIVARIATE,
  433. ;;DENSE, AND OF NEARLY THE SAME SIZE.  OR ABSOLUTELY TREMENDOUS.
  434. ;;(THE CLASSICAL MULTIPLICATION SCHEME IS N**2 WHERE N IS SIZE OF
  435. ;;POLYNOMIAL   (OR N*M FOR DIFFERENT SIZES).  FOR THIS
  436. ;;CASE, N IS APPX. THE SIZE OF LARGER.
  437.  
  438. (DEFMFUN $FASTTIMES (X Y)
  439.   (COND ((AND (NOT (ATOM X)) (NOT (ATOM Y))
  440.           (EQUAL (CAR X) (CAR Y)) (EQUAL (CAAR X) 'MRAT)
  441.           (EQUAL (CDDR X) 1) (EQUAL (CDDR Y) 1))
  442.      (CONS (CAR X)(CONS (FPTIMES (CADR X)(CADR Y))1)))
  443.     (T (merror "Use FASTTIMES only on CRE polynomials with same varlists"))))
  444.  
  445. (DEFUN FPTIMES (X Y) (COND ((OR (PZEROP X) (PZEROP Y)) (PZERO))
  446.                ((PCOEFP X) (PCTIMES X Y))
  447.                ((PCOEFP Y) (PCTIMES Y X))
  448.                ((EQ (CAR X) (CAR Y))
  449.                 (COND((OR(UNIVAR(CDR X))(UNIVAR(CDR Y)))
  450.                   (CONS (CAR X) (PTIMES1 (CDR X) (CDR Y))))
  451.                  (T(CONS (CAR X) (FPTIMES1 (CDR X)(CDR Y))))))
  452.                ((POINTERGP (CAR X) (CAR Y))
  453.                 (CONS (CAR X) (PCTIMES1 Y (CDR X))))
  454.                (T (CONS (CAR Y) (PCTIMES1 X (CDR Y))))))
  455.      
  456. (DEFUN FPTIMES1 (F G)
  457.        (PROG (A B C D)
  458.          (COND ((OR (NULL F) (NULL G)) (RETURN NIL))
  459.            ((NULL (CDDR F))
  460.             (RETURN (LSFT (PCTIMES1 (CADR F) G) (CAR F))))
  461.            ((NULL (CDDR G))
  462.             (RETURN (LSFT (PCTIMES1 (CADR G) F) (CAR G)))))
  463.          (SETQ D (LSH (ADD1 (MAX (CAR F) (CAR G))) -1))
  464.          (SETQ F (HALFSPLIT F D) G (HALFSPLIT G D))
  465.          (SETQ A (FPTIMES1 (CAR F) (CAR G)))
  466.          (SETQ B
  467.            (FPTIMES1 (PPLUS1 (CAR F) (CDR F)) (PPLUS1 (CAR G) (CDR G))))
  468.          (SETQ C (FPTIMES1 (CDR F) (CDR G)))
  469.          (SETQ B (PDIFFER1 (PDIFFER1 B A) C))
  470.          (RETURN (PPLUS1 (LSFT A (LSH D 1)) (PPLUS1 (LSFT B D) C)))))
  471.  
  472. (DEFUN HALFSPLIT (P D)
  473.        (DO ((A) (P P (CDDR P)))
  474.        ((OR (NULL P) (< (CAR P) D)) (CONS (NREVERSE A) P))
  475.        (SETQ A (CONS (CADR P) (CONS (f- (CAR P) D) A)))))
  476.  
  477. (DEFUN LSFT (P N)
  478.        (DO ((Q P (CDDR (RPLACA Q (f+ (CAR Q) N))))) ((NULL Q)))
  479.        P)
  480.  
  481.  
  482. (declare-top (SPLITFILE RATWT)
  483.      (SPECIAL WTSOFAR XWEIGHT $RATWTLVL V *X* *I*)
  484.      (FIXNUM *I* XWEIGHT WTSOFAR XWT (PWEIGHT NOTYPE))) 
  485.  
  486. ;;; TO TRUNCATE ON E, DO RATWEIGHT(E,1);
  487. ;;;THEN DO RATWTLVL:N.  ALL POWERS >N GO TO 0.
  488.  
  489. (DEFMFUN $RATWEIGHT N 
  490.  (COND ((ODDP N) (MERROR "RATWEIGHT takes an even number of arguments.")))
  491.  (DO ((*I* 1 (f+ *I* 2))) ((> *I* N))
  492.      (RPLACD (or (zl-ASSOC (ARG *I*) *RATWEIGHTS)
  493.          (CAR (SETQ *RATWEIGHTS (CONS (LIST (ARG *I*)) *RATWEIGHTS))))
  494.          (ARG (f1+ *I*)))) 
  495.  (SETQ $RATWEIGHTS (CONS '(MLIST SIMP) (DOT2L *RATWEIGHTS)))
  496.  (COND ((= N 0) $RATWEIGHTS) (T (CONS '(MLIST) (LISTIFY N)))))
  497.  
  498. (DEFUN PWEIGHT (X) (OR (GET X '$RATWEIGHT) 0.)) 
  499.  
  500. (DEFUN WTPTIMES (X Y WTSOFAR) 
  501.        (COND ((OR (PZEROP X) (PZEROP Y) (> WTSOFAR $RATWTLVL))
  502.           (PZERO))
  503.          ((PCOEFP X) (WTPCTIMES X Y))
  504.          ((PCOEFP Y) (WTPCTIMES Y X))
  505.          ((EQ (CAR X) (CAR Y))
  506.           (PALGSIMP (CAR X)
  507.              (WTPTIMES1 (CDR X)
  508.                 (CDR Y)
  509.                 (PWEIGHT (CAR X)))
  510.             (ALG X)))
  511.          ((POINTERGP (CAR X) (CAR Y))
  512.           (PSIMP (CAR X)
  513.              (WTPCTIMES1 Y (CDR X) (PWEIGHT (CAR X)))))
  514.          (T (PSIMP (CAR Y)
  515.                (WTPCTIMES1 X (CDR Y) (PWEIGHT (CAR Y))))))) 
  516.  
  517. (DEFUN WTPTIMES1 (*X* Y XWEIGHT) 
  518.        (PROG (U* V)
  519.          (declare (special v))
  520.          (SETQ V (SETQ U* (WTPTIMES2 Y)))
  521.     A    (SETQ *X* (CDDR *X*))
  522.          (COND ((NULL *X*) (RETURN U*)))
  523.          (WTPTIMES3 Y)
  524.          (GO A)))
  525.  
  526.  
  527. (DEFUN WTPTIMES2 (Y) 
  528.        (COND ((NULL Y) NIL)
  529.          (T ((LAMBDA (II) (DECLARE (FIXNUM II))
  530.                  (COND ((> II $RATWTLVL) (WTPTIMES2 (CDDR Y)))
  531.                    (T (PCOEFADD (f+ (CAR *X*) (CAR Y))
  532.                         (WTPTIMES (CADR *X*) (CADR Y) II)
  533.                         (WTPTIMES2 (CDDR Y))))))
  534.          (f+ (f* XWEIGHT (f+ (CAR *X*) (CAR Y))) WTSOFAR))))) 
  535.  
  536. (DEFUN WTPTIMES3 (Y) 
  537.        (PROG ((E 0) U C)
  538.     (DECLARE (FIXNUM E) (special v))
  539.  
  540.     A1   (COND ((NULL Y) (RETURN NIL)))
  541.          (SETQ E (f+ (CAR *X*) (CAR Y)))
  542.          (SETQ C (WTPTIMES (CADR Y) (CADR *X*) (f+ WTSOFAR (f* XWEIGHT E))))
  543.          (COND ((PZEROP C) (SETQ Y (CDDR Y)) (GO A1))
  544.            ((OR (NULL V) (> E (CAR V))) (SETQ U* (SETQ V (PPLUS1 U* (LIST E C)))) (SETQ Y (CDDR Y)) (GO A1))
  545.            ((EQN E (CAR V))
  546.             (SETQ C (PPLUS C (CADR V)))
  547.             (COND ((PZEROP C) (SETQ U* (SETQ V (PDIFFER1 U* (LIST (CAR V) (CADR V)))))) (T (RPLACA (CDR V) C)))
  548.             (SETQ Y (CDDR Y))
  549.             (GO A1)))
  550.     A    (COND ((AND (CDDR V) (> (CADDR V) E)) (SETQ V (CDDR V)) (GO A)))
  551.          (SETQ U (CDR V))
  552.     B    (COND ((OR (NULL (CDR U)) (< (CADR U) E)) (RPLACD U (CONS E (CONS C (CDR U)))) (GO E)))
  553.          (COND ((PZEROP (SETQ C (PPLUS (CADDR U) C))) (RPLACD U (CDDDR U)) (GO D)) (T (RPLACA (CDDR U) C)))
  554.     E    (SETQ U (CDDR U))
  555.     D    (SETQ Y (CDDR Y))
  556.          (COND ((NULL Y) (RETURN NIL))
  557.                ((PZEROP
  558.              (SETQ C (WTPTIMES (CADR *X*) (CADR Y)
  559.                        (f+ WTSOFAR (f* XWEIGHT
  560.                              (SETQ E (f+ (CAR *X*) (CAR Y))))))))
  561.             (GO D)))
  562.     C    (COND ((AND (CDR U) (> (CADR U) E)) (SETQ U (CDDR U)) (GO C)))
  563.          (GO B))) 
  564.  
  565.  
  566. (DEFUN WTPCTIMES (C P) 
  567.          (COND ((PCOEFP P) (CTIMES C P))
  568.                (T (PSIMP (CAR P) (WTPCTIMES1 C (CDR P) (PWEIGHT (CAR P)))))))
  569.  
  570. (DEFUN WTPCTIMES1 (C X XWT) 
  571.   (PROG (CC) 
  572.     (RETURN
  573.      (COND ((NULL X) NIL)
  574.            (T (SETQ CC (WTPTIMES C
  575.                      (CADR X)
  576.                      (f+ WTSOFAR (f* XWT (CAR X)))))
  577.           (COND ((PZEROP CC) (WTPCTIMES1 C (CDDR X) XWT))
  578.             (T (CONS (CAR X)
  579.                  (CONS CC
  580.                        (WTPCTIMES1 C
  581.                            (CDDR X)
  582.                            XWT))))))))))
  583.  
  584. (DEFUN WTPEXPT (X N) (COND ((= N 0) 1) ((= N 1) X) (T (WTPTIMES X (WTPEXPT X (f1- N)) 0))))
  585.  
  586.  
  587. (declare-top (SPLITFILE HORNER))
  588.  
  589. (DEFMFUN $HORNER NARGS
  590.  (DECLARE (FIXNUM NARGS))
  591.  (IF (= NARGS 0) (WNA-ERR '$HORNER))
  592.  (LET (($RATFAC NIL) (VARLIST (CDR $RATVARS)) GENVAR (X NIL)
  593.        (ARG1 (TAYCHK2RAT (ARG 1)))
  594.        (L (CDR (LISTIFY NARGS))))
  595.       (COND ((MBAGP ARG1)
  596.          (CONS (CAR ARG1)
  597.            (MAPCAR #'(LAMBDA (U) (APPLY '$HORNER (CONS U L)))
  598.                (CDR ARG1))))
  599.         (T (SETQ X (APPLY '$RAT (CONS ARG1 L)))
  600.            (MAPC #'(LAMBDA (Y Z) (PUTPROP Y Z 'DISREP))
  601.              (CADDDR (CAR X))
  602.              (CADDAR X))
  603.            (DIV* (HORNREP (CADR X)) (HORNREP (CDDR X)))))))
  604.  
  605. (DEFUN HORNREP (P) (IF (PCOEFP P) P (HORN+ (CDR P) (GET (CAR P) 'DISREP))))
  606.  
  607. (DEFUN HORN+ (L VAR)
  608.        (PROG (ANS LAST)
  609.          (SETQ ANS (HORNREP (CADR L)))
  610.        A (SETQ LAST (CAR L) L (CDDR L))
  611.          (COND ((NULL L)
  612.             (RETURN (COND ((EQUAL LAST 0) ANS)
  613.                   (T (LIST '(MTIMES)
  614.                        (LIST '(MEXPT) VAR LAST) ANS)))))
  615.            (T (SETQ ANS (LIST '(MPLUS)
  616.                       (HORNREP (CADR L))
  617.                       (LIST '(MTIMES) 
  618.                         (LIST '(MEXPT) VAR (DIFFERENCE LAST (CAR L)))
  619.                         ANS)))))
  620.          (GO A)))
  621.  
  622.  
  623. (declare-top (SPLITFILE PFRAC)
  624.      (SPECIAL Y RISCHPF GENVAR $SAVEFACTORS CHECKFACTORS W
  625.           EXP VAR X $FACTORFLAG $RATFAC
  626.           $KEEPFLOAT RATFORM ROOTFACTOR 
  627.           WHOLEPART PARNUMER VARLIST N))
  628.  
  629. (declare-top(*lexpr partfrac))
  630.  
  631. (DEFMFUN $PARTFRAC (EXP VAR) 
  632.        (COND ((AND (NOT (ATOM EXP)) (MEMQ (CAAR EXP) '(MEQUAL MLIST $MATRIX)))
  633.           (CONS (CAR EXP) (MAPCAR (FN (U) ($PARTFRAC U VAR)) (CDR EXP))))
  634.          ((AND (ATOM VAR) (NOT (AMONG VAR EXP))) EXP)
  635.          (T (LET (($SAVEFACTORS T) (CHECKFACTORS ()) (VARLIST (LIST VAR))
  636.               $RATFAC $ALGEBRAIC RATFORM GENVAR)
  637.              (DESETQ (RATFORM . EXP) (TAYCHK2RAT EXP))
  638.              (SETQ VAR (CAADR (RATF VAR)))
  639.              (SETQ EXP (PARTFRAC EXP VAR))
  640.              (SETQ EXP (CONS (CAR EXP)        ;FULL DECOMP?
  641.                      (MAPCAN #'PARTFRACA (CDR EXP))))
  642.              (ADD2* (DISREP (CAR EXP))
  643.                 (CONS '(MPLUS) 
  644.                   (MAPCAR
  645.                    (FN (L)
  646.                        (LET (((COEF POLY EXP) L))
  647.                         (LIST '(MTIMES)
  648.                           (DISREP  COEF)
  649.                           (LIST '(MEXPT)
  650.                             (DISREP POLY)
  651.                             (MINUS EXP)))))
  652.                    (CDR EXP))))))))
  653.  
  654. (defun partfraca (llist)
  655.        (let (((coef poly exp) llist))
  656.         (do ((nc (ratdivide coef poly) (ratdivide (car nc) poly))
  657.          (n exp (f1- n))
  658.          (ans))
  659.         ((rzerop (car nc)) (cons (list (cdr nc) poly n) ans))
  660.         (push (list (cdr nc) poly n) ans))))
  661.  
  662. (defun partfrac (rat var &OPTIONAL facdenom)        
  663.    (let* (((wholepart frpart) (pdivide (car rat) (cdr rat)))
  664.       ((num . denom) (ratqu frpart (cdr rat))))
  665.       (cond ((pzerop num) (cons wholepart nil))
  666.         ((or (pcoefp denom) (pointergp var (car denom))) (cons rat nil))
  667.         (t (let (((content bpart) (oldcontent denom)))
  668.             (do ((factor (or facdenom (pfactor bpart)) (cddr factor))
  669.              (apart) (y) (parnumer))
  670.             ((null factor) (cons wholepart parnumer))
  671.             (cond
  672.              ((zerop (pdegree (car factor) var)))
  673.              (t (setq apart (pexpt (car factor) (cadr factor))
  674.                   bpart (pquotient bpart apart)
  675.                   y (bprog apart bpart)
  676.                   frpart (cdr (ratdivide (ratti num (cdr y) t)
  677.                              apart)))
  678.                 (push (list (ratqu frpart content)
  679.                     (car factor)
  680.                     (cadr factor))
  681.                   parnumer)
  682.                 (desetq (num . content)
  683.                     (cdr (ratdivide (ratqu (ratti num (car y) t)
  684.                                content)
  685.                             bpart)))))))))))
  686.  
  687. (declare-top(unspecial exp f n ss v var w xv y
  688.             *a* *chk *l *max *p
  689.             *res u* *var* *x* *y*))
  690.  
  691. ;; $RATDIFF TAKES DERIVATIVES FAST.  IT ASSUMES THAT THE
  692. ;; ONLY ENTITY WHICH DEPENDS ON X IS X ITSELF.
  693. ;; THAT IS, DEPENDENCIES DECLARED EXPLICITLY OR IMPLICITLY ARE
  694. ;; TOTALLY IGNORED.  RATDIFF(F(X),X) IS 0.  RATDIFF(Y,X) IS 0.
  695. ;; ANY OTHER USAGE MUST GO THROUGH $DIFF.
  696. ;; FURTHERMORE, X IS ASSUMED TO BE AN ATOM OR A SINGLE ITEM ON
  697. ;; VARLIST.  E.G. X MIGHT BE SIN(U), BUT NOT 2*SIN(U).
  698.  
  699. (DECLARE-TOP(SPLITFILE RATDIF) (SPECIAL VARLIST GENVAR X))
  700.  
  701. (DEFMFUN $RATDIFF (P X)
  702.        (IF ($RATP P)
  703.        (SETQ P (MINIMIZE-VARLIST
  704.             (IF (MEMQ 'TRUNC (CDAR P)) ($TAYTORAT P) P))))
  705.        (LET ((FORMFLAG ($RATP P)) (VARLIST) (GENVAR))
  706.         (NEWVAR X) (NEWVAR P)
  707.         (OR (ANDMAPC #'(lambda (EXP)
  708.                  (OR (ALIKE1 X EXP) (FREE EXP X)))
  709.              VARLIST)
  710.         (MERROR "RATDIFF variable is embedded in kernel"))
  711.         (SETQ P (RATF P))
  712.         (SETQ X (CAADR (RATF X)))
  713.         (SETQ P (CONS (CAR P) (RATDERIVATIVE (CDR P) X)))
  714.         (IF FORMFLAG P ($RATDISREP P))))
  715.  
  716. (DECLARE-TOP(UNSPECIAL X))
  717.  
  718. (declare-top (SPLITFILE PFET) (SPECIAL $PFEFORMAT VARLIST $FACTORFLAG M V DOSIMP))
  719.  
  720. (DEFMFUN $PFET (M) 
  721.        (PROG (LISTOV $PFEFORMAT VARLIST $FACTORFLAG) 
  722.          (SETQ $PFEFORMAT T)
  723.          (NEWVAR M)
  724.          (SETQ LISTOV VARLIST)
  725.          (MAPC #'(LAMBDA (R) (SETQ M (PFET1 M R)))
  726.            LISTOV)
  727.          (SETQ M (SIMPLIFY M))
  728.          (SETQ M (COND ((ATOM M) M)
  729.                ((EQ (CAAR M) 'MPLUS)
  730.                 (CONS '(MPLUS)
  731.                   (MAPCAR #'$RATEXPAND (CDR M))))
  732.                (T ($RATEXPAND M))))
  733.          (RETURN (COND ((ATOM M) M)
  734.                ((EQ (CAAR M) 'MPLUS)
  735.                 (CONS '(MPLUS)
  736.                   (MAPCAR #'SSSQFR (CDR M))))
  737.                (T (SSSQFR M))))))
  738.  
  739. (DEFUN SSSQFR (X) (LET ((DOSIMP T)) (SIMPLIFY ($SQFR X))))
  740.  
  741. (DEFUN PFET1 (M V) 
  742.        (COND ((ATOM M) M)
  743.          ((EQ (CAAR M) 'MPLUS)
  744.           (CONS '(MPLUS)
  745.             (MAPCAR #'(LAMBDA (S) ($PARTFRAC S V))
  746.                 (CDR M))))
  747.          (T ($PARTFRAC M V))))
  748.  
  749. (DECLARE-TOP(UNSPECIAL M V))
  750.