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 13(macsyma-module rat3a) 14 15;; This is the new Rational Function Package Part 1. 16;; It includes most of the coefficient and polynomial routines required 17;; by the rest of the functions. Certain specialized functions are found 18;; elsewhere. Most low-level utility programs are also in this file. 19 20(declare-top (unspecial coef var exp p y)) 21 22;;These do not need to be special for this file and they 23;;slow it down on lispm. We also eliminated the special 24;;from ptimes2--wfs 25 26;; Global variables referenced throughout the rational function package. 27 28(defmvar modulus nil "Global switch for doing modular arithmetic") 29 30;; CQUOTIENT 31;; 32;; Calculate the quotient of two coefficients, which should be numbers. If 33;; MODULUS is non-nil, we try to take the reciprocal of A with respect to the 34;; modulus (using CRECIP) and then multiply by B. Note that this fails if B 35;; divides A as an integer, but B is not a unit in the ring of integers modulo 36;; MODULUS. For example, 37;; 38;; (let ((modulus 20)) (cquotient 10 5)) => ERROR 39;; 40;; If MODULUS is nil, then we work over the ring of integers when A and B are 41;; integers, and raise a RAT-ERROR if A is not divisible by B. If either A or B 42;; is a float then the division is done in floating point. Floats can get as far 43;; as the rational function code if $KEEPFLOAT is true. 44(defun cquotient (a b) 45 (cond ((equal a 0) 0) 46 ((null modulus) 47 (if (or (floatp a) (floatp b) 48 (zerop (rem a b))) 49 (/ a b) 50 (rat-error "CQUOTIENT: quotient is not exact"))) 51 (t (ctimes a (crecip b))))) 52 53;; ALG 54;; 55;; Get any value stored on the tellrat property of (car l). Returns NIL if L 56;; turns out not to be a list or if $ALGEBRAIC is false. 57(defun alg (l) 58 (unless (atom l) (algv (car l)))) 59 60;; PACOEFP 61;; 62;; Return T if either X is a bare coefficient or X is a polynomial whose main 63;; variable has a declared value as an algebraic integer. Otherwise return NIL. 64(defun pacoefp (x) 65 (and (or (pcoefp x) (alg x)) 66 T)) 67 68;; LEADTERM 69;; 70;; Return the leading term of POLY as a polynomial itself. 71(defun leadterm (poly) 72 (if (pcoefp poly) 73 poly 74 (make-poly (p-var poly) (p-le poly) (p-lc poly)))) 75 76;; CRECIP 77;; 78;; Takes the inverse of an integer N mod MODULUS. If there is a modulus then the 79;; result is constrained to lie in (-modulus/2, modulus/2] 80;; 81;; This just uses the extended Euclidean algorithm: once you have found a,b such 82;; that a*n + b*modulus = 1 then a must be the reciprocal you're after. 83;; 84;; When MODULUS is greater than 2^15, we use exactly the same algorithm in 85;; CRECIP-GENERAL, but it can't use fixnum arithmetic. Note: There's no 86;; particular reason to target 32 bits except that trying to work out the right 87;; types on the fly looks complicated and this lisp compiler, at least, uses 32 88;; bit words. Since we have to take a product, we constrain the types to 16 bit 89;; numbers. 90(defun crecip (n) 91 ;; Punt on anything complicated 92 (unless (and modulus (typep modulus '(unsigned-byte 15))) 93 (return-from crecip (crecip-general n))) 94 95 ;; And make sure that -MODULUS < N < MODULUS 96 (unless (<= (- modulus) n modulus) 97 (merror "N is out of range [-MODULUS, MODULUS] in crecip.")) 98 99 ;; N in [0, MODULUS] 100 (when (minusp n) (setf n (+ n modulus))) 101 102 ;; The mod-copy parameter stores a copy of MODULUS on the stack, which is 103 ;; useful because the lisp implementation doesn't know that the special 104 ;; variable MODULUS is still an (unsigned-byte 15) when we get to the end 105 ;; (since it can't tell that our function calls don't change it behind our 106 ;; backs, I guess) 107 (let ((mod modulus) (remainder n) (a 1) (b 0) 108 (mod-copy modulus)) 109 ;; On SBCL in 2013 at least, the compiler doesn't spot that MOD and 110 ;; REMAINDER are unsigned and bounded above by MODULUS, a 16-bit integer. So 111 ;; we have to tell it. Also, the lisp implementation can't really be 112 ;; expected to know that Bezout coefficients are bounded by the modulus and 113 ;; remainder, so we have to tell it that too. 114 (declare (type (unsigned-byte 15) mod mod-copy remainder) 115 (type (signed-byte 16) a b)) 116 117 (loop 118 until (= remainder 1) 119 120 when (zerop remainder) do 121 (merror (intl:gettext "CRECIP: ~M does not have an inverse with modulus=~M") 122 n modulus) 123 doing 124 (multiple-value-bind (quot rem) 125 (truncate mod remainder) 126 (setf mod remainder 127 remainder rem) 128 (psetf a (- b (* a quot)) 129 b a)) 130 131 finally 132 ;; Since this isn't some general purpose Euclidean algorithm, but 133 ;; instead is trying to find a modulo inverse, we need to ensure that 134 ;; the Bezout coefficient we found (called A) is actually in [0, 135 ;; MODULUS). 136 ;; 137 ;; The general code calls CMOD here, but that doesn't know about the 138 ;; types of A and MODULUS, so we do it by hand, special-casing the easy 139 ;; case of modulus=2. 140 (return 141 (if (= mod-copy 2) 142 (logand a 1) 143 (let ((nn (mod a mod-copy))) 144 ;; nn here is in [0, modulus) 145 (if (<= (* 2 nn) mod-copy) 146 nn 147 (- nn mod-copy)))))))) 148 149;; CRECIP-GENERAL 150;; 151;; The general algorithm for CRECIP, valid when the modulus is any integer. See 152;; CRECIP for more details. 153(defun crecip-general (n) 154 ;; We assume that |n| < modulus, so n+modulus is always positive 155 (let ((mod modulus) 156 (remainder (if (minusp n) (+ n modulus) n)) 157 (a 1) (b 0)) 158 (loop 159 until (= remainder 1) 160 161 when (zerop remainder) do 162 (merror (intl:gettext "CRECIP: ~M does not have an inverse with modulus=~M") 163 n modulus) 164 doing 165 (let ((quotient (truncate mod remainder))) 166 (psetf mod remainder 167 remainder (- mod (* quotient remainder))) 168 (psetf a (- b (* a quotient)) 169 b a)) 170 171 finally (return (cmod a))))) 172 173;; CEXPT 174;; 175;; Raise an coefficient to a positive integral power. BASE should be a 176;; number. POW should be a non-negative integer. 177(defun cexpt (base pow) 178 (unless (typep pow '(integer 0)) 179 (error "CEXPT only defined for non-negative integral exponents.")) 180 (if (not modulus) 181 (expt base pow) 182 (do ((pow (ash pow -1) (ash pow -1)) 183 (s (if (oddp pow) base 1))) 184 ((zerop pow) s) 185 (setq base (rem (* base base) modulus)) 186 (when (oddp pow) (setq s (rem (* s base) modulus)))))) 187 188;; CMOD 189;; 190;; When MODULUS is null, this is the identity. Otherwise, it normalises N, which 191;; should be a number, to lie in the range (-modulus/2, modulus/2]. 192(defun cmod (n) 193 (declare (type number n)) 194 (if (not modulus) 195 n 196 (let ((rem (mod n modulus))) 197 (if (<= (* 2 rem) modulus) 198 rem 199 (- rem modulus))))) 200 201(defun cplus (a b) (cmod (+ a b))) 202(defun ctimes (a b) (cmod (* a b))) 203(defun cdifference (a b) (cmod (- a b))) 204 205;; SET-MODULUS 206;; 207;; Set the base in which the rational function package works. This does 208;; sanity-checking on the value chosen and is probably the way you should set 209;; the global value. 210;; 211;; Valid values for M are either a positive integer or NULL. 212(defun set-modulus (m) 213 (if (or (null m) (typep m '(integer 1))) 214 (setq modulus m) 215 (error "modulus must be a positive integer or nil")) 216 (values)) 217 218;; PCOEFADD 219;; 220;; Prepend a term to an existing polynomial. EXPONENT should be the exponent of 221;; the term to add; COEFF should be its coefficient; REMAINDER is a list of 222;; polynomial terms. The function returns polynomial terms that correspond to 223;; adding the given term. 224;; 225;; The function doesn't check that EXPONENT is higher than the highest exponent 226;; in REMAINDER, so you have to do this yourself. 227(defun pcoefadd (exponent coeff remainder) 228 (if (pzerop coeff) 229 remainder 230 (cons exponent (cons coeff remainder)))) 231 232;; PPLUS 233;; 234;; Add together two polynomials. 235(defun pplus (x y) 236 (cond ((pcoefp x) (pcplus x y)) 237 ((pcoefp y) (pcplus y x)) 238 ((eq (p-var x) (p-var y)) 239 (psimp (p-var x) (ptptplus (p-terms y) (p-terms x)))) 240 ((pointergp (p-var x) (p-var y)) 241 (psimp (p-var x) (ptcplus y (p-terms x)))) 242 (t (psimp (p-var y) (ptcplus x (p-terms y)))))) 243 244;; PTPTPLUS 245;; 246;; Add together two lists of polynomial terms. 247(defun ptptplus (x y) 248 (cond ((ptzerop x) y) 249 ((ptzerop y) x) 250 ((= (pt-le x) (pt-le y)) 251 (pcoefadd (pt-le x) 252 (pplus (pt-lc x) (pt-lc y)) 253 (ptptplus (pt-red x) (pt-red y)))) 254 ((> (pt-le x) (pt-le y)) 255 (cons (pt-le x) (cons (pt-lc x) (ptptplus (pt-red x) y)))) 256 (t (cons (pt-le y) (cons (pt-lc y) (ptptplus x (pt-red y))))))) 257 258;; PCPLUS 259;; 260;; Add a coefficient to a polynomial 261(defun pcplus (c p) 262 (if (pcoefp p) 263 (cplus p c) 264 (psimp (p-var p) 265 (ptcplus c (p-terms p))))) 266 267;; PTCPLUS 268;; 269;; Add a coefficient to a list of terms. C should be a used as a coefficient; 270;; TERMS is a list of a polynomial's terms. Note that we don't assume that C is 271;; a number: it might be a polynomial in a variable that isn't the main variable 272;; of the polynomial. 273(defun ptcplus (c terms) 274 (cond 275 ;; Adding zero doesn't do anything. 276 ((pzerop c) terms) 277 ;; Adding to zero, you just get the coefficient. 278 ((null terms) (list 0 c)) 279 ;; If terms are from a constant polynomial, we can just add C to its leading 280 ;; coefficient (which might not be a number in the multivariate case, so you 281 ;; have to use PPLUS) 282 ((zerop (pt-le terms)) 283 (pcoefadd 0 (pplus c (pt-lc terms)) nil)) 284 ;; If TERMS is a polynomial with degree > 0, recurse. 285 (t 286 (cons (pt-le terms) (cons (pt-lc terms) (ptcplus c (pt-red terms))))))) 287 288;; PDIFFERENCE 289;; 290;; Compute the difference of two polynomials 291(defun pdifference (x y) 292 (cond 293 ;; If Y is a coefficient, it's a number, so we can just add -Y to X using 294 ;; pcplus. If, however, X is the coefficient, we have to negate all the 295 ;; coefficients in Y, so defer to a utility function. 296 ((pcoefp x) (pcdiffer x y)) 297 ((pcoefp y) (pcplus (cminus y) x)) 298 ;; If X and Y have the same variable, work down their lists of terms. 299 ((eq (p-var x) (p-var y)) 300 (psimp (p-var x) (ptptdiffer (p-terms x) (p-terms y)))) 301 ;; Treat Y as a coefficient in the main variable of X. 302 ((pointergp (p-var x) (p-var y)) 303 (psimp (p-var x) (ptcdiffer-minus (p-terms x) y))) 304 ;; Treat X as a coefficient in the main variable of Y. 305 (t (psimp (p-var y) (ptcdiffer x (p-terms y)))))) 306 307;; PTPTDIFFER 308;; 309;; Compute the difference of two lists of polynomial terms (assumed to represent 310;; two polynomials in the same variable). 311(defun ptptdiffer (x y) 312 (cond 313 ((ptzerop x) (ptminus y)) 314 ((ptzerop y) x) 315 ((= (pt-le x) (pt-le y)) 316 (pcoefadd (pt-le x) 317 (pdifference (pt-lc x) (pt-lc y)) 318 (ptptdiffer (pt-red x) (pt-red y)))) 319 ((> (pt-le x) (pt-le y)) 320 (cons (pt-le x) (cons (pt-lc x) (ptptdiffer (pt-red x) y)))) 321 (t (cons (pt-le y) (cons (pminus (pt-lc y)) 322 (ptptdiffer x (pt-red y))))))) 323;; PCDIFFER 324;; 325;; Subtract the polynomial P from the coefficient C to form c - p. 326(defun pcdiffer (c p) 327 (if (pcoefp p) 328 (cdifference c p) 329 (psimp (p-var p) (ptcdiffer c (p-terms p))))) 330 331;; PTCDIFFER 332;; 333;; Subtract a polynomial represented by the list of terms, TERMS, from the 334;; coefficient C. 335(defun ptcdiffer (c terms) 336 (cond 337 ;; Unlike in the plus case or in PTCDIFFER-MINUS, we don't have a shortcut 338 ;; if C=0. However, if TERMS is null then we are calculating C-0, which is 339 ;; easy: 340 ((null terms) 341 (if (pzerop c) nil (list 0 c))) 342 ;; If the leading exponent is zero (in the main variable), then we can 343 ;; subtract the coefficients. Of course, these might actually be polynomials 344 ;; in other variables, so do this using pdifference. 345 ((zerop (pt-le terms)) 346 (pcoefadd 0 (pdifference c (pt-lc terms)) nil)) 347 ;; Otherwise we have to negate the leading coefficient (using pminus of 348 ;; course, because it might be a polynomial in other variables) and recurse. 349 (t 350 (cons (pt-le terms) 351 (cons (pminus (pt-lc terms)) (ptcdiffer c (pt-red terms))))))) 352 353;; PTCDIFFER-MINUS 354;; 355;; Subtract a coefficient, C, from a polynomial represented by a list of terms, 356;; TERMS, to form "p-c". This is the same as PTCDIFFER but the opposite sign (we 357;; don't implement it by (pminus (ptcdiffer c terms)) because that would require 358;; walking the polynomial twice) 359(defun ptcdiffer-minus (terms c) 360 (cond 361 ;; We're subtracting zero from a polynomial, which is easy! 362 ((pzerop c) terms) 363 ;; We're subtracting a coefficient from zero, which just comes out as the 364 ;; negation of the coefficient (compute it using pminus) 365 ((null terms) (list 0 (pminus c))) 366 ;; If the leading exponent is zero, subtract the coefficients just like in 367 ;; PTCDIFFER. 368 ((zerop (pt-le terms)) 369 (pcoefadd 0 (pdifference (pt-lc terms) c) nil)) 370 ;; Otherwise recurse. 371 (t 372 (cons (pt-le terms) 373 (cons (pt-lc terms) (ptcdiffer-minus (pt-red terms) c)))))) 374 375;; PCSUB 376;; 377;; Substitute values for variables in the polynomial P. VARS and VALS should be 378;; list of variables to substitute for and values to substitute, respectively. 379;; 380;; The code assumes that if VAR1 precedes VAR2 in the list then (POINTERGP VAR1 381;; VAR2). As such, VAR1 won't appear in the coefficients of a polynomial whose 382;; main variable is VAR2. 383(defun pcsub (p vals vars) 384 (cond 385 ;; Nothing to substitute, or P has no variables in it. 386 ((or (null vals) (pcoefp p)) p) 387 ;; The first variable in our list is the main variable of P. 388 ((eq (p-var p) (first vars)) 389 (ptcsub (p-terms p) (first vals) 390 (cdr vals) (cdr vars))) 391 ;; If the first var should appear before the main variable of P, we know it 392 ;; doesn't appear in any of the coefficients, so can (tail-)recurse on vals 393 ;; + vars. 394 ((pointergp (car vars) (p-var p)) 395 (pcsub p (cdr vals) (cdr vars))) 396 ;; Else, the main variable shouldn't get clobbered, but maybe we should 397 ;; replace variables in the coefficients. 398 (t (psimp (p-var p) (ptcsub-args (p-terms p) vals vars))))) 399 400;; PCSUBST 401;; 402;; Substitute VAL for VAR in a polynomial. Like PCSUB, but with only a single 403;; var to be substituted. 404;; 405;; (The logic of this function is exactly the same as PCSUB, but is marginally 406;; simpler because there are no more vars afterwards. Presumably, it was 407;; thought worth separating this case out from PCSUB to avoid spurious 408;; consing. I'm not convinced. RJS) 409(defun pcsubst (p val var) 410 (cond ((pcoefp p) p) 411 ((eq (p-var p) var) (ptcsub (cdr p) val nil nil)) 412 ((pointergp var (p-var p)) p) 413 (t (psimp (car p) (ptcsub-args (cdr p) (list val) (list var)))))) 414 415;; PTCSUB 416;; 417;; Substitute a constant, VAL, for the main variable in TERMS, which represent 418;; the terms of a polynomial. The coefficients might themselves be polynomials 419;; and, if so, we might substitute values for them too. To do so, pass VALS and 420;; VARS, with the same ordering requirements as in PCSUB. 421(defun ptcsub (terms val vals vars) 422 (if (eql val 0) 423 ;; If we're substituting 0 for the value, then we just extract the 424 ;; constant term. 425 (pcsub (ptterm terms 0) vals vars) 426 ;; Otherwise, walk through the polynomial using Horner's scheme to 427 ;; evaluate it. Because the polynomial is sparse, you can't just multiply 428 ;; by VAL every step, and instead have to keep track of the jump in 429 ;; exponents, which is what the LAST-LE variable does. 430 (do ((terms (pt-red terms) (pt-red terms)) 431 (ans (pcsub (pt-lc terms) vals vars) 432 (pplus (ptimes ans (pexpt val (- last-le (pt-le terms)))) 433 (pcsub (pt-lc terms) vals vars))) 434 (last-le (pt-le terms) (pt-le terms))) 435 ((null terms) 436 (ptimes ans (pexpt val last-le)))))) 437 438;; PTCSUB-ARGS 439;; 440;; Substitute values for vars in TERMS, which should be the terms of some 441;; polynomial. Unlike PTCSUB, we assume that the main variable of the polynomial 442;; isn't being substituted. VARS and VALS should be ordered as in PCSUB. 443(defun ptcsub-args (terms vals vars) 444 (loop 445 for (exp coef) on terms by #'cddr 446 unless (pzerop (setq coef (pcsub coef vals vars))) 447 nconc (list exp coef))) 448 449;; PCSUBSTY 450;; 451;; Like PCSUB, but with arguments in a different order and with a special case 452;; that you can pass atoms for VALS and VARS, in which case they will be treated 453;; as one-element lists. The big difference with PCSUB is that we don't assume 454;; that VARS and VALS come pre-sorted, and sort them here. 455(defun pcsubsty (vals vars p) 456 (cond 457 ;; If there is nothing to do, the answer is just P. 458 ((null vars) p) 459 ;; When there's only one variable, we don't need to do any sorting, so skip 460 ;; it and call PCSUB directly. 461 ((atom vars) (pcsub p (list vals) (list vars))) 462 ;; Otherwise, call PCSUB with a sorted list of VARS and VALS. 463 (t 464 (let ((pairs (sort (mapcar #'cons vars vals) #'pointergp :key #'car))) 465 (pcsub p (mapcar #'cdr pairs) (mapcar #'car pairs)))))) 466 467;; PDERIVATIVE 468;; 469;; Compute the derivative of the polynomial P with respect to the variable VARI. 470(defun pderivative (p vari) 471 (cond 472 ;; The derivative of a constant is zero. 473 ((pcoefp p) 0) 474 ;; If we have the same variable, do the differentiation term-by-term. 475 ((eq vari (p-var p)) 476 (psimp (p-var p) (ptderivative (p-terms p)))) 477 ;; If VARI > (P-VAR P) then we know it doesn't occur in any of the 478 ;; coefficients either, so return zero. This test comes after the one above 479 ;; because we expect more univariate polynomials and eq is cheaper than 480 ;; pointergp. 481 ((pointergp vari (p-var p)) 0) 482 ;; The other possibility is that (P-VAR P) > VARI, so the coefficients might 483 ;; need differentiating. 484 (t 485 (psimp (p-var p) (ptderivative-coeffs (p-terms p) vari))))) 486 487;; PTDERIVATIVE 488;; 489;; Formally differentiate TERMS, which is a list of the terms of some 490;; polynomial, with respect to that polynomial's main variable. 491(defun ptderivative (terms) 492 (if (or (null terms) (zerop (pt-le terms))) 493 ;; Zero or constant polynomials -> 0 494 nil 495 ;; Recurse, adding up "k . x^(k-1)" each time. 496 (pcoefadd (1- (pt-le terms)) 497 (pctimes (cmod (pt-le terms)) (pt-lc terms)) 498 (ptderivative (pt-red terms))))) 499 500;; PTDERIVATIVE-COEFFS 501;; 502;; Differentiate TERMS, which is a list of the terms of some polynomial, with 503;; respect to the variable VARI. We assume that VARI is not the main variable of 504;; the polynomial, but it might crop up in the coefficients. 505(defun ptderivative-coeffs (terms vari) 506 (and terms 507 ;; Recurse down the list of terms, calling PDERIVATIVE to actually 508 ;; differentiate each coefficient, then PTDERIVATIVE-COEFFS to do the rest. 509 (pcoefadd (pt-le terms) 510 (pderivative (pt-lc terms) vari) 511 (ptderivative-coeffs (pt-red terms) vari)))) 512 513;; PDIVIDE 514;; 515;; Polynomial division with remainder. X and Y should be polynomials. If V 516;; denotes the main variable of X, then we are carrying out the division in a 517;; ring of polynomials over Q where all variables that occur after V have been 518;; formally inverted. This is a Euclidean ring, and PDIVIDE implements division 519;; with remainder in this ring. 520;; 521;; The result is a list of two elements (Q R). Each is a rational function (a 522;; cons pair of polynomials), representing an element of F[V]. 523(defun pdivide (x y) 524 (cond 525 ((pzerop y) (rat-error "PDIVIDE: Quotient by zero")) 526 ;; If Y is a coefficient, it doesn't matter what X is: we can always do the 527 ;; division. 528 ((pacoefp y) (list (ratreduce x y) (rzero))) 529 ;; If X is a coefficient but Y isn't then the quotient must be zero 530 ((pacoefp x) (list (rzero) (cons x 1))) 531 ;; If neither is a coefficient then compare the variables. If V is greater 532 ;; than the main variable of Y, then Y is invertible in F[V]. 533 ((pointergp (p-var x) (p-var y)) (list (ratreduce x y) (rzero))) 534 ;; If we've got to here, V might occur in the coefficients of Y, but it 535 ;; needn't be the main variable. 536 (t 537 (do* ((lcy (cons (p-lc y) 1)) 538 (q (rzero)) 539 (r (cons x 1)) 540 (k (- (pdegree x (p-var y)) (p-le y)) 541 (- (pdegree (car r) (p-var y)) (p-le y)))) 542 543 ;; k is the degree of the numerator of the remainder minus the degree 544 ;; of y, both in the leading variable of y. For there to be further 545 ;; factors of y to subtract from q, this must be non-negative. 546 ((minusp k) (list q r)) 547 548 ;; Divide the leading coefficient of r (which means the leading term of 549 ;; the numerator, divided by the denominator) by the leading coefficient 550 ;; of y. 551 ;; 552 ;; The quotient gets added to q and gets multiplied back up by y and the 553 ;; result is subtracted from r. 554 (let* ((lcr (cons (p-lc (car r)) (cdr r))) 555 (quot (ratquotient lcr lcy)) 556 (quot-simp (cons (psimp (p-var y) (list k (car quot))) 557 (cdr quot)))) 558 (setf q (ratplus q quot-simp) 559 r (ratplus r (rattimes (cons (pminus y) 1) quot-simp t)))))))) 560 561;; PEXPT 562;; 563;; Polynomial exponentiation. Raise the polynomial P to the power N (which 564;; should be an integer) 565(defun pexpt (p n) 566 (cond 567 ;; p^0 = 1; p^1 = p 568 ((= n 0) 1) 569 ((= n 1) p) 570 ;; p^(-n) = 1/p^n 571 ((minusp n) (pquotient 1 (pexpt p (- n)))) 572 ;; When p is a coefficient, we can the simpler cexpt (which expects n >= 0, 573 ;; guaranteed by the previous clause) 574 ((pcoefp p) (cexpt p n)) 575 ;; If the main variable of P is an algebraic integer, calculate the power by 576 ;; repeated squaring (which will correctly take the remainder wrt the 577 ;; minimal polynomial for the variable) 578 ((alg p) (pexptsq p n)) 579 ;; If p is a monomial in the main variable, we're doing something like 580 ;; (x^2(y+1))^n, which is x^2n (y+1)^n, exponentiate the coefficient by 581 ;; recursion and just multiply the exponent. The call to PCOEFADD is just to 582 ;; ensure that we get zero if the coefficient raises to the power 583 ;; zero. (Possible when the coefficient is an algebraic integer) 584 ((null (p-red p)) 585 (psimp (p-var p) 586 (pcoefadd (* n (p-le p)) (pexpt (p-lc p) n) nil))) 587 ;; In the general case, expand using the binomial theorem. Write the 588 ;; calculation as 589 ;; 590 ;; (b + rest)^n = sum (binomial (n,k) * rest^k * b^(n-k), k, 0, n) 591 ;; 592 ;; We pre-compute a list of descending powers of B and use the formula 593 ;; 594 ;; binomial(n,k)/binomial(n,k-1) = (n+1-k) / k 595 ;; 596 ;; to keep track of the binomial coefficient. 597 (t 598 (let ((descending-powers (p-descending-powers 599 (make-poly (p-var p) (p-le p) (p-lc p)) n)) 600 (rest (psimp (p-var p) (p-red p)))) 601 (do* ((b-list descending-powers (rest b-list)) 602 (k 0 (1+ k)) 603 (n-choose-k 1 (truncate (* n-choose-k (- (1+ n) k)) k)) 604 (rest-pow 1 (case k 605 (1 rest) 606 (2 (pexpt rest 2)) 607 (t (ptimes rest rest-pow)))) 608 (sum (first descending-powers) 609 (pplus sum 610 (if b-list 611 (ptimes (pctimes (cmod n-choose-k) rest-pow) 612 (first b-list)) 613 (pctimes (cmod n-choose-k) rest-pow))))) 614 ((> k n) sum)))))) 615 616;; P-DESCENDING-POWERS 617;; 618;; Return a list of the powers of the polynomial P in descending order, starting 619;; with P^N and ending with P. 620(defun p-descending-powers (p n) 621 (let ((lst (list p))) 622 (dotimes (i (1- n)) (push (ptimes p (car lst)) lst)) 623 lst)) 624 625;; PMINUSP 626;; 627;; Returns true if the coefficient of the leading monomial of the polynomial is 628;; negative. Note that this depends on the variable ordering (for example, 629;; consider x-y). 630;; 631;; (pminusp '(y 1 -1 0 (x 1 1))) => T but 632;; (pminusp '(x 1 1 0 (y 1 -1))) => NIL 633(defun pminusp (p) 634 (if (realp p) (minusp p) 635 (pminusp (p-lc p)))) 636 637;; PMINUS 638;; 639;; Unary negation for polynomials. 640(defun pminus (p) 641 (if (pcoefp p) (cminus p) 642 (cons (p-var p) (ptminus (p-terms p))))) 643 644;; PTMINUS 645;; 646;; Negate a list of polynomial terms. 647(defun ptminus (x) 648 (loop for (exp coef) on x by #'cddr 649 nconc (list exp (pminus coef)))) 650 651;; PMOD 652;; 653;; Reduce a polynomial modulo the current value of MODULUS. 654(defun pmod (p) 655 (if (pcoefp p) (cmod p) 656 (psimp (car p) 657 (loop for (exp coef) on (p-terms p) by #'cddr 658 unless (pzerop (setq coef (pmod coef))) 659 nconc (list exp coef))))) 660 661;; PQUOTIENT 662;; 663;; Calculate x/y in the polynomial ring over the integers. Y should divide X 664;; without remainder. 665(defun pquotient (x y) 666 (cond ((pcoefp x) 667 (cond ((pzerop x) (pzero)) 668 ((pcoefp y) (cquotient x y)) 669 ((alg y) (paquo x y)) 670 (t (rat-error "PQUOTIENT: Quotient by a polynomial of higher degree (case 1)")))) 671 672 ((pcoefp y) 673 (cond ((pzerop y) (rat-error "PQUOTIENT: Quotient by zero")) 674 (modulus (pctimes (crecip y) x)) 675 (t (pcquotient x y)))) 676 677 ;; If (alg y) is true, then y is a polynomial in some variable that 678 ;; itself has a minimum polynomial. Moreover, the $algebraic flag must 679 ;; be true. We first try to compute an exact quotient ignoring that 680 ;; minimal polynomial, by binding $algebraic to nil. If that fails, we 681 ;; try to invert y and then multiply the results together. 682 ((alg y) (or (let ($algebraic) 683 (ignore-rat-err (pquotient x y))) 684 (patimes x (rainv y)))) 685 686 ;; If the main variable of Y comes after the main variable of X, Y must 687 ;; be free of that variable, so must divide each coefficient in X. Thus 688 ;; we can use PCQUOTIENT. 689 ((pointergp (p-var x) (p-var y)) (pcquotient x y)) 690 691 ;; Either Y contains a variable that is not in X, or they have the same 692 ;; main variable and Y has a higher degree. There can't possibly be an 693 ;; exact quotient. 694 ((pointergp (p-var y) (p-var x)) 695 (rat-error "PQUOTIENT: Quotient by a polynomial of higher degree (case 2a)")) 696 ((> (p-le y) (p-le x)) 697 (rat-error "PQUOTIENT: Quotient by a polynomial of higher degree (case 2b)")) 698 699 ;; If we got to here then X and Y have the same main variable and Y has 700 ;; a degree less than or equal to that of X. We can now forget about the 701 ;; main variable and work on the terms, with PTPTQUOTIENT. 702 (t 703 (psimp (p-var x) (ptptquotient (p-terms x) (p-terms y)))))) 704 705;; PCQUOTIENT 706;; 707;; Divide the polynomial P by Q. Q should be either a coefficient (so that 708;; (pcoefp q) => T), or should be a polynomial in a later variable than the main 709;; variable of P. Either way, Q is free of the main variable of P. The division 710;; is done at each coefficient. 711(defun pcquotient (p q) 712 (psimp (p-var p) 713 (loop 714 for (exp coef) on (p-terms p) by #'cddr 715 nconc (list exp (pquotient coef q))))) 716 717;; PTPTQUOTIENT 718;; 719;; Exactly divide two polynomials in the same variable, represented here by the 720;; list of their terms. 721(defun ptptquotient (u v) 722 ;; The algorithm is classic long division. You notice that if X/Y = Q then X = 723 ;; QY, so lc(X) = lc(Q)lc(Y) (where lc(Q)=Q when Q is a bare coefficient). Now 724 ;; divide again in the ring of coefficients to see that lc(X)/lc(Y) = 725 ;; lc(Q). Of course, you also know that le(Q) = le(X) - le(Y). 726 ;; 727 ;; Once you know lc(Q), you can subtract Y * lc(Q)*(var^le(Q)) from X and 728 ;; repeat. You know that you'll remove the leading term, so the algorithm will 729 ;; always terminate. To do the subtraction, use PTPT-SUBTRACT-POWERED-PRODUCT. 730 (do ((q-terms nil) 731 (u u (ptpt-subtract-powered-product (pt-red u) (pt-red v) 732 (first q-terms) (second q-terms)))) 733 ((ptzerop u) 734 (nreverse q-terms)) 735 ;; If B didn't divide A after all, then eventually we'll end up with the 736 ;; remainder in u, which has lower degree than that of B. 737 (when (< (pt-le u) (pt-le v)) 738 (rat-error "PTPTQUOTIENT: Polynomial quotient is not exact")) 739 (let ((le-q (- (pt-le u) (pt-le v))) 740 (lc-q (pquotient (pt-lc u) (pt-lc v)))) 741 ;; We've calculated the leading exponent and coefficient of q. Push them 742 ;; backwards onto q-terms (which holds the terms in reverse order). 743 (setf q-terms (cons lc-q (cons le-q q-terms)))))) 744 745;; PTPT-SUBTRACT-POWERED-PRODUCT 746;; 747;; U and V are the terms of two polynomials, A and B, in the same variable, x. Q 748;; is free of x. This function computes the terms of A - x^k * B * Q. This 749;; rather specialised function is used to update a numerator when doing 750;; polynomial long division. 751(defun ptpt-subtract-powered-product (u v q k) 752 (cond 753 ;; A - x^k * 0 * Q = A 754 ((null v) u) 755 ;; 0 - x^k * B * Q = x^k * B * (- Q) 756 ((null u) (pcetimes1 v k (pminus q))) 757 (t 758 ;; hipow is the highest exponent in x^k*B*Q. 759 (let ((hipow (+ (pt-le v) k))) 760 (cond 761 ;; If hipow is greater than the highest exponent in A, we have to 762 ;; prepend the first coefficient, which will be Q * lc(B). We can then 763 ;; recurse to this function to sort out the rest of the sum. 764 ((> hipow (pt-le u)) 765 (pcoefadd hipow 766 (ptimes q (pminus (pt-lc v))) 767 (ptpt-subtract-powered-product u (pt-red v) q k))) 768 ;; If hipow is equal to the highest exponent in A, we can just subtract 769 ;; the two leading coefficients and recurse to sort out the rest. 770 ((= hipow (pt-le u)) 771 (pcoefadd hipow 772 (pdifference (pt-lc u) (ptimes q (pt-lc v))) 773 (ptpt-subtract-powered-product (pt-red u) (pt-red v) q k))) 774 ;; If hipow is lower than the highest exponent in A then keep the first 775 ;; term of A and recurse. 776 (t 777 (list* (pt-le u) (pt-lc u) 778 (ptpt-subtract-powered-product (pt-red u) v q k)))))))) 779 780(defun algord (var) 781 (and $algebraic (get var 'algord))) 782 783;; PSIMP 784;; 785;; Return a "simplified" polynomial whose main variable is VAR and whose terms 786;; are given by X. 787;; 788;; If the polynomial is free of X, the result is the zero'th order coefficient: 789;; either a polynomial in later variables or a number. PSIMP also deals with 790;; reordering variables when $ALGEBRAIC is true, behaviour which is triggered by 791;; the ALGORD property on the main variable. 792(defun psimp (var x) 793 (cond ((ptzerop x) 0) 794 ((atom x) x) 795 ((zerop (pt-le x)) (pt-lc x)) 796 ((algord var) 797 ;; Fix wrong alg ordering: We deal with the case that the main variable 798 ;; of a coefficient should precede VAR. 799 (do ((p x (cddr p)) (sum 0)) 800 ((null p) 801 (if (pzerop sum) 802 (cons var x) 803 (pplus sum (p-delete-zeros var x)))) 804 ;; We only need to worry about the wrong ordering if a coefficient is 805 ;; a polynomial in another variable, and that variable should precede 806 ;; VAR. 807 (unless (or (pcoefp (pt-lc p)) 808 (pointergp var (p-var (pt-lc p)))) 809 (setq sum (pplus sum 810 (if (zerop (pt-le p)) (pt-lc p) 811 (ptimes (make-poly var (pt-le p) 1) 812 (pt-lc p))))) 813 ;; When we finish, we'll call PPLUS to add SUM and the remainder of 814 ;; X, and this line zeroes out this term in X (through P) to avoid 815 ;; double counting. The term will be deleted by the call to 816 ;; P-DELETE-ZEROS later. 817 (setf (pt-lc p) 0)))) 818 819 (t 820 (cons var x)))) 821 822;; P-DELETE-ZEROS 823;; 824;; Destructively operate on X, deleting any terms that have a zero coefficient. 825(defun p-delete-zeros (var x) 826 ;; The idea is that P always points one before the term in which we're 827 ;; interested. When that term has zero coefficient, it is trimmed from P by 828 ;; replacing the cdr. Consing NIL to the front of X allows us to throw away 829 ;; the first term if necessary. 830 (do ((p (setq x (cons nil x)))) 831 ((null (cdr p)) 832 ;; Switch off $algebraic so that we can recurse to PSIMP without any fear 833 ;; of an infinite recursion - PSIMP only calls this function when (ALGORD 834 ;; VAR) is true, and that only happens when $algebraic is true. 835 (let (($algebraic)) (psimp var (cdr x)))) 836 (if (pzerop (pt-lc (cdr p))) 837 (setf (cdr p) (pt-red (cdr p))) 838 (setq p (cddr p))))) 839 840;; PTTERM 841;; 842;; Given X representing the terms of a polynomial in a variable z, return the 843;; coefficient of z^n. 844(defun ptterm (x n) 845 (do ((x x (pt-red x))) 846 ((ptzerop x) (pzero)) 847 (cond ((< (pt-le x) n) (return (pzero))) 848 ((= (pt-le x) n) (return (pt-lc x)))))) 849 850(defun ptimes (x y) 851 (cond ((pcoefp x) (if (pzerop x) (pzero) (pctimes x y))) 852 ((pcoefp y) (if (pzerop y) (pzero) (pctimes y x))) 853 ((eq (p-var x) (p-var y)) 854 (palgsimp (p-var x) (ptimes1 (p-terms x) (p-terms y)) (alg x))) 855 ((pointergp (p-var x) (p-var y)) 856 (psimp (p-var x) (pctimes1 y (p-terms x)))) 857 (t (psimp (p-var y) (pctimes1 x (p-terms y)))))) 858 859(defun ptimes1 (x y-orig &aux uuu ) 860 (do ((vvv (setq uuu (pcetimes1 y-orig (pt-le x) (pt-lc x)))) 861 (x (pt-red x) (pt-red x))) 862 ((ptzerop x) uuu) 863 (let ((y y-orig) (xe (pt-le x)) (xc (pt-lc x))) 864 (prog (e u c) 865 a1 (cond ((null y) (return nil))) 866 (setq e (+ xe (car y))) 867 (setq c (ptimes (cadr y) xc)) 868 (cond ((pzerop c) (setq y (cddr y)) (go a1)) 869 ((or (null vvv) (> e (car vvv))) 870 (setq uuu (setq vvv (ptptplus uuu (list e c)))) 871 (setq y (cddr y)) (go a1)) 872 ((= e (car vvv)) 873 (setq c (pplus c (cadr vvv))) 874 (cond ((pzerop c) 875 (setq uuu (setq vvv (ptptdiffer uuu (list (car vvv) (cadr vvv)))))) 876 (t (rplaca (cdr vvv) c))) 877 (setq y (cddr y)) 878 (go a1))) 879 a 880 (cond ((and (cddr vvv) (> (caddr vvv) e)) 881 (setq vvv (cddr vvv)) (go a))) 882 (setq u (cdr vvv )) 883 b (cond ((or (null (cdr u)) (< (cadr u) e)) 884 (rplacd u (cons e (cons c (cdr u)))) (go e))) 885 (cond ((pzerop (setq c (pplus (caddr u) c))) 886 (rplacd u (cdddr u)) (go d)) 887 (t (rplaca (cddr u) c))) 888 e (setq u (cddr u)) 889 d (setq y (cddr y)) 890 (cond ((null y) (return nil))) 891 (setq e (+ xe (car y))) 892 (setq c (ptimes (cadr y) xc)) 893 c (cond ((and (cdr u) (> (cadr u) e)) (setq u (cddr u)) (go c))) 894 (go b)))) 895 uuu) 896 897(defun pcetimes1 (y e c) ;C*V^E*Y 898 (loop for (exp coef) on y by #'cddr 899 unless (pzerop (setq coef (ptimes c coef))) 900 nconc (list (+ e exp) coef))) 901 902(defun pctimes (c p) 903 (if (pcoefp p) (ctimes c p) 904 (psimp (p-var p) (pctimes1 c (p-terms p))))) 905 906(defun pctimes1 (c terms) 907 (loop for (exp coef) on terms by #'cddr 908 unless (pzerop (setq coef (ptimes c coef))) 909 nconc (list exp coef))) 910 911(defun leadalgcoef (p) 912 (cond ((pacoefp p) p) 913 (t (leadalgcoef (p-lc p))) )) 914 915(defun painvmod (q) 916 (cond ((pcoefp q) (crecip q)) 917 (t (paquo (list (car q) 0 1) q )))) 918 919(defun palgsimp (var p tell) ;TELL=(N X) -> X^(1/N) 920 (psimp var (cond ((or (null tell) (null p) 921 (< (car p) (car tell))) p) 922 ((null (cddr tell)) (pasimp1 p (car tell) (cadr tell))) 923 (t (pgcd1 p tell)) ))) 924 925(defun pasimp1 (p deg kernel) ;assumes deg>=(car p) 926 (do ((a p (pt-red a)) 927 (b p a)) 928 ((or (null a) (< (pt-le a) deg)) 929 (rplacd (cdr b) nil) 930 (ptptplus (pctimes1 kernel p) a)) 931 (rplaca a (- (pt-le a) deg)))) 932 933(defun monize (p) 934 (cond ((pcoefp p) (if (pzerop p) p 1)) 935 (t (cons (p-var p) (pmonicize (copy-list (p-terms p))))))) 936 937(defun pmonicize (p) ;CLOBBERS POLY 938 (cond ((equal (pt-lc p) 1) p) 939 (t (pmon1 (painvmod (leadalgcoef (pt-lc p))) p) p))) 940 941(defun pmon1 (mult l) 942 (cond (l (pmon1 mult (pt-red l)) 943 (setf (pt-lc l) (ptimes mult (pt-lc l)))))) 944 945(defun pmonz (poly &aux lc) ;A^(N-1)*P(X/A) 946 (setq poly (pabs poly)) 947 (cond ((equal (setq lc (p-lc poly)) 1) poly) 948 (t (do ((p (p-red poly) (pt-red p)) 949 (p1 (make-poly (p-var poly) (p-le poly) 1)) 950 (mult 1) 951 (deg (1- (p-le poly)) (pt-le p))) 952 ((null p) p1) 953 (setq mult (ptimes mult (pexpt lc (- deg (pt-le p))))) 954 (nconc p1 (list (pt-le p) (ptimes mult (pt-lc p)))))))) 955 956;; THESE ARE ROUTINES FOR MANIPULATING ALGEBRAIC NUMBERS 957 958(defun algnormal (p) (car (rquotient p (leadalgcoef p)))) 959 960(defun algcontent (p) 961 (destructuring-let* ((lcf (leadalgcoef p)) 962 ((prim . denom) (rquotient p lcf))) 963 (list (ratreduce lcf denom) prim))) 964 965(defun rquotient (p q &aux algfac* a e) ;FINDS PSEUDO QUOTIENT IF PSEUDOREM=0 966 (cond ((equal p q) (cons 1 1)) 967 ((pcoefp q) (ratreduce p q)) 968 ((setq a (testdivide p q)) (cons a 1)) 969 ((alg q) (rattimes (cons p 1) (rainv q) t)) 970 (t (cond ((alg (setq a (leadalgcoef q))) 971 (setq a (rainv a)) 972 (setq p (ptimes p (car a))) 973 (setq q (ptimes q (car a))) 974 (setq a (cdr a)) )) 975 (cond ((minusp (setq e (+ 1 (- (cadr q)) (pdegree p (car q))))) 976 (rat-error "RQUOTIENT: Quotient by a polynomial of higher degree"))) 977 (setq a (pexpt a e)) 978 (ratreduce (or (testdivide (ptimes a p) q) 979 (prog2 (setq a (pexpt (p-lc q) e)) 980 (pquotient (ptimes a p) q))) 981 a)) )) 982 983(defun patimes (x r) (pquotientchk (ptimes x (car r)) (cdr r))) 984 985(defun paquo (x y) (patimes x (rainv y))) 986 987(defun mpget (var) 988 (cond ((null (setq var (alg var))) nil) 989 ((cddr var) var) 990 (t (list (car var) 1 0 (pminus (cadr var)))))) 991 992 993(defun rainv (q) 994 (cond ((pcoefp q) 995 (cond (modulus (cons (crecip q) 1)) 996 (t (cons 1 q)))) 997 (t (let ((var (car q)) (p (mpget q))) 998 (declare (special var)) ;who uses this? --gsb 999 (cond ((null p) (cons 1 q)) 1000 (t (setq p (car (let ($ratalgdenom) 1001 (bprog q (cons var p))))) 1002 (rattimes (cons (car p) 1) (rainv (cdr p)) t))))))) 1003 1004(defun pexptsq (p n) 1005 (do ((n (ash n -1) (ash n -1)) 1006 (s (if (oddp n) p 1))) 1007 ((zerop n) s) 1008 (setq p (ptimes p p)) 1009 (and (oddp n) (setq s (ptimes s p))) )) 1010 1011;; THIS IS THE END OF THE NEW RATIONAL FUNCTION PACKAGE PART 1. 1012