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 / rat3a.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  25.9 KB  |  822 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 1980 Massachusetts Institute of Technology         ;;;
  9. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10.  
  11. (in-package "MAXIMA")
  12. (macsyma-module rat3a)
  13.  
  14. ;; This is the new Rational Function Package Part 1.
  15. ;; It includes most of the coefficient and polynomial routines required
  16. ;; by the rest of the functions.  Certain specialized functions are found
  17. ;; elsewhere.  Most low-level utility programs are also in this file.
  18.  
  19. (DECLARE-TOP(SPECIAL U* *A* *VAR* *X* V*)
  20.      (GENPREFIX A_1))
  21.  
  22. (declare-top (unspecial coef var exp p y ))  
  23. ;;These do not need to be special for this file and they
  24. ;;slow it down on lispm. We also eliminated the special
  25. ;;from ptimes2--wfs
  26.  
  27. (LOAD-MACSYMA-MACROS RATMAC)
  28.  
  29. ;; Global variables referenced throughout the rational function package.
  30.  
  31. (DEFMVAR MODULUS NIL "Global switch for doing modular arithmetic")
  32. (DEFMVAR HMODULUS NIL "Half of MODULUS")
  33. (DEFMVAR ERRRJFFLAG NIL "Controls action of ERRRJF (MAXIMA-ERROR or throw)")
  34.  
  35. ;;(DEFsubst POINTERGP (A B) (> (VALGET A) (VALGET B))) ;;better a macro in ratmac.
  36.  
  37. (DEFMACRO BCTIMES (&REST L)
  38.       `(REMAINDER (TIMES . ,L) MODULUS))
  39.  
  40. (DEFUN CQUOTIENT (A B)
  41.        (COND ((EQN A 0) 0)
  42.          ((NULL MODULUS)
  43.           (COND ((EQN 0 (CREMAINDER A B)) (QUOTIENT A B))
  44.             (T (ERRRJF "quotient is not exact"))))
  45.          (T (CTIMES A (CRECIP B)))))
  46.  
  47. (DEFUN ALG (L) (AND $ALGEBRAIC (NOT (ATOM L)) (GET (CAR L) 'TELLRAT)))
  48.  
  49. (DEFUN PACOEFP (X) (OR (PCOEFP X) (ALG X)))
  50.  
  51. (DEFUN LEADTERM (POLY)
  52.        (COND ((PCOEFP POLY) POLY)
  53.          (T (MAKE-POLY (P-VAR POLY) (P-LE POLY) (P-LC POLY)))))
  54.  
  55. ;; (DEFUN LCF (POLY)
  56. ;;        (COND ((PCOEFP POLY) POLY)
  57. ;;          (T (P-LC POLY))))
  58.  
  59. ;; (DEFUN LEXP (POLY) (COND ((PCOEFP POLY) 0) (T (P-LE POLY))))
  60.  
  61. (DEFUN CREMAINDER (A B)
  62.     (COND ((OR MODULUS (FLOATP A) (FLOATP B)) 0) 
  63.           ((REMAINDER A B))))
  64.  
  65. ;(DEFUN CBEXPT (P N)
  66. ;  (DO ((N (ASH N -1) (ASH N -1)) (S (IF (LOGTEST N 1) P 1)))
  67. ;      ((ZEROP N) S)
  68. ;    (SETQ P (BCTIMES P P))
  69. ;    (AND (ODDP N) (SETQ S (BCTIMES S P)))))
  70.  
  71. (DEFUN CBEXPT (P N)
  72.        (DO ((N (QUOTIENT N 2) (QUOTIENT N 2))
  73.         (S (COND ((ODDP N) P) (T 1))))
  74.        ((ZEROP N) S)
  75.        (SETQ P (BCTIMES P P))
  76.        (AND (ODDP N) (SETQ S (BCTIMES S P)))))
  77.  
  78.  
  79. #+PDP10
  80. (DEFUN CDIFFERENCE (X Y) (CPLUS X (CMINUS Y)))
  81.  
  82. ;; Coefficient Arithmetic -- coefficients are assumed to be something
  83. ;; that is NUMBERP in lisp.  If MODULUS is non-NIL, then all coefficients
  84. ;; are assumed to be less than its value.  Some functions use LOGAND
  85. ;; when MODULUS = 2.  This will not work with bignum coefficients.
  86.  
  87. ;; The following five routines have been hand coded in RAT;RATLAP >.  Please
  88. ;; do not delete these definitions.  They are needed for Macsymas which run
  89. ;; on non-PDP10 systems
  90.  
  91. #-(or PDP10 cl)
  92. (PROGN 'COMPILE
  93.  
  94. (DEFUN CMOD (N)
  95.   (COND ((NULL MODULUS) N)
  96.     ((BIGP MODULUS)
  97.      (SETQ N (REMAINDER N MODULUS))
  98.      (COND ((LESSP N 0)
  99.         (IF (LESSP (TIMES 2 N) (MINUS MODULUS))
  100.             (SETQ N (PLUS N MODULUS))))
  101.            ((GREATERP (TIMES 2 N) MODULUS) (SETQ N (DIFFERENCE N MODULUS)))) 
  102.      N)
  103.     ((BIGP N)
  104.      (SETQ N (REMAINDER N MODULUS))
  105.      (COND ((f< N 0)
  106.         (IF (< N (f- (LSH MODULUS -1))) (SETQ N (f+ N MODULUS))))
  107.            ((> N (LSH MODULUS -1)) (SETQ N (f- N MODULUS))))
  108.      N)
  109.     ((= MODULUS 2) (LOGAND N 1))
  110.     (T (LET ((NN (fixnum-remainder N MODULUS)))
  111.          (DECLARE (FIXNUM NN))
  112.          (COND ((f< NN 0)
  113.             (AND (f< NN (f- (ASH MODULUS -1))) (SETQ NN (f+ NN MODULUS))))
  114.            ((f> NN (ASH MODULUS -1)) (SETQ NN (f- NN MODULUS))))
  115.          NN))))
  116.  
  117. (DEFMFUN CPLUS (A B)
  118.   (COND ((NULL MODULUS) (PLUS A B))
  119.     ((BIGP MODULUS) (CMOD (PLUS A B)))
  120.     ((= MODULUS 2) (LOGAND 1 (LOGXOR A B)))
  121.     (T (LET ((NN (fixnum-remainder (f+ A B) MODULUS)))
  122.          (DECLARE (FIXNUM NN))
  123.          ;; This is an efficiency hack which assumes that MINUSP is
  124.          ;; faster than <.  The first clause could simply be
  125.          ;; ((< NN (f- (LSH MODULUS -1))) (f+ NN MODULUS))
  126.          (COND ((MINUSP NN)
  127.             (IF (f< NN (f- (ASH MODULUS -1))) (f+ NN MODULUS) NN))
  128.            ((f> NN (ASH MODULUS -1)) (f- NN MODULUS))
  129.            (T NN))))))
  130.  
  131. (DEFUN CDIFFERENCE (A B)
  132.   (COND ((NULL MODULUS) (DIFFERENCE A B))
  133.     ((BIGP MODULUS) (CMOD (DIFFERENCE A B)))
  134.     ((= MODULUS 2) (LOGAND 1 (LOGXOR A B)))
  135.  
  136.     (T (LET ((NN (fixnum-remainder (f- A B) MODULUS)))
  137.         ;;Note: the above code is invalid since the difference
  138.         ;;of two fixnum's may well be a bignum.
  139.         ;;now if they were already in the -Modulus/2 <=x < Modulus/2
  140.         ;;range then this would be ok.
  141.          (DECLARE (FIXNUM NN))
  142.          ;; This is an efficiency hack which assumes that MINUSP is
  143.          ;; faster than <.  The first clause could simply be
  144.          ;; ((< NN (f- (ASH MODULUS -1))) (f+ NN MODULUS))
  145.          (COND ((MINUSP NN)
  146.             (IF (f< NN (f- (ASH MODULUS -1))) (f+ NN MODULUS) NN))
  147.            ((f> NN (ASH MODULUS -1)) (f- NN MODULUS))
  148.            (T NN))))))
  149.  
  150. ;; (*DBLREM A B M) computes (A * B)/M in one swell foop so as not to
  151. ;; cons up a bignum intermediate for A*B.  MODULUS is known to be a fixnum
  152. ;; at the point where the call occurs, and thus A and B are known to be
  153. ;; fixnums.
  154.  
  155. #-cl
  156. (DECLARE-TOP(FIXNUM (*DBLREM FIXNUM FIXNUM FIXNUM)))
  157.  
  158. #+NIL
  159. (DEFUN *DBLREM (A B M)
  160.   ;(REMAINDER (TIMES A B) M)
  161.   ;(si:%dblrem a b c d) is (\ (f+ (f* a b) c) d).
  162.   ;This should probably be open-compiled...  Later, when we know it works.
  163.   ; Anyway, there is no interpreted si:%dblrem function, only the compiler's
  164.   ; hack.
  165.   (si:%dblrem a b 0 m))
  166.  
  167. #+CL
  168. (DEFUN *DBLREM (A B M)
  169.   ;; remind me to put the required primitives into NIL. -GJC
  170.   (REMAINDER (TIMES A B) M))
  171.  
  172. #+NIL
  173. (DEFUN *DBLREM (A B M)
  174.   ;; This generates the EMUL and EDIV instructions.
  175.   (COMPILER:%DBLREM A B 0 M))
  176.  
  177. (DEFMFUN CTIMES (A B)
  178.   (COND ((NULL MODULUS) (TIMES A B))
  179.     ((BIGP MODULUS) (CMOD (TIMES A B)))
  180.     ((= MODULUS 2) (LOGAND 1 A B))
  181.     (T (LET ((NN (*DBLREM A B MODULUS)))
  182.          (DECLARE (FIXNUM NN))
  183.          (COND ((MINUSP NN)
  184.             (IF (f< NN (f- (ASH MODULUS -1))) (f+ NN MODULUS) NN))
  185.            ((f> NN (ASH MODULUS -1)) (f- NN MODULUS))
  186.            (T NN)))))))
  187.  
  188. ;; Takes the inverse of an integer N mod P.  Solve N*X + P*Y = 1 
  189. ;; I suspect that N is guaranteed to be less than P, since in the case
  190. ;; where P is a fixnum, N is also assumed to be one.
  191.  
  192. (DEFUN INVMOD (N MODULUS) (CRECIP N))
  193.  
  194. #-lispm
  195. (DEFMFUN CRECIP (N)
  196.   (COND
  197.     ;; Have to use bignum arithmetic if modulus is a bignum
  198.     ((BIGP MODULUS) 
  199.      (PROG (A1 A2 Y1 Y2 Q (big-n N))
  200.        (IF (MINUSP Big-n) (SETQ Big-n (PLUS Big-n MODULUS)))
  201.        (SETQ A1 MODULUS A2 Big-n)
  202.        (SETQ Y1 0 Y2 1)
  203.        (GO STEP3)
  204.     STEP2 (SETQ Q (QUOTIENT A1 A2))
  205.        (PSETQ A1 A2 A2 (DIFFERENCE A1 (TIMES A2 Q)))
  206.        (PSETQ Y1 Y2 Y2 (DIFFERENCE Y1 (TIMES Y2 Q)))
  207.     STEP3 (COND ((ZEROP A2) (MERROR "Inverse of zero divisor?"))
  208.             ((NOT (EQUAL A2 1)) (GO STEP2)))
  209.        (RETURN (CMOD Y2))))
  210.     ;; Here we can use fixnum arithmetic
  211.     (T (PROG ((A1 0) (A2 0) (Y1 0) (Y2 0) (Q 0) (NN 0) (PP 0))
  212.          (DECLARE (FIXNUM A1 A2 Y1 Y2 Q NN PP))
  213.          (SETQ NN N PP MODULUS)
  214.          (COND ((MINUSP NN) (SETQ NN (f+ NN PP))))
  215.          (SETQ A1 PP A2 NN)
  216.          (SETQ Y1 0 Y2 1)
  217.          (GO STEP3)
  218.       STEP2 (SETQ Q (// A1 A2))
  219.          (PSETQ A1 A2 A2 (fixnum-remainder A1 A2))
  220.          (PSETQ Y1 Y2 Y2 (f- Y1 (f* Y2 Q)))
  221.       STEP3 (COND ((f= A2 0) (MERROR "Inverse of zero divisor?"))
  222.               ((NOT (f= A2 1)) (GO STEP2)))
  223.          ;; Is there any reason why this can't be (RETURN (CMOD Y2)) ? -cwh
  224.          (RETURN  (cmod y2)
  225. ;                    (COND ((= PP 2) (LOGAND 1 Y2))
  226. ;               (T (LET ((NN (fixnum-remainder Y2 PP)))
  227. ;                (DECLARE (FIXNUM NN))
  228. ;                (COND ((MINUSP NN)
  229. ;                       (AND (f< NN (f- (ASH PP -1)))
  230. ;                        (SETQ NN (f+ NN PP))))
  231. ;                      ((f> NN (ASH PP -1))
  232. ;                       (SETQ NN (f- NN PP))))
  233. ;                NN)))
  234.               )
  235.          ))))
  236.  
  237. (DEFUN CEXPT (N E)
  238.   (COND    ((NULL MODULUS) (EXPT N E))
  239.     (T (CMOD (CBEXPT N E)))))
  240.  
  241. ;; End of non-PDP10 definitions of CMOD, CPLUS, CDIFFERENCE, CTIMES, CRECIP, CEXPT.
  242.  
  243.  
  244.  
  245. ;; Notice how much cleaner things are on the Lisp Machine.
  246. ;; @@@@@ This definition may be better for NIL also, as long as the
  247. ;; arithmetic is hacked to not use the shadowed maclisp-fixnum functions,
  248. ;; = changes to eql, and (< n 0) -> (minusp n).
  249. #+lispm
  250. (DEFUN CRECIP (N &AUX (A1 MODULUS) A2 (Y1 0) (Y2 1) Q)
  251.   (SETQ A2 (IF (< N 0) (+ N MODULUS) N))
  252.   (DO ()
  253.       ((= A2 1) (CMOD Y2))
  254.       (IF (= A2 0)
  255.       (ERROR "Inverse of zero divisor -- ~D modulo ~D" N MODULUS))
  256.     (SETQ Q (// A1 A2))
  257.     (PSETQ A1 A2 A2 (- A1 (* A2 Q)))
  258.     (PSETQ Y1 Y2 Y2 (- Y1 (* Y2 Q)))))
  259.  
  260. ;;the following definitions are ok for 3600 and more transparent
  261. ;;and quicker.  Note for kcl, we provide some c definitions.
  262.  
  263. #+cl
  264. #-kcl
  265. (progn
  266. (defmacro mcmod (n) ;;; gives answers from -modulus/2 ...0 1 2 +modulus/2
  267.   `(let ((.n. (mod ,n modulus)))
  268.      (cond ((<= (times 2 .n.) modulus) .n.)
  269.        (t 
  270.         (difference .n. modulus)))))
  271.  
  272.  
  273. (defun cplus (a b)
  274.   (cond ((null modulus)(plus a b))
  275.     (t (mcmod (plus a b)))))
  276.  
  277. (defun cmod (n)
  278.   (cond ((null modulus ) n)
  279.         (t (mcmod n))))
  280.  
  281. (defun ctimes (a b)
  282.   (cond ((null modulus) (times a b))
  283.      (t (mcmod (times a b)))))
  284.  
  285. (defun cdifference (a b)
  286.   (cond ((null modulus) (- a b))
  287.      (t (mcmod (- a b)))))
  288.  
  289. )
  290.  
  291.  
  292.  
  293.  
  294. (DEFUN SETQMODULUS (M)
  295.    (COND ((NUMBERP M)
  296.       (COND ((GREATERP M 0)
  297.          (SETQ HMODULUS (QUOTIENT M 2))
  298.          (SETQ MODULUS M))
  299.         (T (MERROR "Modulus must be a number > 0"))))
  300.      (T (SETQ HMODULUS (SETQ MODULUS NIL)))))
  301.  
  302. (DEFMFUN PCOEFADD (E C X)
  303.   (COND ((PZEROP C) X)
  304.     ((BIGP E) (MERROR "Exponent out of range"))
  305.     (T (CONS E (CONS C X)))))
  306.  
  307. (DEFMFUN PPLUS (X Y)
  308.   (COND ((PCOEFP X) (PCPLUS X Y))
  309.     ((PCOEFP Y) (PCPLUS Y X))
  310.     ((EQ (P-VAR X) (P-VAR Y))
  311.      (PSIMP (P-VAR X) (PPLUS1 (P-TERMS Y) (P-TERMS X))))
  312.     ((POINTERGP (P-VAR X) (P-VAR Y))
  313.      (PSIMP (P-VAR X) (PCPLUS1 Y (P-TERMS X))))
  314.     (T (PSIMP (P-VAR Y) (PCPLUS1 X (P-TERMS Y))))))
  315.  
  316. (DEFUN PPLUS1 (X Y)
  317.        (COND ((PTZEROP X) Y)
  318.          ((PTZEROP Y) X)
  319.          ((= (PT-LE X) (PT-LE Y))
  320.           (PCOEFADD (PT-LE X)
  321.             (PPLUS (PT-LC X) (PT-LC Y))
  322.             (PPLUS1 (PT-RED X) (PT-RED Y))))
  323.          ((f> (PT-LE X) (PT-LE Y))
  324.           (CONS (PT-LE X) (CONS (PT-LC X) (PPLUS1 (PT-RED X) Y))))
  325.          (T (CONS (PT-LE Y) (CONS (PT-LC Y) (PPLUS1 X (PT-RED Y))))))) 
  326.  
  327. (DEFUN PCPLUS (C P) (COND ((PCOEFP P) (CPLUS P C))
  328.               (T (PSIMP (P-VAR P) (PCPLUS1 C (P-TERMS P))))))
  329.  
  330. (DEFUN PCPLUS1 (C X)
  331.        (COND ((NULL X)
  332.           (COND ((PZEROP C) NIL) (T (CONS 0 (CONS C NIL)))))
  333.          ((PZEROP (CAR X)) (PCOEFADD 0 (PPLUS C (CADR X)) NIL))
  334.          (T (CONS (CAR X) (CONS (CADR X) (PCPLUS1 C (CDDR X)))))))
  335.      
  336.  
  337. (DEFMFUN PDIFFERENCE (X Y)
  338.   (COND ((PCOEFP X) (PCDIFFER X Y))
  339.     ((PCOEFP Y) (PCPLUS (CMINUS Y) X))
  340.     ((EQ (P-VAR X) (P-VAR Y))
  341.      (PSIMP (P-VAR X) (PDIFFER1 (P-TERMS X) (P-TERMS Y))))
  342.     ((POINTERGP (P-VAR X) (P-VAR Y))
  343.      (PSIMP (P-VAR X) (PCDIFFER2 (P-TERMS X) Y)))
  344.     (T (PSIMP (P-VAR Y) (PCDIFFER1 X (P-TERMS Y))))))
  345.  
  346. (DEFUN PDIFFER1 (X Y)
  347.        (COND ((PTZEROP X) (PMINUS1 Y))
  348.          ((PTZEROP Y) X)
  349.          ((= (PT-LE X) (PT-LE Y))
  350.           (PCOEFADD (PT-LE X)
  351.             (PDIFFERENCE (PT-LC X) (PT-LC Y))
  352.             (PDIFFER1 (PT-RED X) (PT-RED Y))))
  353.          ((f> (PT-LE X) (PT-LE Y))
  354.           (CONS (PT-LE X) (CONS (PT-LC X) (PDIFFER1 (PT-RED X) Y))))
  355.          (T (CONS (PT-LE Y) (CONS (PMINUS (PT-LC Y))
  356.                       (PDIFFER1 X (PT-RED Y)))))))
  357.  
  358. (DEFUN PCDIFFER (C P)
  359.        (COND ((PCOEFP P) (CDIFFERENCE C P))
  360.          (T (PSIMP (CAR P) (PCDIFFER1 C (CDR P))))))     
  361.  
  362. (DEFUN PCDIFFER1 (C X)
  363.        (COND ((NULL X) (COND ((PZEROP C) NIL) (T (LIST 0 C))))
  364.          ((PZEROP (CAR X))
  365.           (PCOEFADD 0 (PDIFFERENCE C (CADR X)) NIL))
  366.          (T (CONS (CAR X)
  367.               (CONS (PMINUS (CADR X)) (PCDIFFER1 C (CDDR X)))))))
  368.  
  369. (DEFUN PCDIFFER2 (X C)
  370.        (COND ((NULL X) (COND ((PZEROP C) NIL) (T (LIST 0 (PMINUS C)))))
  371.          ((PZEROP (CAR X)) (PCOEFADD 0 (PDIFFERENCE (CADR X) C) NIL))
  372.          (T (CONS (CAR X) (CONS (CADR X) (PCDIFFER2 (CDDR X) C))))))
  373.  
  374. (DEFUN PCSUBSTY (VALS VARS P)                ;list of vals for vars
  375.        (COND ((NULL VARS) P)
  376.          ((ATOM VARS) (PCSUB P (LIST VALS) (LIST VARS)))    ;one val hack
  377.          (T (SETQ VARS (SORTCAR (MAPCAR #'CONS VARS VALS) #'POINTERGP))
  378.         (PCSUB P (MAPCAR (FUNCTION CDR) VARS)
  379.                  (MAPCAR (FUNCTION CAR) VARS)))))
  380.  
  381. (DEFUN PCSUBST (P VAL VAR)                ;one val for one var
  382.   (COND ((PCOEFP P) P)
  383.     ((EQ (CAR P) VAR) (PCSUB1 (CDR P) VAL () ()))
  384.     ((POINTERGP VAR (CAR P)) P)
  385.     (T (PSIMP (CAR P) (PCSUB2 (CDR P) (NCONS VAL) (NCONS VAR)))))) 
  386.  
  387. (DEFUN PCSUB1 (A VAL VALS VARS)
  388.        (if (EQUAL VAL 0) (PCSUB (PTERM A 0) VALS VARS)
  389.        (do ((p (pt-red a) (pt-red p))
  390.         (ans (pcsub (pt-lc a) vals vars)
  391.              (pplus (ptimes ans
  392.                     (pexpt val (f- ld (pt-le p))))
  393.                 (pcsub (pt-lc p) vals vars)))
  394.         (ld (pt-le a) (pt-le p)))
  395.            ((null p) (ptimes ans (pexpt val ld))))))
  396.  
  397. (DEFUN PCSUB (P VALS VARS)
  398.        (COND ((NULL VALS) P)
  399.          ((PCOEFP P) P)
  400.          ((EQ (p-var P) (CAR VARS)) (PCSUB1 (p-terms P) (car vals)
  401.                           (cdr VALS) (cdr VARS)))
  402.          ((POINTERGP (CAR VARS) (p-var P))
  403.           (PCSUB P (CDR VALS) (CDR VARS)))
  404.          (T (PSIMP (p-var P) (PCSUB2 (p-terms P) VALS VARS)))))
  405.  
  406. (DEFUN PCSUB2 (terms VALS VARS)
  407.   (sloop for (exp coef) on terms by 'pt-red
  408.     unless (pzerop (setq coef (pcsub coef vals vars)))
  409.     nconc (list exp coef)))
  410.  
  411.  
  412.  
  413. ;    THIS DERIVATIVE PROGRAM ASSUMES NO DEPENDENCIES HIDDEN ON VARLIST
  414. ;    OR ELSWHERE.  COMPARE TO RATDX.
  415.  
  416. ;(DEFMFUN PDERIVATIVE (P *VAR*)
  417. ;  (IF (PCOEFP P) (CDERIVATIVE P *VAR*)
  418. ;      (PSIMP (P-VAR P) (COND ((EQ *VAR* (P-VAR P))
  419. ;                  (PDERIVATIVE2 (P-TERMS P)))
  420. ;                 ((POINTERGP *VAR* (P-VAR P))
  421. ;                  (PTZERO))
  422. ;                 (T (PDERIVATIVE3 (P-TERMS P)))))))
  423. ;
  424. ;(DEFUN PDERIVATIVE2 (X)
  425. ;  (COND ((NULL X) NIL)
  426. ;    ((ZEROP (PT-LE X)) NIL)
  427. ;    (T (PCOEFADD (f1- (PT-LE X))
  428. ;             (PCTIMES (CMOD (PT-LE X)) (PT-LC X))
  429. ;             (PDERIVATIVE2 (PT-RED X))))))
  430. ;
  431. ;(DEFUN PDERIVATIVE3 (X)
  432. ;       (COND ((NULL X) NIL)
  433. ;         (T (PCOEFADD
  434. ;         (PT-LE X)
  435. ;         (PDERIVATIVE (PT-LC X) *VAR*)
  436. ;         (PDERIVATIVE3 (PT-RED X))))))
  437.  
  438. ;;replaced above by version without special bindings.
  439. (DEFMFUN PDERIVATIVE (P VARI)
  440.   (IF (PCOEFP P) (CDERIVATIVE P VARI)
  441.       (PSIMP (P-VAR P) (COND ((EQ VARI (P-VAR P))
  442.                   (PDERIVATIVE2 (P-TERMS P)))
  443.                  ((POINTERGP VARI (P-VAR P))
  444.                   (PTZERO))
  445.                  (T (PDERIVATIVE3 (P-TERMS P) VARI))))))
  446.  
  447. (DEFUN PDERIVATIVE2 (X)
  448.   (COND ((NULL X) NIL)
  449.     ((ZEROP (PT-LE X)) NIL)
  450.     (T (PCOEFADD (f1- (PT-LE X))
  451.              (PCTIMES (CMOD (PT-LE X)) (PT-LC X))
  452.              (PDERIVATIVE2 (PT-RED X))))))
  453.  
  454. (DEFUN PDERIVATIVE3 (X VARI)
  455.        (COND ((NULL X) NIL)
  456.          (T (PCOEFADD
  457.          (PT-LE X)
  458.          (PDERIVATIVE (PT-LC X) VARI)
  459.          (PDERIVATIVE3 (PT-RED X) VARI)))))
  460.  
  461. (DEFMFUN PDIVIDE (X Y)
  462.        (COND ((PZEROP Y) (ERRRJF "Quotient by zero"))
  463.          ((PACOEFP Y) (LIST (RATREDUCE X Y) (RZERO)))
  464.          ((PACOEFP X) (LIST (RZERO) (CONS X 1)))
  465.          ((POINTERGP (CAR X) (CAR Y)) (LIST (RATREDUCE X Y) (RZERO)))
  466.          (T (PDIVIDE1 X Y))))
  467.  
  468. (DEFUN PDIVIDE1 (U V)
  469.   (PROG (K INC LCU LCV Q R)
  470.     (SETQ LCV (CONS (CADDR V) 1))
  471.     (SETQ Q (RZERO))
  472.     (SETQ R (CONS U 1) )
  473.    A    (SETQ K (f- (PDEGREE (CAR R) (P-VAR V)) (P-LE V)))
  474.         (IF (MINUSP K) (RETURN (LIST Q R)))
  475.     (SETQ LCU (CONS (P-LC (CAR R)) (CDR R)))
  476.     (SETQ INC (RATQUOTIENT LCU LCV))
  477.     (SETQ INC (CONS (PSIMP (CAR V) (LIST K (CAR INC)))
  478.             (CDR INC)))
  479.     (SETQ Q (RATPLUS Q INC))
  480.     (SETQ R (RATPLUS R  (RATTIMES (CONS (PMINUS V) 1) INC T)))
  481.     (GO A)))
  482.      
  483.  
  484. (DEFMFUN PEXPT (P N)
  485.   (COND ((= N 0) 1)
  486.     ((= N 1) P)
  487.     ((MINUSP N) (PQUOTIENT 1 (PEXPT P (f- N))))
  488.     ((PCOEFP P) (CEXPT P N))
  489.     ((ALG P) (PEXPTSQ P N))
  490.     ((NULL (P-RED P))
  491.      (PSIMP (P-VAR P)
  492.         (PCOEFADD (f* N (P-LE P)) (PEXPT (P-LC P) N) NIL)))
  493.     (T (LET ((*A* (f1+ N))
  494.          (*X* (PSIMP (P-VAR P) (P-RED P)))
  495.          (B (MAKE-POLY (P-VAR P) (P-LE P) (P-LC P))))
  496.          (DO ((BL (LIST (PEXPT B 2) B)
  497.               (CONS (PTIMES B (CAR BL)) BL))
  498.           (M 2 (f1+ M)))
  499.          ((= M N) (PPLUS (CAR BL)
  500.                  (PEXPT2 1 1 *X* (CDR BL)))))))))
  501.  
  502. (DEFUN NXTBINCOEF (M NOM) (QUOTIENT (TIMES NOM (DIFFERENCE *A* M)) M))
  503.  
  504. (DEFUN PEXPT2 (M NOM B L)
  505.   (IF (NULL L) B
  506.       (PPLUS (PTIMES (PCTIMES (CMOD (SETQ NOM (NXTBINCOEF M NOM))) B) (CAR L))
  507.          (PEXPT2 (f1+ M)
  508.              NOM
  509.              (IF (EQ *X* B) (PEXPT B 2)
  510.              (PTIMES *X* B))
  511.              (CDR L)))))
  512.  
  513. (DEFMFUN PMINUSP (P) (IF (NUMBERP P) (MINUSP P)
  514.              (PMINUSP (P-LC P))))
  515.  
  516. (DEFMFUN PMINUS (P) (IF (PCOEFP P) (CMINUS P)
  517.             (CONS (P-VAR P) (PMINUS1 (P-TERMS P)))))
  518.  
  519. (DEFUN PMINUS1 (X)
  520.   (sloop for (exp coef) on x by 'pt-red
  521.     nconc (list exp (pminus coef))))
  522.  
  523. (DEFMFUN PMOD (P)
  524.   (IF (PCOEFP P) (CMOD P)
  525.       (PSIMP (CAR P)
  526.          (sloop for (exp coef) on (p-terms p) by 'pt-red
  527.            unless (pzerop (setq coef (pmod coef)))
  528.            nconc (list exp coef)))))
  529.  
  530. (DEFMFUN PQUOTIENT (X Y)
  531.  (COND ((PCOEFP X)
  532.     (COND ((PZEROP X) (PZERO))
  533.           ((PCOEFP Y) (CQUOTIENT X Y))
  534.           ((ALG Y) (PAQUO X Y))
  535.           (T (ERRRJF "Quotient by a polynomial of higher degree"))))
  536.        ((PCOEFP Y) (COND ((PZEROP Y) (ERRRJF "Quotient by zero"))
  537.              (MODULUS (PCTIMES (CRECIP Y) X))
  538.              (T (PCQUOTIENT X Y))))
  539.        ((ALG Y) (OR (LET ((ERRRJFFLAG T) $ALGEBRAIC)
  540.               (CATCH 'RATERR (PQUOTIENT X Y)))
  541.              (PATIMES X (RAINV Y))))
  542.        ((POINTERGP (P-VAR X) (P-VAR Y)) (PCQUOTIENT X Y))
  543.        ((OR (POINTERGP (P-VAR Y) (P-VAR X)) (f> (P-LE Y) (P-LE X)))
  544.     (ERRRJF "Quotient by a polynomial of higher degree"))
  545.        (T (PSIMP (P-VAR X) (PQUOTIENT1 (P-TERMS X) (P-TERMS Y))))))
  546.  
  547. (DEFUN PCQUOTIENT (P Q)
  548.   (psimp (P-VAR p) (pcquotient1 (p-terms p) q)))
  549.  
  550. (DEFUN PCQUOTIENT1 (P1 Q)
  551.   (sloop for (exp coef) on P1 by 'PT-RED
  552.     nconc (list exp (pquotient coef Q))))
  553.  
  554. (DECLARE-TOP(SPECIAL K Q*)
  555.      (FIXNUM K I))
  556.  
  557. (DEFUN PQUOTIENT1 (U V &AUX Q* (k 0))
  558.   (declare (fixnum k))
  559.   (sloop do (setq  k (f- (pt-le u) (pt-le v)))
  560.     when (minusp k) do (ERRRJF "Polynomial quotient is not exact")
  561.     nconc (list k (setq q* (pquotient (pt-lc u) (pt-lc v))))
  562.     until (ptzerop (setq u (pquotient2 (pt-red u) (pt-red v))))))
  563.  
  564. (DEFUN PQUOTIENT2 (X Y &AUX (I 0))            ;X-v^k*Y*Q*
  565.   (COND ((NULL Y) X)
  566.     ((NULL X) (PCETIMES1 Y K (PMINUS Q*)))
  567.     ((MINUSP (SETQ I (f- (PT-LE X) (PT-LE Y) K)))
  568.      (PCOEFADD (f+ (PT-LE Y) K)
  569.            (PTIMES Q* (PMINUS (PT-LC Y)))
  570.            (PQUOTIENT2 X (PT-RED Y))))
  571.     ((ZEROP I) (PCOEFADD (PT-LE X)
  572.                  (PDIFFERENCE (PT-LC X) (PTIMES Q* (PT-LC Y)))
  573.                  (PQUOTIENT2 (PT-RED X) (PT-RED Y))))
  574.     (T (CONS (PT-LE X) (CONS (PT-LC X) (PQUOTIENT2 (PT-RED X) Y))))))
  575.  
  576. (DECLARE-TOP(UNSPECIAL K Q*)
  577.      (NOTYPE K I))
  578.  
  579. (defun algord (var) (and $algebraic (get var 'algord)))
  580.  
  581. (DEFUN PSIMP (VAR X)
  582.        (COND ((PTZEROP X) 0)
  583.          ((ATOM X) X)
  584.          ((ZEROP (PT-LE X)) (PT-LC X))
  585.          ((algord var)                ;wrong alg ordering
  586.           (do ((p x (cddr p)) (sum 0))
  587.           ((null p) (cond ((pzerop sum) (cons var x))
  588.                   (t (pplus sum (psimp2 var x)))))
  589.         (unless (or (pcoefp (cadr p)) (pointergp var (caadr p)))
  590.             (setq sum (pplus sum
  591.                      (if (zerop (pt-le p)) (pt-lc p)
  592.                          (ptimes 
  593.                           (make-poly var (pt-le p) 1)
  594.                           (pt-lc p)))))
  595.             (setf (pt-lc p) 0))))
  596.          (T (CONS VAR X))))
  597.  
  598. (defun psimp1 (var x) (let ($algebraic) (psimp var x)))
  599.  
  600. (defun psimp2 (var x)
  601.        (do ((p (setq x (cons nil x)) ))
  602.        ((null (cdr p)) (psimp1 var (cdr x)))
  603.        (cond ((pzerop (caddr p)) (rplacd p (cdddr p)))
  604.          (t (setq p (cddr p))))))
  605.  
  606. (DEFUN PTERM  (P N)
  607.        (do ((p p (pt-red p)))
  608.        ((ptzerop p) (pzero))
  609.        (cond ((f< (pt-le p) n) (return (pzero)))
  610.          ((f= (pt-le p) n) (return (pt-lc p))))))
  611.  
  612.  
  613. (DEFMFUN PTIMES (X Y)
  614.   (COND ((PCOEFP X) (IF (PZEROP X) (PZERO) (PCTIMES X Y)))
  615.     ((PCOEFP Y) (IF (PZEROP Y) (PZERO) (PCTIMES Y X)))
  616.     ((EQ (P-VAR X) (P-VAR Y))
  617.      (PALGSIMP (P-VAR X) (PTIMES1 (P-TERMS X) (P-TERMS Y)) (ALG X)))
  618.     ((POINTERGP (P-VAR X) (P-VAR Y))
  619.      (PSIMP (P-VAR X) (PCTIMES1 Y (P-TERMS X))))
  620.     (T (PSIMP (P-VAR Y) (PCTIMES1 X (P-TERMS Y))))))
  621. ;
  622. ;(DEFUN PTIMES1 (X Y &AUX U*)
  623. ;  (DO ((V* (SETQ U* (PCETIMES1 Y (PT-LE X) (PT-LC X))))
  624. ;       (X (PT-RED X) (PT-RED X)))
  625. ;      ((PTZEROP X) U*)
  626. ;    (PTIMES2 Y (PT-LE X) (PT-LC X))))
  627. ;
  628. ;
  629. ;
  630. ;(DEFUN PTIMES2 (Y XE XC) 
  631. ;  (PROG (E U C) 
  632. ;     A1 (COND ((NULL Y) (RETURN NIL)))
  633. ;    (SETQ E (f+ XE (CAR Y)))
  634. ;    (SETQ C (PTIMES (CADR Y) XC))
  635. ;    (COND ((PZEROP C) (SETQ Y (CDDR Y)) (GO A1))
  636. ;          ((OR (NULL V*) (> E (CAR V*)))
  637. ;           (SETQ U* (SETQ V* (PPLUS1 U* (LIST E C))))
  638. ;           (SETQ Y (CDDR Y)) (GO A1))
  639. ;          ((= E (CAR V*))
  640. ;           (SETQ C (PPLUS C (CADR V*)))
  641. ;           (COND ((PZEROP C)
  642. ;              (SETQ U* (SETQ V* (PDIFFER1 U* (LIST (CAR V*) (CADR V*))))))
  643. ;             (T (RPLACA (CDR V*) C)))
  644. ;           (SETQ Y (CDDR Y))
  645. ;           (GO A1)))
  646. ;     A  (COND ((AND (CDDR V*) (> (CADDR V*) E)) (SETQ V* (CDDR V*)) (GO A)))
  647. ;    (SETQ U (CDR V*))
  648. ;     B  (COND ((OR (NULL (CDR U)) (< (CADR U) E))
  649. ;           (RPLACD U (CONS E (CONS C (CDR U)))) (GO E)))
  650. ;    (COND ((PZEROP (SETQ C (PPLUS (CADDR U) C)))
  651. ;           (RPLACD U (CDDDR U)) (GO D))
  652. ;          (T (RPLACA (CDDR U) C)))
  653. ;     E  (SETQ U (CDDR U))
  654. ;     D  (SETQ Y (CDDR Y))
  655. ;    (COND ((NULL Y) (RETURN NIL)))
  656. ;    (SETQ E (f+ XE (CAR Y)))
  657. ;    (SETQ C (PTIMES (CADR Y) XC))
  658. ;     C  (COND ((AND (CDR U) (> (CADR U) E)) (SETQ U (CDDR U)) (GO C)))
  659. ;    (GO B)))
  660.  
  661.  
  662. (DEFUN   PTIMES1 (X Y-ORIG &AUX UUU  )
  663.   (DO ((VVV (SETQ UUU (PCETIMES1 Y-ORIG (PT-LE X) (PT-LC X))))
  664.        (X (PT-RED X) (PT-RED X)))
  665.       ((PTZEROP X) UUU)
  666.     (let ((y Y-ORIG) (xe (PT-LE X)) (xc (PT-LC X)))
  667.   (PROG (E U C) 
  668.      A1 (COND ((NULL Y) (RETURN NIL)))
  669.     (SETQ E (f+ XE (CAR Y)))
  670.     (SETQ C (PTIMES (CADR Y) XC))
  671.     (COND ((PZEROP C) (SETQ Y (CDDR Y)) (GO A1))
  672.           ((OR (NULL VVV) (f> E (CAR VVV)))
  673.            (SETQ UUU (SETQ VVV (PPLUS1 UUU (LIST E C))))
  674.            (SETQ Y (CDDR Y)) (GO A1))
  675.           ((f= E (CAR VVV))
  676.            (SETQ C (PPLUS C (CADR VVV)))
  677.            (COND ((PZEROP C)
  678.               (SETQ UUU (SETQ VVV (PDIFFER1 UUU (LIST (CAR VVV) (CADR VVV))))))
  679.              (T (RPLACA (CDR VVV) C)))
  680.            (SETQ Y (CDDR Y))
  681.            (GO A1)))
  682.      A  
  683.     (COND ((AND (CDDR VVV) (f> (CADDR VVV) E))
  684.            (SETQ VVV (CDDR VVV)) (GO A)))
  685.     (SETQ U (CDR VVV ))
  686.      B  (COND ((OR (NULL (CDR U)) (f< (CADR U) E))
  687.            (RPLACD U (CONS E (CONS C (CDR U)))) (GO E)))
  688.     (COND ((PZEROP (SETQ C (PPLUS (CADDR U) C)))
  689.            (RPLACD U (CDDDR U)) (GO D))
  690.           (T (RPLACA (CDDR U) C)))
  691.      E  (SETQ U (CDDR U))
  692.      D  (SETQ Y (CDDR Y))
  693.     (COND ((NULL Y) (RETURN NIL)))
  694.     (SETQ E (f+ XE (CAR Y)))
  695.     (SETQ C (PTIMES (CADR Y) XC))
  696.      C  (COND ((AND (CDR U) (f> (CADR U) E)) (SETQ U (CDDR U)) (GO C)))
  697.     (GO B))))
  698.   uuu)
  699.  
  700. (DEFUN PCETIMES1 (Y E C)                ;C*V^E*Y
  701.   (SLOOP FOR (EXP COEF) ON Y BY 'PT-RED
  702.     UNLESS (PZEROP (SETQ COEF (PTIMES C COEF)))
  703.     NCONC (LIST (f+ E EXP) COEF)))
  704.  
  705. (DEFUN PCTIMES (C P)
  706.   (IF (PCOEFP P) (CTIMES C P)
  707.       (PSIMP (P-VAR P) (PCTIMES1 C (P-TERMS P)))))
  708.  
  709. (DEFUN PCTIMES1 (C terms)
  710.   (sloop for (exp coef) on terms by 'pt-red
  711.     unless (pzerop (setq coef (PTIMES C coef)))
  712.     nconc (list exp coef)))
  713.  
  714. (DEFUN LEADALGCOEF (P)
  715.        (COND ((PACOEFP P) P)
  716.          (T (LEADALGCOEF (P-LC P))) ))
  717.  
  718. (DEFUN PAINVMOD (Q)
  719.        (COND ((PCOEFP Q) (CRECIP Q))
  720.          (T (PAQUO (LIST (CAR Q) 0 1) Q ))))
  721.  
  722. (DEFUN PALGSIMP (VAR P TELL)    ;TELL=(N X) -> X^(1/N)
  723.        (PSIMP VAR (COND ((OR (NULL TELL) (NULL P)
  724.                  (f< (CAR P) (CAR TELL))) P)
  725.             ((NULL (CDDR TELL)) (PASIMP1 P (CAR TELL) (CADR TELL)))
  726.             (T (PGCD1 P TELL)) )))
  727.  
  728. (DEFUN PASIMP1 (P DEG KERNEL)                ;assumes deg>=(car p)
  729.        (DO ((A P (PT-RED A))
  730.         (B P A))
  731.        ((OR (NULL A) (< (PT-LE A) DEG))
  732.         (RPLACD (CDR B) NIL)
  733.         (PPLUS1 (PCTIMES1 KERNEL P) A))
  734.      (RPLACA A (f- (PT-LE A) DEG))))
  735.  
  736. (DEFMFUN MONIZE (P) 
  737.        (COND ((PCOEFP P) (IF (PZEROP P) P 1))
  738.          (T (CONS (P-VAR P) (PMONICIZE (COPY1 (P-TERMS P)))))))
  739.  
  740. (DEFUN PMONICIZE (P)                    ;CLOBBERS POLY
  741.        (COND ((EQN (PT-LC P) 1) P)
  742.          (T (PMON1 (PAINVMOD (LEADALGCOEF (PT-LC P))) P) P)))
  743.  
  744. (DEFUN PMON1 (MULT L)
  745.   (COND (L (PMON1 MULT (PT-RED L))
  746.        (SETF (PT-LC L) (PTIMES MULT (PT-LC L))))))
  747.  
  748. (DEFUN PMONZ (POLY &AUX LC)                ;A^(N-1)*P(X/A)
  749.   (SETQ POLY (PABS POLY))       
  750.   (COND ((EQUAL (SETQ LC (P-LC POLY)) 1) POLY)
  751.     (T (DO ((P (P-RED POLY) (PT-RED P))
  752.         (P1 (MAKE-POLY (P-VAR POLY) (P-LE POLY) 1))
  753.         (MULT 1)
  754.         (DEG (f1- (P-LE POLY)) (PT-LE P)))
  755.            ((NULL P) P1)
  756.          (SETQ MULT (PTIMES MULT (PEXPT LC (f- DEG (PT-LE P)))))
  757.          (NCONC P1 (LIST (PT-LE P) (PTIMES MULT (PT-LC P))))))))
  758.  
  759. ;    THESE ARE ROUTINES FOR MANIPULATING ALGEBRAIC NUMBERS
  760.  
  761. (DEFUN ALGNORMAL (P) (CAR (RQUOTIENT P (LEADALGCOEF P))))
  762.  
  763. (DEFUN ALGCONTENT (P)
  764.        (LET* ((LCF (LEADALGCOEF P))
  765.           ((PRIM . DENOM) (RQUOTIENT P LCF)))
  766.          (LIST (RATREDUCE LCF DENOM) PRIM)))
  767.  
  768. (DEFUN RQUOTIENT (P Q &aux ALGFAC* A E)  ;FINDS PSEUDO QUOTIENT IF PSEUDOREM=0
  769.    (COND ((EQUAL P Q) (CONS 1 1))
  770.      ((PCOEFP Q) (RATREDUCE P Q))
  771.      ((SETQ A (TESTDIVIDE P Q)) (CONS A 1))
  772.      ((ALG Q) (RATTIMES (CONS P 1) (RAINV Q) T))
  773.      (T (COND ((ALG (SETQ A (LEADALGCOEF Q)))
  774.            (SETQ A (RAINV A))
  775.            (SETQ P (PTIMES P (CAR A)))
  776.            (SETQ Q (PTIMES Q (CAR A)))
  777.            (SETQ A (CDR A)) ))
  778.         (COND ((MINUSP (SETQ E (f+ 1 (f- (CADR Q)) (PDEGREE P (CAR Q)))))
  779.            (ERRRJF "Quotient by a polynomial of higher degree")))
  780.         (SETQ A (PEXPT A E))
  781.         (RATREDUCE (OR (TESTDIVIDE (PTIMES A P) Q)
  782.                (PROG2 (SETQ A (PEXPT (P-LC Q) E))
  783.                   (PQUOTIENT (PTIMES A P) Q)))
  784.                A)) ))
  785.  
  786. (DEFUN PATIMES (X R) (PQUOTIENTCHK (PTIMES X (CAR R)) (CDR R)))
  787.  
  788. (DEFUN PAQUO (X Y) (PATIMES X (RAINV Y)))
  789.  
  790. (DEFUN MPGET (VAR)
  791.        (COND ((NULL (SETQ VAR (ALG VAR))) NIL)
  792.          ((CDDR VAR) VAR)
  793.          (T (LIST (CAR VAR) 1 0 (PMINUS (CADR VAR))))))
  794.  
  795.  
  796. (DEFUN RAINV (Q)
  797.   (COND ((PCOEFP Q)
  798.      (COND (MODULUS (CONS (CRECIP Q) 1))
  799.            (T (CONS 1 Q))))
  800.     (T (LET ((VAR (CAR Q)) (P (MPGET Q)))
  801.          (declare (special var))    ;who uses this? --gsb
  802.          (COND ((NULL P) (CONS 1 Q))
  803.            (T (SETQ P (CAR (LET ($RATALGDENOM)
  804.                      (BPROG Q (CONS VAR P)))))
  805.               (RATTIMES (CONS (CAR P) 1) (RAINV (CDR P)) T)))))))
  806.  
  807. (DEFUN PEXPTSQ (P N)
  808.   (DO ((N (QUOTIENT N 2) (QUOTIENT N 2))
  809.        (S (COND ((ODDP N) P) (T 1))))
  810.       ((ZEROP N) S)
  811.     (SETQ P (PTIMES P P))
  812.     (AND (ODDP N) (SETQ S (PTIMES S P))) ))
  813.  
  814. ;    THIS IS THE END OF THE NEW RATIONAL FUNCTION PACKAGE PART 1.
  815.  
  816. ;; Someone should determine which of the variables special in this file
  817. ;; (and others) are really special and which are used just for dynamic
  818. ;; scoping.  -- cwh
  819.  
  820. #-NIL ;Minimize complaints.
  821. (DECLARE-TOP(UNSPECIAL V* *A* U* *VAR*))
  822.