1 ;;; calc-poly.el --- polynomial functions for Calc
3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
4 ;; 2005, 2006, 2007, 2008, 2009, 2010 Free Software Foundation, Inc.
6 ;; Author: David Gillespie <daveg@synaptics.com>
7 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
9 ;; This file is part of GNU Emacs.
11 ;; GNU Emacs is free software: you can redistribute it and/or modify
12 ;; it under the terms of the GNU General Public License as published by
13 ;; the Free Software Foundation, either version 3 of the License, or
14 ;; (at your option) any later version.
16 ;; GNU Emacs is distributed in the hope that it will be useful,
17 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
18 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 ;; GNU General Public License for more details.
21 ;; You should have received a copy of the GNU General Public License
22 ;; along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>.
28 ;; This file is autoloaded from calc-ext.el.
33 (defun calcFunc-pcont (expr &optional var
)
34 (cond ((Math-primp expr
)
35 (cond ((Math-zerop expr
) 1)
36 ((Math-messy-integerp expr
) (math-trunc expr
))
37 ((Math-objectp expr
) expr
)
38 ((or (equal expr var
) (not var
)) 1)
41 (math-mul (calcFunc-pcont (nth 1 expr
) var
)
42 (calcFunc-pcont (nth 2 expr
) var
)))
44 (math-div (calcFunc-pcont (nth 1 expr
) var
)
45 (calcFunc-pcont (nth 2 expr
) var
)))
46 ((and (eq (car expr
) '^
) (Math-natnump (nth 2 expr
)))
47 (math-pow (calcFunc-pcont (nth 1 expr
) var
) (nth 2 expr
)))
48 ((memq (car expr
) '(neg polar
))
49 (calcFunc-pcont (nth 1 expr
) var
))
51 (let ((p (math-is-polynomial expr var
)))
53 (let ((lead (nth (1- (length p
)) p
))
54 (cont (math-poly-gcd-list p
)))
55 (if (math-guess-if-neg lead
)
59 ((memq (car expr
) '(+ - cplx sdev
))
60 (let ((cont (calcFunc-pcont (nth 1 expr
) var
)))
63 (let ((c2 (calcFunc-pcont (nth 2 expr
) var
)))
64 (if (and (math-negp cont
)
65 (if (eq (car expr
) '-
) (math-posp c2
) (math-negp c2
)))
66 (math-neg (math-poly-gcd cont c2
))
67 (math-poly-gcd cont c2
))))))
71 (defun calcFunc-pprim (expr &optional var
)
72 (let ((cont (calcFunc-pcont expr var
)))
73 (if (math-equal-int cont
1)
75 (math-poly-div-exact expr cont var
))))
77 (defun math-div-poly-const (expr c
)
78 (cond ((memq (car-safe expr
) '(+ -
))
80 (math-div-poly-const (nth 1 expr
) c
)
81 (math-div-poly-const (nth 2 expr
) c
)))
82 (t (math-div expr c
))))
84 (defun calcFunc-pdeg (expr &optional var
)
86 '(neg (var inf var-inf
))
88 (or (math-polynomial-p expr var
)
89 (math-reject-arg expr
"Expected a polynomial"))
90 (math-poly-degree expr
))))
92 (defun math-poly-degree (expr)
93 (cond ((Math-primp expr
)
94 (if (eq (car-safe expr
) 'var
) 1 0))
96 (math-poly-degree (nth 1 expr
)))
98 (+ (math-poly-degree (nth 1 expr
))
99 (math-poly-degree (nth 2 expr
))))
101 (- (math-poly-degree (nth 1 expr
))
102 (math-poly-degree (nth 2 expr
))))
103 ((and (eq (car expr
) '^
) (natnump (nth 2 expr
)))
104 (* (math-poly-degree (nth 1 expr
)) (nth 2 expr
)))
105 ((memq (car expr
) '(+ -
))
106 (max (math-poly-degree (nth 1 expr
))
107 (math-poly-degree (nth 2 expr
))))
110 (defun calcFunc-plead (expr var
)
111 (cond ((eq (car-safe expr
) '*)
112 (math-mul (calcFunc-plead (nth 1 expr
) var
)
113 (calcFunc-plead (nth 2 expr
) var
)))
114 ((eq (car-safe expr
) '/)
115 (math-div (calcFunc-plead (nth 1 expr
) var
)
116 (calcFunc-plead (nth 2 expr
) var
)))
117 ((and (eq (car-safe expr
) '^
) (math-natnump (nth 2 expr
)))
118 (math-pow (calcFunc-plead (nth 1 expr
) var
) (nth 2 expr
)))
124 (let ((p (math-is-polynomial expr var
)))
126 (nth (1- (length p
)) p
)
133 ;;; Polynomial quotient, remainder, and GCD.
134 ;;; Originally by Ove Ewerlid (ewerlid@mizar.DoCS.UU.SE).
135 ;;; Modifications and simplifications by daveg.
137 (defvar math-poly-modulus
1)
139 ;;; Return gcd of two polynomials
140 (defun calcFunc-pgcd (pn pd
)
141 (if (math-any-floats pn
)
142 (math-reject-arg pn
"Coefficients must be rational"))
143 (if (math-any-floats pd
)
144 (math-reject-arg pd
"Coefficients must be rational"))
145 (let ((calc-prefer-frac t
)
146 (math-poly-modulus (math-poly-modulus pn pd
)))
147 (math-poly-gcd pn pd
)))
149 ;;; Return only quotient to top of stack (nil if zero)
151 ;; calc-poly-div-remainder is a local variable for
152 ;; calc-poly-div (in calc-alg.el), but is used by
153 ;; calcFunc-pdiv, which is called by calc-poly-div.
154 (defvar calc-poly-div-remainder
)
156 (defun calcFunc-pdiv (pn pd
&optional base
)
157 (let* ((calc-prefer-frac t
)
158 (math-poly-modulus (math-poly-modulus pn pd
))
159 (res (math-poly-div pn pd base
)))
160 (setq calc-poly-div-remainder
(cdr res
))
163 ;;; Return only remainder to top of stack
164 (defun calcFunc-prem (pn pd
&optional base
)
165 (let ((calc-prefer-frac t
)
166 (math-poly-modulus (math-poly-modulus pn pd
)))
167 (cdr (math-poly-div pn pd base
))))
169 (defun calcFunc-pdivrem (pn pd
&optional base
)
170 (let* ((calc-prefer-frac t
)
171 (math-poly-modulus (math-poly-modulus pn pd
))
172 (res (math-poly-div pn pd base
)))
173 (list 'vec
(car res
) (cdr res
))))
175 (defun calcFunc-pdivide (pn pd
&optional base
)
176 (let* ((calc-prefer-frac t
)
177 (math-poly-modulus (math-poly-modulus pn pd
))
178 (res (math-poly-div pn pd base
)))
179 (math-add (car res
) (math-div (cdr res
) pd
))))
182 ;;; Multiply two terms, expanding out products of sums.
183 (defun math-mul-thru (lhs rhs
)
184 (if (memq (car-safe lhs
) '(+ -
))
186 (math-mul-thru (nth 1 lhs
) rhs
)
187 (math-mul-thru (nth 2 lhs
) rhs
))
188 (if (memq (car-safe rhs
) '(+ -
))
190 (math-mul-thru lhs
(nth 1 rhs
))
191 (math-mul-thru lhs
(nth 2 rhs
)))
192 (math-mul lhs rhs
))))
194 (defun math-div-thru (num den
)
195 (if (memq (car-safe num
) '(+ -
))
197 (math-div-thru (nth 1 num
) den
)
198 (math-div-thru (nth 2 num
) den
))
202 ;;; Sort the terms of a sum into canonical order.
203 (defun math-sort-terms (expr)
204 (if (memq (car-safe expr
) '(+ -
))
206 (sort (math-sum-to-list expr
)
207 (function (lambda (a b
) (math-beforep (car a
) (car b
))))))
210 (defun math-list-to-sum (lst)
212 (list (if (cdr (car lst
)) '-
'+)
213 (math-list-to-sum (cdr lst
))
216 (math-neg (car (car lst
)))
219 (defun math-sum-to-list (tree &optional neg
)
220 (cond ((eq (car-safe tree
) '+)
221 (nconc (math-sum-to-list (nth 1 tree
) neg
)
222 (math-sum-to-list (nth 2 tree
) neg
)))
223 ((eq (car-safe tree
) '-
)
224 (nconc (math-sum-to-list (nth 1 tree
) neg
)
225 (math-sum-to-list (nth 2 tree
) (not neg
))))
226 (t (list (cons tree neg
)))))
228 ;;; Check if the polynomial coefficients are modulo forms.
229 (defun math-poly-modulus (expr &optional expr2
)
230 (or (math-poly-modulus-rec expr
)
231 (and expr2
(math-poly-modulus-rec expr2
))
234 (defun math-poly-modulus-rec (expr)
235 (if (and (eq (car-safe expr
) 'mod
) (Math-natnump (nth 2 expr
)))
236 (list 'mod
1 (nth 2 expr
))
237 (and (memq (car-safe expr
) '(+ -
* /))
238 (or (math-poly-modulus-rec (nth 1 expr
))
239 (math-poly-modulus-rec (nth 2 expr
))))))
242 ;;; Divide two polynomials. Return (quotient . remainder).
243 (defvar math-poly-div-base nil
)
244 (defun math-poly-div (u v
&optional math-poly-div-base
)
245 (if math-poly-div-base
246 (math-do-poly-div u v
)
247 (math-do-poly-div (calcFunc-expand u
) (calcFunc-expand v
))))
249 (defun math-poly-div-exact (u v
&optional base
)
250 (let ((res (math-poly-div u v base
)))
253 (math-reject-arg (list 'vec u v
) "Argument is not a polynomial"))))
255 (defun math-do-poly-div (u v
)
256 (cond ((math-constp u
)
258 (cons (math-div u v
) 0)
263 (if (memq (car-safe u
) '(+ -
))
264 (math-add-or-sub (math-poly-div-exact (nth 1 u
) v
)
265 (math-poly-div-exact (nth 2 u
) v
)
270 (cons math-poly-modulus
0))
271 ((and (math-atomic-factorp u
) (math-atomic-factorp v
))
272 (cons (math-simplify (math-div u v
)) 0))
274 (let ((base (or math-poly-div-base
275 (math-poly-div-base u v
)))
278 (null (setq vp
(math-is-polynomial v base nil
'gen
))))
280 (setq up
(math-is-polynomial u base nil
'gen
)
281 res
(math-poly-div-coefs up vp
))
282 (cons (math-build-polynomial-expr (car res
) base
)
283 (math-build-polynomial-expr (cdr res
) base
)))))))
285 (defun math-poly-div-rec (u v
)
286 (cond ((math-constp u
)
291 (if (memq (car-safe u
) '(+ -
))
292 (math-add-or-sub (math-poly-div-rec (nth 1 u
) v
)
293 (math-poly-div-rec (nth 2 u
) v
)
296 ((Math-equal u v
) math-poly-modulus
)
297 ((and (math-atomic-factorp u
) (math-atomic-factorp v
))
298 (math-simplify (math-div u v
)))
302 (let ((base (math-poly-div-base u v
))
305 (null (setq vp
(math-is-polynomial v base nil
'gen
))))
307 (setq up
(math-is-polynomial u base nil
'gen
)
308 res
(math-poly-div-coefs up vp
))
309 (math-add (math-build-polynomial-expr (car res
) base
)
310 (math-div (math-build-polynomial-expr (cdr res
) base
)
313 ;;; Divide two polynomials in coefficient-list form. Return (quot . rem).
314 (defun math-poly-div-coefs (u v
)
315 (cond ((null v
) (math-reject-arg nil
"Division by zero"))
316 ((< (length u
) (length v
)) (cons nil u
))
322 (let ((qk (math-poly-div-rec (math-simplify (car urev
))
326 (if (or q
(not (math-zerop qk
)))
327 (setq q
(cons qk q
)))
328 (while (setq up
(cdr up
) vp
(cdr vp
))
329 (setcar up
(math-sub (car up
) (math-mul-thru qk
(car vp
)))))
330 (setq urev
(cdr urev
))
332 (while (and urev
(Math-zerop (car urev
)))
333 (setq urev
(cdr urev
)))
334 (cons q
(nreverse (mapcar 'math-simplify urev
)))))
336 (cons (list (math-poly-div-rec (car u
) (car v
)))
339 ;;; Perform a pseudo-division of polynomials. (See Knuth section 4.6.1.)
340 ;;; This returns only the remainder from the pseudo-division.
341 (defun math-poly-pseudo-div (u v
)
343 ((< (length u
) (length v
)) u
)
344 ((or (cdr u
) (cdr v
))
345 (let ((urev (reverse u
))
351 (while (setq up
(cdr up
) vp
(cdr vp
))
352 (setcar up
(math-sub (math-mul-thru (car vrev
) (car up
))
353 (math-mul-thru (car urev
) (car vp
)))))
354 (setq urev
(cdr urev
))
357 (setcar up
(math-mul-thru (car vrev
) (car up
)))
359 (while (and urev
(Math-zerop (car urev
)))
360 (setq urev
(cdr urev
)))
361 (nreverse (mapcar 'math-simplify urev
))))
364 ;;; Compute the GCD of two multivariate polynomials.
365 (defun math-poly-gcd (u v
)
366 (cond ((Math-equal u v
) u
)
370 (calcFunc-gcd u
(calcFunc-pcont v
))))
374 (calcFunc-gcd v
(calcFunc-pcont u
))))
376 (let ((base (math-poly-gcd-base u v
)))
380 (math-build-polynomial-expr
381 (math-poly-gcd-coefs (math-is-polynomial u base nil
'gen
)
382 (math-is-polynomial v base nil
'gen
))
384 (calcFunc-gcd (calcFunc-pcont u
) (calcFunc-pcont u
)))))))
386 (defun math-poly-div-list (lst a
)
390 (math-mul-list lst a
)
391 (mapcar (function (lambda (x) (math-poly-div-exact x a
))) lst
))))
393 (defun math-mul-list (lst a
)
397 (mapcar 'math-neg lst
)
399 (mapcar (function (lambda (x) (math-mul x a
))) lst
)))))
401 ;;; Run GCD on all elements in a list.
402 (defun math-poly-gcd-list (lst)
403 (if (or (memq 1 lst
) (memq -
1 lst
))
404 (math-poly-gcd-frac-list lst
)
405 (let ((gcd (car lst
)))
406 (while (and (setq lst
(cdr lst
)) (not (eq gcd
1)))
408 (setq gcd
(math-poly-gcd gcd
(car lst
)))))
409 (if lst
(setq lst
(math-poly-gcd-frac-list lst
)))
412 (defun math-poly-gcd-frac-list (lst)
413 (while (and lst
(not (eq (car-safe (car lst
)) 'frac
)))
414 (setq lst
(cdr lst
)))
416 (let ((denom (nth 2 (car lst
))))
417 (while (setq lst
(cdr lst
))
418 (if (eq (car-safe (car lst
)) 'frac
)
419 (setq denom
(calcFunc-lcm denom
(nth 2 (car lst
))))))
420 (list 'frac
1 denom
))
423 ;;; Compute the GCD of two monovariate polynomial lists.
424 ;;; Knuth section 4.6.1, algorithm C.
425 (defun math-poly-gcd-coefs (u v
)
426 (let ((d (math-poly-gcd (math-poly-gcd-list u
)
427 (math-poly-gcd-list v
)))
428 (g 1) (h 1) (z 0) hh r delta ghd
)
429 (while (and u v
(Math-zerop (car u
)) (Math-zerop (car v
)))
430 (setq u
(cdr u
) v
(cdr v
) z
(1+ z
)))
432 (setq u
(math-poly-div-list u d
)
433 v
(math-poly-div-list v d
)))
435 (setq delta
(- (length u
) (length v
)))
437 (setq r u u v v r delta
(- delta
)))
438 (setq r
(math-poly-pseudo-div u v
))
441 v
(math-poly-div-list r
(math-mul g
(math-pow h delta
)))
442 g
(nth (1- (length u
)) u
)
444 (math-mul (math-pow g delta
) (math-pow h
(- 1 delta
)))
445 (math-poly-div-exact (math-pow g delta
)
446 (math-pow h
(1- delta
))))))
449 (math-mul-list (math-poly-div-list v
(math-poly-gcd-list v
)) d
)))
450 (if (math-guess-if-neg (nth (1- (length v
)) v
))
451 (setq v
(math-mul-list v -
1)))
452 (while (>= (setq z
(1- z
)) 0)
457 ;;; Return true if is a factor containing no sums or quotients.
458 (defun math-atomic-factorp (expr)
459 (cond ((eq (car-safe expr
) '*)
460 (and (math-atomic-factorp (nth 1 expr
))
461 (math-atomic-factorp (nth 2 expr
))))
462 ((memq (car-safe expr
) '(+ -
/))
464 ((memq (car-safe expr
) '(^ neg
))
465 (math-atomic-factorp (nth 1 expr
)))
468 ;;; Find a suitable base for dividing a by b.
469 ;;; The base must exist in both expressions.
470 ;;; The degree in the numerator must be higher or equal than the
471 ;;; degree in the denominator.
472 ;;; If the above conditions are not met the quotient is just a remainder.
473 ;;; Return nil if this is the case.
475 (defun math-poly-div-base (a b
)
477 (and (setq a-base
(math-total-polynomial-base a
))
478 (setq b-base
(math-total-polynomial-base b
))
481 (let ((maybe (assoc (car (car a-base
)) b-base
)))
483 (if (>= (nth 1 (car a-base
)) (nth 1 maybe
))
484 (throw 'return
(car (car a-base
))))))
485 (setq a-base
(cdr a-base
)))))))
487 ;;; Same as above but for gcd algorithm.
488 ;;; Here there is no requirement that degree(a) > degree(b).
489 ;;; Take the base that has the highest degree considering both a and b.
490 ;;; ("a^20+b^21+x^3+a+b", "a+b^2+x^5+a^22+b^10") --> (a 22)
492 (defun math-poly-gcd-base (a b
)
494 (and (setq a-base
(math-total-polynomial-base a
))
495 (setq b-base
(math-total-polynomial-base b
))
497 (while (and a-base b-base
)
498 (if (> (nth 1 (car a-base
)) (nth 1 (car b-base
)))
499 (if (assoc (car (car a-base
)) b-base
)
500 (throw 'return
(car (car a-base
)))
501 (setq a-base
(cdr a-base
)))
502 (if (assoc (car (car b-base
)) a-base
)
503 (throw 'return
(car (car b-base
)))
504 (setq b-base
(cdr b-base
)))))))))
506 ;;; Sort a list of polynomial bases.
507 (defun math-sort-poly-base-list (lst)
508 (sort lst
(function (lambda (a b
)
509 (or (> (nth 1 a
) (nth 1 b
))
510 (and (= (nth 1 a
) (nth 1 b
))
511 (math-beforep (car a
) (car b
))))))))
513 ;;; Given an expression find all variables that are polynomial bases.
514 ;;; Return list in the form '( (var1 degree1) (var2 degree2) ... ).
516 ;; The variable math-poly-base-total-base is local to
517 ;; math-total-polynomial-base, but is used by math-polynomial-p1,
518 ;; which is called by math-total-polynomial-base.
519 (defvar math-poly-base-total-base
)
521 (defun math-total-polynomial-base (expr)
522 (let ((math-poly-base-total-base nil
))
523 (math-polynomial-base expr
'math-polynomial-p1
)
524 (math-sort-poly-base-list math-poly-base-total-base
)))
526 ;; The variable math-poly-base-top-expr is local to math-polynomial-base
527 ;; in calc-alg.el, but is used by math-polynomial-p1 which is called
528 ;; by math-polynomial-base.
529 (defvar math-poly-base-top-expr
)
531 (defun math-polynomial-p1 (subexpr)
532 (or (assoc subexpr math-poly-base-total-base
)
533 (memq (car subexpr
) '(+ -
* / neg
))
534 (and (eq (car subexpr
) '^
) (natnump (nth 2 subexpr
)))
535 (let* ((math-poly-base-variable subexpr
)
536 (exponent (math-polynomial-p math-poly-base-top-expr subexpr
)))
538 (setq math-poly-base-total-base
(cons (list subexpr exponent
)
539 math-poly-base-total-base
)))))
542 ;; The variable math-factored-vars is local to calcFunc-factors and
543 ;; calcFunc-factor, but is used by math-factor-expr and
544 ;; math-factor-expr-part, which are called (directly and indirectly) by
545 ;; calcFunc-factor and calcFunc-factors.
546 (defvar math-factored-vars
)
548 ;; The variable math-fact-expr is local to calcFunc-factors,
549 ;; calcFunc-factor and math-factor-expr, but is used by math-factor-expr-try
550 ;; and math-factor-expr-part, which are called (directly and indirectly) by
551 ;; calcFunc-factor, calcFunc-factors and math-factor-expr.
552 (defvar math-fact-expr
)
554 ;; The variable math-to-list is local to calcFunc-factors and
555 ;; calcFunc-factor, but is used by math-accum-factors, which is
556 ;; called (indirectly) by calcFunc-factors and calcFunc-factor.
557 (defvar math-to-list
)
559 (defun calcFunc-factors (math-fact-expr &optional var
)
560 (let ((math-factored-vars (if var t nil
))
562 (calc-prefer-frac t
))
564 (setq var
(math-polynomial-base math-fact-expr
)))
565 (let ((res (math-factor-finish
566 (or (catch 'factor
(math-factor-expr-try var
))
568 (math-simplify (if (math-vectorp res
)
570 (list 'vec
(list 'vec res
1)))))))
572 (defun calcFunc-factor (math-fact-expr &optional var
)
573 (let ((math-factored-vars nil
)
575 (calc-prefer-frac t
))
576 (math-simplify (math-factor-finish
578 (let ((math-factored-vars t
))
579 (or (catch 'factor
(math-factor-expr-try var
)) math-fact-expr
))
580 (math-factor-expr math-fact-expr
))))))
582 (defun math-factor-finish (x)
585 (if (eq (car x
) 'calcFunc-Fac-Prot
)
586 (math-factor-finish (nth 1 x
))
587 (cons (car x
) (mapcar 'math-factor-finish
(cdr x
))))))
589 (defun math-factor-protect (x)
590 (if (memq (car-safe x
) '(+ -
))
591 (list 'calcFunc-Fac-Prot x
)
594 (defun math-factor-expr (math-fact-expr)
595 (cond ((eq math-factored-vars t
) math-fact-expr
)
596 ((or (memq (car-safe math-fact-expr
) '(* / ^ neg
))
597 (assq (car-safe math-fact-expr
) calc-tweak-eqn-table
))
598 (cons (car math-fact-expr
) (mapcar 'math-factor-expr
(cdr math-fact-expr
))))
599 ((memq (car-safe math-fact-expr
) '(+ -
))
600 (let* ((math-factored-vars math-factored-vars
)
601 (y (catch 'factor
(math-factor-expr-part math-fact-expr
))))
607 (defun math-factor-expr-part (x) ; uses "expr"
608 (if (memq (car-safe x
) '(+ -
* / ^ neg
))
609 (while (setq x
(cdr x
))
610 (math-factor-expr-part (car x
)))
611 (and (not (Math-objvecp x
))
612 (not (assoc x math-factored-vars
))
613 (> (math-factor-contains math-fact-expr x
) 1)
614 (setq math-factored-vars
(cons (list x
) math-factored-vars
))
615 (math-factor-expr-try x
))))
617 ;; The variable math-fet-x is local to math-factor-expr-try, but is
618 ;; used by math-factor-poly-coefs, which is called by math-factor-expr-try.
621 (defun math-factor-expr-try (math-fet-x)
622 (if (eq (car-safe math-fact-expr
) '*)
623 (let ((res1 (catch 'factor
(let ((math-fact-expr (nth 1 math-fact-expr
)))
624 (math-factor-expr-try math-fet-x
))))
625 (res2 (catch 'factor
(let ((math-fact-expr (nth 2 math-fact-expr
)))
626 (math-factor-expr-try math-fet-x
)))))
628 (throw 'factor
(math-accum-factors (or res1
(nth 1 math-fact-expr
)) 1
629 (or res2
(nth 2 math-fact-expr
))))))
630 (let* ((p (math-is-polynomial math-fact-expr math-fet-x
30 'gen
))
631 (math-poly-modulus (math-poly-modulus math-fact-expr
))
634 (setq res
(math-factor-poly-coefs p
))
635 (throw 'factor res
)))))
637 (defun math-accum-factors (fac pow facs
)
639 (if (math-vectorp fac
)
641 (while (setq fac
(cdr fac
))
642 (setq facs
(math-accum-factors (nth 1 (car fac
))
643 (* pow
(nth 2 (car fac
)))
646 (if (and (eq (car-safe fac
) '^
) (natnump (nth 2 fac
)))
647 (setq pow
(* pow
(nth 2 fac
))
651 (or (math-vectorp facs
)
652 (setq facs
(if (eq facs
1) '(vec)
653 (list 'vec
(list 'vec facs
1)))))
655 (while (and (setq found
(cdr found
))
656 (not (equal fac
(nth 1 (car found
))))))
659 (setcar (cdr (cdr (car found
))) (+ pow
(nth 2 (car found
))))
661 ;; Put constant term first.
662 (if (and (cdr facs
) (Math-ratp (nth 1 (nth 1 facs
))))
663 (cons 'vec
(cons (nth 1 facs
) (cons (list 'vec fac pow
)
665 (cons 'vec
(cons (list 'vec fac pow
) (cdr facs
))))))))
666 (math-mul (math-pow fac pow
) facs
)))
668 (defun math-factor-poly-coefs (p &optional square-free
) ; uses "x"
673 ;; Strip off multiples of math-fet-x.
674 ((Math-zerop (car p
))
676 (while (and p
(Math-zerop (car p
)))
677 (setq z
(1+ z
) p
(cdr p
)))
679 (setq p
(math-factor-poly-coefs p square-free
))
680 (setq p
(math-sort-terms (math-factor-expr (car p
)))))
681 (math-accum-factors math-fet-x z
(math-factor-protect p
))))
683 ;; Factor out content.
684 ((and (not square-free
)
685 (not (eq 1 (setq t1
(math-mul (math-poly-gcd-list p
)
686 (if (math-guess-if-neg
687 (nth (1- (length p
)) p
))
689 (math-accum-factors t1
1 (math-factor-poly-coefs
690 (math-poly-div-list p t1
) 'cont
)))
692 ;; Check if linear in math-fet-x.
695 (math-add (math-factor-protect
697 (math-factor-expr (car p
))))
698 (math-mul math-fet-x
(math-factor-protect
700 (math-factor-expr (nth 1 p
))))))))
702 ;; If symbolic coefficients, use FactorRules.
704 (while (and pp
(or (Math-ratp (car pp
))
705 (and (eq (car (car pp
)) 'mod
)
706 (Math-integerp (nth 1 (car pp
)))
707 (Math-integerp (nth 2 (car pp
))))))
710 (let ((res (math-rewrite
711 (list 'calcFunc-thecoefs math-fet-x
(cons 'vec p
))
712 '(var FactorRules var-FactorRules
))))
713 (or (and (eq (car-safe res
) 'calcFunc-thefactors
)
715 (math-vectorp (nth 2 res
))
718 (while (setq vec
(cdr vec
))
719 (setq facs
(math-accum-factors (car vec
) 1 facs
)))
721 (math-build-polynomial-expr p math-fet-x
))))
723 ;; Check if rational coefficients (i.e., not modulo a prime).
724 ((eq math-poly-modulus
1)
726 ;; Check if there are any squared terms, or a content not = 1.
727 (if (or (eq square-free t
)
728 (equal (setq t1
(math-poly-gcd-coefs
729 p
(setq t2
(math-poly-deriv-coefs p
))))
732 ;; We now have a square-free polynomial with integer coefs.
733 ;; For now, we use a kludgey method that finds linear and
734 ;; quadratic terms using floating-point root-finding.
735 (if (setq t1
(let ((calc-symbolic-mode nil
))
736 (math-poly-all-roots nil p t
)))
737 (let ((roots (car t1
))
738 (csign (if (math-negp (nth (1- (length p
)) p
)) -
1 1))
743 (let ((coef0 (car (car roots
)))
744 (coef1 (cdr (car roots
))))
745 (setq expr
(math-accum-factors
747 (let ((den (math-lcm-denoms
749 (setq scale
(math-div scale den
))
752 (math-mul den
(math-pow math-fet-x
2))
753 (math-mul (math-mul coef1 den
)
755 (math-mul coef0 den
)))
756 (let ((den (math-lcm-denoms coef0
)))
757 (setq scale
(math-div scale den
))
758 (math-add (math-mul den math-fet-x
)
759 (math-mul coef0 den
))))
762 (setq expr
(math-accum-factors
765 (math-build-polynomial-expr
766 (math-mul-list (nth 1 t1
) scale
)
768 (math-build-polynomial-expr p math-fet-x
)) ; can't factor it.
770 ;; Separate out the squared terms (Knuth exercise 4.6.2-34).
771 ;; This step also divides out the content of the polynomial.
772 (let* ((cabs (math-poly-gcd-list p
))
773 (csign (if (math-negp (nth (1- (length p
)) p
)) -
1 1))
774 (t1s (math-mul-list t1 csign
))
776 (v (car (math-poly-div-coefs p t1s
)))
777 (w (car (math-poly-div-coefs t2 t1s
))))
779 (not (math-poly-zerop
780 (setq t2
(math-poly-simplify
782 w
1 (math-poly-deriv-coefs v
) -
1)))))
783 (setq t1
(math-poly-gcd-coefs v t2
)
785 v
(car (math-poly-div-coefs v t1
))
786 w
(car (math-poly-div-coefs t2 t1
))))
788 t2
(math-accum-factors (math-factor-poly-coefs v t
)
791 (setq t2
(math-accum-factors (math-factor-poly-coefs
796 (math-accum-factors (math-mul cabs csign
) 1 t2
))))
798 ;; Factoring modulo a prime.
799 ((and (= (length (setq temp
(math-poly-gcd-coefs
800 p
(math-poly-deriv-coefs p
))))
804 (setq temp
(nthcdr (nth 2 math-poly-modulus
) temp
)
805 p
(cons (car temp
) p
)))
806 (and (setq temp
(math-factor-poly-coefs p
))
807 (math-pow temp
(nth 2 math-poly-modulus
))))
809 (math-reject-arg nil
"*Modulo factorization not yet implemented")))))
811 (defun math-poly-deriv-coefs (p)
814 (while (setq p
(cdr p
))
815 (setq dp
(cons (math-mul (car p
) n
) dp
)
819 (defun math-factor-contains (x a
)
822 (if (memq (car-safe x
) '(+ -
* / neg
))
824 (while (setq x
(cdr x
))
825 (setq sum
(+ sum
(math-factor-contains (car x
) a
))))
827 (if (and (eq (car-safe x
) '^
)
829 (* (math-factor-contains (nth 1 x
) a
) (nth 2 x
))
836 ;;; Merge all quotients and expand/simplify the numerator
837 (defun calcFunc-nrat (expr)
838 (if (math-any-floats expr
)
839 (setq expr
(calcFunc-pfrac expr
)))
840 (if (or (math-vectorp expr
)
841 (assq (car-safe expr
) calc-tweak-eqn-table
))
842 (cons (car expr
) (mapcar 'calcFunc-nrat
(cdr expr
)))
843 (let* ((calc-prefer-frac t
)
844 (res (math-to-ratpoly expr
))
845 (num (math-simplify (math-sort-terms (calcFunc-expand (car res
)))))
846 (den (math-simplify (math-sort-terms (calcFunc-expand (cdr res
)))))
847 (g (math-poly-gcd num den
)))
849 (let ((num2 (math-poly-div num g
))
850 (den2 (math-poly-div den g
)))
851 (and (eq (cdr num2
) 0) (eq (cdr den2
) 0)
852 (setq num
(car num2
) den
(car den2
)))))
853 (math-simplify (math-div num den
)))))
855 ;;; Returns expressions (num . denom).
856 (defun math-to-ratpoly (expr)
857 (let ((res (math-to-ratpoly-rec expr
)))
858 (cons (math-simplify (car res
)) (math-simplify (cdr res
)))))
860 (defun math-to-ratpoly-rec (expr)
861 (cond ((Math-primp expr
)
863 ((memq (car expr
) '(+ -
))
864 (let ((r1 (math-to-ratpoly-rec (nth 1 expr
)))
865 (r2 (math-to-ratpoly-rec (nth 2 expr
))))
866 (if (equal (cdr r1
) (cdr r2
))
867 (cons (list (car expr
) (car r1
) (car r2
)) (cdr r1
))
869 (cons (list (car expr
)
870 (math-mul (car r1
) (cdr r2
))
874 (cons (list (car expr
)
876 (math-mul (car r2
) (cdr r1
)))
878 (let ((g (math-poly-gcd (cdr r1
) (cdr r2
))))
879 (let ((d1 (and (not (eq g
1)) (math-poly-div (cdr r1
) g
)))
880 (d2 (and (not (eq g
1)) (math-poly-div
881 (math-mul (car r1
) (cdr r2
))
883 (if (and (eq (cdr d1
) 0) (eq (cdr d2
) 0))
884 (cons (list (car expr
) (car d2
)
885 (math-mul (car r2
) (car d1
)))
886 (math-mul (car d1
) (cdr r2
)))
887 (cons (list (car expr
)
888 (math-mul (car r1
) (cdr r2
))
889 (math-mul (car r2
) (cdr r1
)))
890 (math-mul (cdr r1
) (cdr r2
)))))))))))
892 (let* ((r1 (math-to-ratpoly-rec (nth 1 expr
)))
893 (r2 (math-to-ratpoly-rec (nth 2 expr
)))
894 (g (math-mul (math-poly-gcd (car r1
) (cdr r2
))
895 (math-poly-gcd (cdr r1
) (car r2
)))))
897 (cons (math-mul (car r1
) (car r2
))
898 (math-mul (cdr r1
) (cdr r2
)))
899 (cons (math-poly-div-exact (math-mul (car r1
) (car r2
)) g
)
900 (math-poly-div-exact (math-mul (cdr r1
) (cdr r2
)) g
)))))
902 (let* ((r1 (math-to-ratpoly-rec (nth 1 expr
)))
903 (r2 (math-to-ratpoly-rec (nth 2 expr
))))
904 (if (and (eq (cdr r1
) 1) (eq (cdr r2
) 1))
905 (cons (car r1
) (car r2
))
906 (let ((g (math-mul (math-poly-gcd (car r1
) (car r2
))
907 (math-poly-gcd (cdr r1
) (cdr r2
)))))
909 (cons (math-mul (car r1
) (cdr r2
))
910 (math-mul (cdr r1
) (car r2
)))
911 (cons (math-poly-div-exact (math-mul (car r1
) (cdr r2
)) g
)
912 (math-poly-div-exact (math-mul (cdr r1
) (car r2
))
914 ((and (eq (car expr
) '^
) (integerp (nth 2 expr
)))
915 (let ((r1 (math-to-ratpoly-rec (nth 1 expr
))))
916 (if (> (nth 2 expr
) 0)
917 (cons (math-pow (car r1
) (nth 2 expr
))
918 (math-pow (cdr r1
) (nth 2 expr
)))
919 (cons (math-pow (cdr r1
) (- (nth 2 expr
)))
920 (math-pow (car r1
) (- (nth 2 expr
)))))))
921 ((eq (car expr
) 'neg
)
922 (let ((r1 (math-to-ratpoly-rec (nth 1 expr
))))
923 (cons (math-neg (car r1
)) (cdr r1
))))
927 (defun math-ratpoly-p (expr &optional var
)
928 (cond ((equal expr var
) 1)
929 ((Math-primp expr
) 0)
930 ((memq (car expr
) '(+ -
))
931 (let ((p1 (math-ratpoly-p (nth 1 expr
) var
))
933 (and p1
(setq p2
(math-ratpoly-p (nth 2 expr
) var
))
936 (let ((p1 (math-ratpoly-p (nth 1 expr
) var
))
938 (and p1
(setq p2
(math-ratpoly-p (nth 2 expr
) var
))
940 ((eq (car expr
) 'neg
)
941 (math-ratpoly-p (nth 1 expr
) var
))
943 (let ((p1 (math-ratpoly-p (nth 1 expr
) var
))
945 (and p1
(setq p2
(math-ratpoly-p (nth 2 expr
) var
))
947 ((and (eq (car expr
) '^
)
948 (integerp (nth 2 expr
)))
949 (let ((p1 (math-ratpoly-p (nth 1 expr
) var
)))
950 (and p1
(* p1
(nth 2 expr
)))))
952 ((math-poly-depends expr var
) nil
)
956 (defun calcFunc-apart (expr &optional var
)
957 (cond ((Math-primp expr
) expr
)
959 (math-add (calcFunc-apart (nth 1 expr
) var
)
960 (calcFunc-apart (nth 2 expr
) var
)))
962 (math-sub (calcFunc-apart (nth 1 expr
) var
)
963 (calcFunc-apart (nth 2 expr
) var
)))
964 ((and var
(not (math-ratpoly-p expr var
)))
965 (math-reject-arg expr
"Expected a rational function"))
967 (let* ((calc-prefer-frac t
)
968 (rat (math-to-ratpoly expr
))
972 (setq var
(math-polynomial-base den
)))
973 (if (not (math-ratpoly-p expr var
))
974 (math-reject-arg expr
"Expected a rational function")
975 (let* ((qr (math-poly-div num den
))
978 (math-add q
(or (and var
979 (math-expr-contains den var
)
980 (math-partial-fractions r den var
))
981 (math-div r den
)))))))))
984 (defun math-padded-polynomial (expr var deg
)
985 "Return a polynomial as list of coefficients.
986 If EXPR is of the form \"a + bx + cx^2 + ...\" in the variable VAR, return
987 the list (a b c ...) with at least DEG elements, else return NIL."
988 (let ((p (math-is-polynomial expr var deg
)))
989 (append p
(make-list (- deg
(length p
)) 0))))
991 (defun math-partial-fractions (r den var
)
992 "Return R divided by DEN expressed in partial fractions of VAR.
993 All whole factors of DEN have already been split off from R.
994 If no partial fraction representation can be found, return nil."
995 (let* ((fden (calcFunc-factors den var
))
996 (tdeg (math-polynomial-p den var
))
1001 (tz (make-list (1- tdeg
) 0))
1002 (calc-matrix-mode 'scalar
))
1003 (and (not (and (= (length fden
) 2) (eq (nth 2 (nth 1 fden
)) 1)))
1005 (while (setq fp
(cdr fp
))
1006 (let ((rpt (nth 2 (car fp
)))
1007 (deg (math-polynomial-p (nth 1 (car fp
)) var
))
1013 (setq dvar
(append '(vec) lz
'(1) tz
)
1017 dnum
(math-add dnum
(math-mul dvar
1018 (math-pow var deg2
)))
1019 dlist
(cons (and (= deg2
(1- deg
))
1020 (math-pow (nth 1 (car fp
)) rpt
))
1024 (while (setq fpp
(cdr fpp
))
1026 (setq mult
(math-mul mult
1027 (math-pow (nth 1 (car fpp
))
1028 (nth 2 (car fpp
)))))))
1029 (setq dnum
(math-mul dnum mult
)))
1030 (setq eqns
(math-add eqns
(math-mul dnum
1036 (setq eqns
(math-div (cons 'vec
(math-padded-polynomial r var tdeg
))
1042 (cons 'vec
(math-padded-polynomial
1045 (and (math-vectorp eqns
)
1048 (setq eqns
(nreverse eqns
))
1050 (setq num
(cons (car eqns
) num
)
1053 (setq num
(math-build-polynomial-expr
1055 res
(math-add res
(math-div num
(car dlist
)))
1057 (setq dlist
(cdr dlist
)))
1058 (math-normalize res
)))))))
1062 (defun math-expand-term (expr)
1063 (cond ((and (eq (car-safe expr
) '*)
1064 (memq (car-safe (nth 1 expr
)) '(+ -
)))
1065 (math-add-or-sub (list '* (nth 1 (nth 1 expr
)) (nth 2 expr
))
1066 (list '* (nth 2 (nth 1 expr
)) (nth 2 expr
))
1067 nil
(eq (car (nth 1 expr
)) '-
)))
1068 ((and (eq (car-safe expr
) '*)
1069 (memq (car-safe (nth 2 expr
)) '(+ -
)))
1070 (math-add-or-sub (list '* (nth 1 expr
) (nth 1 (nth 2 expr
)))
1071 (list '* (nth 1 expr
) (nth 2 (nth 2 expr
)))
1072 nil
(eq (car (nth 2 expr
)) '-
)))
1073 ((and (eq (car-safe expr
) '/)
1074 (memq (car-safe (nth 1 expr
)) '(+ -
)))
1075 (math-add-or-sub (list '/ (nth 1 (nth 1 expr
)) (nth 2 expr
))
1076 (list '/ (nth 2 (nth 1 expr
)) (nth 2 expr
))
1077 nil
(eq (car (nth 1 expr
)) '-
)))
1078 ((and (eq (car-safe expr
) '^
)
1079 (memq (car-safe (nth 1 expr
)) '(+ -
))
1080 (integerp (nth 2 expr
))
1082 (or (math-known-matrixp (nth 1 (nth 1 expr
)))
1083 (math-known-matrixp (nth 2 (nth 1 expr
)))
1086 (not (eq calc-matrix-mode
'scalar
))
1087 (not (and (math-known-scalarp (nth 1 (nth 1 expr
)))
1088 (math-known-scalarp (nth 2 (nth 1 expr
)))))))
1090 (if (= (nth 2 expr
) 2)
1091 (math-add-or-sub (list '* (nth 1 (nth 1 expr
)) (nth 1 expr
))
1092 (list '* (nth 2 (nth 1 expr
)) (nth 1 expr
))
1093 nil
(eq (car (nth 1 expr
)) '-
))
1094 (math-add-or-sub (list '* (nth 1 (nth 1 expr
))
1095 (list '^
(nth 1 expr
)
1097 (list '* (nth 2 (nth 1 expr
))
1098 (list '^
(nth 1 expr
)
1100 nil
(eq (car (nth 1 expr
)) '-
)))
1101 (if (> (nth 2 expr
) 0)
1102 (or (and (or (> math-mt-many
500000) (< math-mt-many -
500000))
1103 (math-expand-power (nth 1 expr
) (nth 2 expr
)
1107 (list '^
(nth 1 expr
) (1- (nth 2 expr
)))))
1108 (if (< (nth 2 expr
) 0)
1109 (list '/ 1 (list '^
(nth 1 expr
) (- (nth 2 expr
)))))))))
1112 (defun calcFunc-expand (expr &optional many
)
1113 (math-normalize (math-map-tree 'math-expand-term expr many
)))
1115 (defun math-expand-power (x n
&optional var else-nil
)
1116 (or (and (natnump n
)
1117 (memq (car-safe x
) '(+ -
))
1120 (while (memq (car-safe x
) '(+ -
))
1121 (setq terms
(cons (if (eq (car x
) '-
)
1122 (math-neg (nth 2 x
))
1126 (setq terms
(cons x terms
))
1130 (or (math-expr-contains (car p
) var
)
1131 (setq terms
(delq (car p
) terms
)
1132 cterms
(cons (car p
) cterms
)))
1135 (setq terms
(cons (apply 'calcFunc-add cterms
)
1137 (if (= (length terms
) 2)
1141 (setq accum
(list '+ accum
1142 (list '* (calcFunc-choose n i
)
1144 (list '^
(nth 1 terms
) i
)
1145 (list '^
(car terms
)
1154 (setq accum
(list '+ accum
1155 (list '^
(car p1
) 2))
1157 (while (setq p2
(cdr p2
))
1158 (setq accum
(list '+ accum
1169 (setq accum
(list '+ accum
(list '^
(car p1
) 3))
1171 (while (setq p2
(cdr p2
))
1172 (setq accum
(list '+
1178 (list '^
(car p1
) 2)
1183 (list '^
(car p2
) 2))))
1185 (while (setq p3
(cdr p3
))
1186 (setq accum
(list '+ accum
1198 (defun calcFunc-expandpow (x n
)
1199 (math-normalize (math-expand-power x n
)))
1201 (provide 'calc-poly
)
1203 ;; arch-tag: d2566c51-2ccc-45f1-8c50-f3462c2953ff
1204 ;;; calc-poly.el ends here