Add note that the lapack package needs to loaded to get the functions.
[maxima.git] / src / rat3a.lisp
blob54a223c6e25d9b5b10c7bd0c70bebbefdbf3ef18
1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
4 ;;; ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
8 ;;; (c) Copyright 1980 Massachusetts Institute of Technology ;;;
9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
11 (in-package :maxima)
13 (macsyma-module rat3a)
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.
20 (declare-top (unspecial coef var exp p y))
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
26 ;; CQUOTIENT
28 ;; Calculate the quotient of two coefficients, which should be numbers. If
29 ;; MODULUS is non-nil, we try to take the reciprocal of A with respect to the
30 ;; modulus (using CRECIP) and then multiply by B. Note that this fails if B
31 ;; divides A as an integer, but B is not a unit in the ring of integers modulo
32 ;; MODULUS. For example,
34 ;; (let ((modulus 20)) (cquotient 10 5)) => ERROR
36 ;; If MODULUS is nil, then we work over the ring of integers when A and B are
37 ;; integers, and raise a RAT-ERROR if A is not divisible by B. If either A or B
38 ;; is a float then the division is done in floating point. Floats can get as far
39 ;; as the rational function code if $KEEPFLOAT is true.
41 ;; Before May 2023, this code used (rem a b) along with (/ a b) instead of
42 ;; just (floor a b). For Clozure CL the floor method is faster.
43 (defun cquotient (a b)
44 (cond ((eql a 0) 0) ;not sure this is OK--what if b=0 too?
45 ((null modulus)
46 (cond ((or (floatp a) (floatp b)) (/ a b)) ;not sure about floats here!
48 (multiple-value-bind (q r) (floor a b)
49 ;; when the remainder vanishes, return the quotient; else rat-error.
50 (if (eql 0 r) q (rat-error "CQUOTIENT: quotient is not exact"))))))
52 (ctimes a (crecip b)))))
53 ;; ALG
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))))
60 ;; PACOEFP
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))
68 ;; LEADTERM
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))))
76 ;; CRECIP
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]
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.
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)))
95 ;; And make sure that -MODULUS < N < MODULUS
96 (unless (<= (- modulus) n modulus)
97 (merror "N is out of range [-MODULUS, MODULUS] in crecip."))
99 ;; N in [0, MODULUS]
100 (when (minusp n) (setf n (+ n modulus)))
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))
117 (loop
118 until (= remainder 1)
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))
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).
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)
147 (- nn mod-copy))))))))
149 ;; CRECIP-GENERAL
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)
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))
171 finally (return (cmod a)))))
173 ;; CEXPT
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))))))
188 ;; CMOD
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)
196 (let ((rem (mod n modulus)))
197 (if (<= (* 2 rem) modulus)
199 (- rem modulus)))))
201 (defun cplus (a b) (cmod (+ a b)))
202 (defun ctimes (a b) (cmod (* a b)))
203 (defun cdifference (a b) (cmod (- a b)))
205 ;; SET-MODULUS
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.
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))
218 ;; PCOEFADD
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.
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))))
232 ;; PPLUS
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))))))
244 ;; PTPTPLUS
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)))))))
258 ;; PCPLUS
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)))))
267 ;; PTCPLUS
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.
286 (cons (pt-le terms) (cons (pt-lc terms) (ptcplus c (pt-red terms)))))))
288 ;; PDIFFERENCE
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))))))
307 ;; PTPTDIFFER
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
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)))))
331 ;; PTCDIFFER
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.
350 (cons (pt-le terms)
351 (cons (pminus (pt-lc terms)) (ptcdiffer c (pt-red terms)))))))
353 ;; PTCDIFFER-MINUS
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.
372 (cons (pt-le terms)
373 (cons (pt-lc terms) (ptcdiffer-minus (pt-red terms) c))))))
375 ;; PCSUB
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.
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)))))
400 ;; PCSUBST
402 ;; Substitute VAL for VAR in a polynomial. Like PCSUB, but with only a single
403 ;; var to be substituted.
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))))))
415 ;; PTCSUB
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))))))
438 ;; PTCSUB-ARGS
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)))
449 ;; PCSUBSTY
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.
464 (let ((pairs (sort (mapcar #'cons vars vals) #'pointergp :key #'car)))
465 (pcsub p (mapcar #'cdr pairs) (mapcar #'car pairs))))))
467 ;; PDERIVATIVE
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.
485 (psimp (p-var p) (ptderivative-coeffs (p-terms p) vari)))))
487 ;; PTDERIVATIVE
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
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)))))
500 ;; PTDERIVATIVE-COEFFS
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))))
513 ;; PDIVIDE
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.
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.
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))))
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))
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.
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))))))))
561 ;; PEXPT
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
590 ;; (b + rest)^n = sum (binomial (n,k) * rest^k * b^(n-k), k, 0, n)
592 ;; We pre-compute a list of descending powers of B and use the formula
594 ;; binomial(n,k)/binomial(n,k-1) = (n+1-k) / k
596 ;; to keep track of the binomial coefficient.
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))))))
616 ;; P-DESCENDING-POWERS
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))
625 ;; PMINUSP
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).
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))))
637 ;; PMINUS
639 ;; Unary negation for polynomials.
640 (defun pminus (p)
641 (if (pcoefp p) (cminus p)
642 (cons (p-var p) (ptminus (p-terms p)))))
644 ;; PTMINUS
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))))
651 ;; PMOD
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)))))
661 ;; PQUOTIENT
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)"))))
672 ((pcoefp y)
673 (cond ((pzerop y) (rat-error "PQUOTIENT: Quotient by zero"))
674 (modulus (pctimes (crecip y) x))
675 (t (pcquotient x y))))
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))))
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))
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)"))
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.
703 (psimp (p-var x) (ptptquotient (p-terms x) (p-terms y))))))
705 ;; PCQUOTIENT
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)))))
717 ;; PTPTQUOTIENT
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).
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))))))
745 ;; PTPT-SUBTRACT-POWERED-PRODUCT
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)))
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.
777 (list* (pt-le u) (pt-lc u)
778 (ptpt-subtract-powered-product (pt-red u) v q k))))))))
780 (defun algord (var)
781 (and $algebraic (get var 'algord)))
783 ;; PSIMP
785 ;; Return a "simplified" polynomial whose main variable is VAR and whose terms
786 ;; are given by X.
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))))
820 (cons var x))))
822 ;; P-DELETE-ZEROS
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)))))
840 ;; PTTERM
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))))))
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))))))
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)))
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)
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)))
902 (defun pctimes (c p)
903 (if (pcoefp p) (ctimes c p)
904 (psimp (p-var p) (pctimes1 c (p-terms p)))))
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)))
911 (defun leadalgcoef (p)
912 (cond ((pacoefp p) p)
913 (t (leadalgcoef (p-lc p))) ))
915 (defun painvmod (q)
916 (cond ((pcoefp q) (crecip q))
917 (t (paquo (list (car q) 0 1) q ))))
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)) )))
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))))
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)))))))
937 (defun pmonicize (p) ;CLOBBERS POLY
938 (cond ((equal (pt-lc p) 1) p)
939 (t (pmon1 (painvmod (leadalgcoef (pt-lc p))) p) p)))
941 (defun pmon1 (mult l)
942 (cond (l (pmon1 mult (pt-red l))
943 (setf (pt-lc l) (ptimes mult (pt-lc l))))))
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))))))))
956 ;; THESE ARE ROUTINES FOR MANIPULATING ALGEBRAIC NUMBERS
958 (defun algnormal (p) (car (rquotient p (leadalgcoef p))))
960 (defun algcontent (p)
961 (destructuring-let* ((lcf (leadalgcoef p))
962 ((prim . denom) (rquotient p lcf)))
963 (list (ratreduce lcf denom) prim)))
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)) ))
983 (defun patimes (x r) (pquotientchk (ptimes x (car r)) (cdr r)))
985 (defun paquo (x y) (patimes x (rainv y)))
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))))))
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) var))))
1002 (rattimes (cons (car p) 1) (rainv (cdr p)) t)))))))
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))) ))
1011 ;; THIS IS THE END OF THE NEW RATIONAL FUNCTION PACKAGE PART 1.