1 ;;; calc-poly.el --- polynomial functions for Calc
3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc.
5 ;; Author: David Gillespie <daveg@synaptics.com>
6 ;; Maintainer: Jay Belanger <belanger@truman.edu>
8 ;; This file is part of GNU Emacs.
10 ;; GNU Emacs is distributed in the hope that it will be useful,
11 ;; but WITHOUT ANY WARRANTY. No author or distributor
12 ;; accepts responsibility to anyone for the consequences of using it
13 ;; or for whether it serves any particular purpose or works at all,
14 ;; unless he says so in writing. Refer to the GNU Emacs General Public
15 ;; License for full details.
17 ;; Everyone is granted permission to copy, modify and redistribute
18 ;; GNU Emacs, but only under the conditions described in the
19 ;; GNU Emacs General Public License. A copy of this license is
20 ;; supposed to have been given to you along with GNU Emacs so you
21 ;; can know your rights and responsibilities. It should be in a
22 ;; file named COPYING. Among other things, the copyright notice
23 ;; and this notice must be preserved on all copies.
29 ;; This file is autoloaded from calc-ext.el.
34 (defun calc-Need-calc-poly () nil
)
37 (defun calcFunc-pcont (expr &optional var
)
38 (cond ((Math-primp expr
)
39 (cond ((Math-zerop expr
) 1)
40 ((Math-messy-integerp expr
) (math-trunc expr
))
41 ((Math-objectp expr
) expr
)
42 ((or (equal expr var
) (not var
)) 1)
45 (math-mul (calcFunc-pcont (nth 1 expr
) var
)
46 (calcFunc-pcont (nth 2 expr
) var
)))
48 (math-div (calcFunc-pcont (nth 1 expr
) var
)
49 (calcFunc-pcont (nth 2 expr
) var
)))
50 ((and (eq (car expr
) '^
) (Math-natnump (nth 2 expr
)))
51 (math-pow (calcFunc-pcont (nth 1 expr
) var
) (nth 2 expr
)))
52 ((memq (car expr
) '(neg polar
))
53 (calcFunc-pcont (nth 1 expr
) var
))
55 (let ((p (math-is-polynomial expr var
)))
57 (let ((lead (nth (1- (length p
)) p
))
58 (cont (math-poly-gcd-list p
)))
59 (if (math-guess-if-neg lead
)
63 ((memq (car expr
) '(+ - cplx sdev
))
64 (let ((cont (calcFunc-pcont (nth 1 expr
) var
)))
67 (let ((c2 (calcFunc-pcont (nth 2 expr
) var
)))
68 (if (and (math-negp cont
)
69 (if (eq (car expr
) '-
) (math-posp c2
) (math-negp c2
)))
70 (math-neg (math-poly-gcd cont c2
))
71 (math-poly-gcd cont c2
))))))
75 (defun calcFunc-pprim (expr &optional var
)
76 (let ((cont (calcFunc-pcont expr var
)))
77 (if (math-equal-int cont
1)
79 (math-poly-div-exact expr cont var
))))
81 (defun math-div-poly-const (expr c
)
82 (cond ((memq (car-safe expr
) '(+ -
))
84 (math-div-poly-const (nth 1 expr
) c
)
85 (math-div-poly-const (nth 2 expr
) c
)))
86 (t (math-div expr c
))))
88 (defun calcFunc-pdeg (expr &optional var
)
90 '(neg (var inf var-inf
))
92 (or (math-polynomial-p expr var
)
93 (math-reject-arg expr
"Expected a polynomial"))
94 (math-poly-degree expr
))))
96 (defun math-poly-degree (expr)
97 (cond ((Math-primp expr
)
98 (if (eq (car-safe expr
) 'var
) 1 0))
100 (math-poly-degree (nth 1 expr
)))
102 (+ (math-poly-degree (nth 1 expr
))
103 (math-poly-degree (nth 2 expr
))))
105 (- (math-poly-degree (nth 1 expr
))
106 (math-poly-degree (nth 2 expr
))))
107 ((and (eq (car expr
) '^
) (natnump (nth 2 expr
)))
108 (* (math-poly-degree (nth 1 expr
)) (nth 2 expr
)))
109 ((memq (car expr
) '(+ -
))
110 (max (math-poly-degree (nth 1 expr
))
111 (math-poly-degree (nth 2 expr
))))
114 (defun calcFunc-plead (expr var
)
115 (cond ((eq (car-safe expr
) '*)
116 (math-mul (calcFunc-plead (nth 1 expr
) var
)
117 (calcFunc-plead (nth 2 expr
) var
)))
118 ((eq (car-safe expr
) '/)
119 (math-div (calcFunc-plead (nth 1 expr
) var
)
120 (calcFunc-plead (nth 2 expr
) var
)))
121 ((and (eq (car-safe expr
) '^
) (math-natnump (nth 2 expr
)))
122 (math-pow (calcFunc-plead (nth 1 expr
) var
) (nth 2 expr
)))
128 (let ((p (math-is-polynomial expr var
)))
130 (nth (1- (length p
)) p
)
137 ;;; Polynomial quotient, remainder, and GCD.
138 ;;; Originally by Ove Ewerlid (ewerlid@mizar.DoCS.UU.SE).
139 ;;; Modifications and simplifications by daveg.
141 (defvar math-poly-modulus
1)
143 ;;; Return gcd of two polynomials
144 (defun calcFunc-pgcd (pn pd
)
145 (if (math-any-floats pn
)
146 (math-reject-arg pn
"Coefficients must be rational"))
147 (if (math-any-floats pd
)
148 (math-reject-arg pd
"Coefficients must be rational"))
149 (let ((calc-prefer-frac t
)
150 (math-poly-modulus (math-poly-modulus pn pd
)))
151 (math-poly-gcd pn pd
)))
153 ;;; Return only quotient to top of stack (nil if zero)
155 ;; calc-poly-div-remainder is a local variable for
156 ;; calc-poly-div (in calc-alg.el), but is used by
157 ;; calcFunc-pdiv, which is called by calc-poly-div.
158 (defvar calc-poly-div-remainder
)
160 (defun calcFunc-pdiv (pn pd
&optional base
)
161 (let* ((calc-prefer-frac t
)
162 (math-poly-modulus (math-poly-modulus pn pd
))
163 (res (math-poly-div pn pd base
)))
164 (setq calc-poly-div-remainder
(cdr res
))
167 ;;; Return only remainder to top of stack
168 (defun calcFunc-prem (pn pd
&optional base
)
169 (let ((calc-prefer-frac t
)
170 (math-poly-modulus (math-poly-modulus pn pd
)))
171 (cdr (math-poly-div pn pd base
))))
173 (defun calcFunc-pdivrem (pn pd
&optional base
)
174 (let* ((calc-prefer-frac t
)
175 (math-poly-modulus (math-poly-modulus pn pd
))
176 (res (math-poly-div pn pd base
)))
177 (list 'vec
(car res
) (cdr res
))))
179 (defun calcFunc-pdivide (pn pd
&optional base
)
180 (let* ((calc-prefer-frac t
)
181 (math-poly-modulus (math-poly-modulus pn pd
))
182 (res (math-poly-div pn pd base
)))
183 (math-add (car res
) (math-div (cdr res
) pd
))))
186 ;;; Multiply two terms, expanding out products of sums.
187 (defun math-mul-thru (lhs rhs
)
188 (if (memq (car-safe lhs
) '(+ -
))
190 (math-mul-thru (nth 1 lhs
) rhs
)
191 (math-mul-thru (nth 2 lhs
) rhs
))
192 (if (memq (car-safe rhs
) '(+ -
))
194 (math-mul-thru lhs
(nth 1 rhs
))
195 (math-mul-thru lhs
(nth 2 rhs
)))
196 (math-mul lhs rhs
))))
198 (defun math-div-thru (num den
)
199 (if (memq (car-safe num
) '(+ -
))
201 (math-div-thru (nth 1 num
) den
)
202 (math-div-thru (nth 2 num
) den
))
206 ;;; Sort the terms of a sum into canonical order.
207 (defun math-sort-terms (expr)
208 (if (memq (car-safe expr
) '(+ -
))
210 (sort (math-sum-to-list expr
)
211 (function (lambda (a b
) (math-beforep (car a
) (car b
))))))
214 (defun math-list-to-sum (lst)
216 (list (if (cdr (car lst
)) '-
'+)
217 (math-list-to-sum (cdr lst
))
220 (math-neg (car (car lst
)))
223 (defun math-sum-to-list (tree &optional neg
)
224 (cond ((eq (car-safe tree
) '+)
225 (nconc (math-sum-to-list (nth 1 tree
) neg
)
226 (math-sum-to-list (nth 2 tree
) neg
)))
227 ((eq (car-safe tree
) '-
)
228 (nconc (math-sum-to-list (nth 1 tree
) neg
)
229 (math-sum-to-list (nth 2 tree
) (not neg
))))
230 (t (list (cons tree neg
)))))
232 ;;; Check if the polynomial coefficients are modulo forms.
233 (defun math-poly-modulus (expr &optional expr2
)
234 (or (math-poly-modulus-rec expr
)
235 (and expr2
(math-poly-modulus-rec expr2
))
238 (defun math-poly-modulus-rec (expr)
239 (if (and (eq (car-safe expr
) 'mod
) (Math-natnump (nth 2 expr
)))
240 (list 'mod
1 (nth 2 expr
))
241 (and (memq (car-safe expr
) '(+ -
* /))
242 (or (math-poly-modulus-rec (nth 1 expr
))
243 (math-poly-modulus-rec (nth 2 expr
))))))
246 ;;; Divide two polynomials. Return (quotient . remainder).
247 (defvar math-poly-div-base nil
)
248 (defun math-poly-div (u v
&optional math-poly-div-base
)
249 (if math-poly-div-base
250 (math-do-poly-div u v
)
251 (math-do-poly-div (calcFunc-expand u
) (calcFunc-expand v
))))
253 (defun math-poly-div-exact (u v
&optional base
)
254 (let ((res (math-poly-div u v base
)))
257 (math-reject-arg (list 'vec u v
) "Argument is not a polynomial"))))
259 (defun math-do-poly-div (u v
)
260 (cond ((math-constp u
)
262 (cons (math-div u v
) 0)
267 (if (memq (car-safe u
) '(+ -
))
268 (math-add-or-sub (math-poly-div-exact (nth 1 u
) v
)
269 (math-poly-div-exact (nth 2 u
) v
)
274 (cons math-poly-modulus
0))
275 ((and (math-atomic-factorp u
) (math-atomic-factorp v
))
276 (cons (math-simplify (math-div u v
)) 0))
278 (let ((base (or math-poly-div-base
279 (math-poly-div-base u v
)))
282 (null (setq vp
(math-is-polynomial v base nil
'gen
))))
284 (setq up
(math-is-polynomial u base nil
'gen
)
285 res
(math-poly-div-coefs up vp
))
286 (cons (math-build-polynomial-expr (car res
) base
)
287 (math-build-polynomial-expr (cdr res
) base
)))))))
289 (defun math-poly-div-rec (u v
)
290 (cond ((math-constp u
)
295 (if (memq (car-safe u
) '(+ -
))
296 (math-add-or-sub (math-poly-div-rec (nth 1 u
) v
)
297 (math-poly-div-rec (nth 2 u
) v
)
300 ((Math-equal u v
) math-poly-modulus
)
301 ((and (math-atomic-factorp u
) (math-atomic-factorp v
))
302 (math-simplify (math-div u v
)))
306 (let ((base (math-poly-div-base u v
))
309 (null (setq vp
(math-is-polynomial v base nil
'gen
))))
311 (setq up
(math-is-polynomial u base nil
'gen
)
312 res
(math-poly-div-coefs up vp
))
313 (math-add (math-build-polynomial-expr (car res
) base
)
314 (math-div (math-build-polynomial-expr (cdr res
) base
)
317 ;;; Divide two polynomials in coefficient-list form. Return (quot . rem).
318 (defun math-poly-div-coefs (u v
)
319 (cond ((null v
) (math-reject-arg nil
"Division by zero"))
320 ((< (length u
) (length v
)) (cons nil u
))
326 (let ((qk (math-poly-div-rec (math-simplify (car urev
))
330 (if (or q
(not (math-zerop qk
)))
331 (setq q
(cons qk q
)))
332 (while (setq up
(cdr up
) vp
(cdr vp
))
333 (setcar up
(math-sub (car up
) (math-mul-thru qk
(car vp
)))))
334 (setq urev
(cdr urev
))
336 (while (and urev
(Math-zerop (car urev
)))
337 (setq urev
(cdr urev
)))
338 (cons q
(nreverse (mapcar 'math-simplify urev
)))))
340 (cons (list (math-poly-div-rec (car u
) (car v
)))
343 ;;; Perform a pseudo-division of polynomials. (See Knuth section 4.6.1.)
344 ;;; This returns only the remainder from the pseudo-division.
345 (defun math-poly-pseudo-div (u v
)
347 ((< (length u
) (length v
)) u
)
348 ((or (cdr u
) (cdr v
))
349 (let ((urev (reverse u
))
355 (while (setq up
(cdr up
) vp
(cdr vp
))
356 (setcar up
(math-sub (math-mul-thru (car vrev
) (car up
))
357 (math-mul-thru (car urev
) (car vp
)))))
358 (setq urev
(cdr urev
))
361 (setcar up
(math-mul-thru (car vrev
) (car up
)))
363 (while (and urev
(Math-zerop (car urev
)))
364 (setq urev
(cdr urev
)))
365 (nreverse (mapcar 'math-simplify urev
))))
368 ;;; Compute the GCD of two multivariate polynomials.
369 (defun math-poly-gcd (u v
)
370 (cond ((Math-equal u v
) u
)
374 (calcFunc-gcd u
(calcFunc-pcont v
))))
378 (calcFunc-gcd v
(calcFunc-pcont u
))))
380 (let ((base (math-poly-gcd-base u v
)))
384 (math-build-polynomial-expr
385 (math-poly-gcd-coefs (math-is-polynomial u base nil
'gen
)
386 (math-is-polynomial v base nil
'gen
))
388 (calcFunc-gcd (calcFunc-pcont u
) (calcFunc-pcont u
)))))))
390 (defun math-poly-div-list (lst a
)
394 (math-mul-list lst a
)
395 (mapcar (function (lambda (x) (math-poly-div-exact x a
))) lst
))))
397 (defun math-mul-list (lst a
)
401 (mapcar 'math-neg lst
)
403 (mapcar (function (lambda (x) (math-mul x a
))) lst
)))))
405 ;;; Run GCD on all elements in a list.
406 (defun math-poly-gcd-list (lst)
407 (if (or (memq 1 lst
) (memq -
1 lst
))
408 (math-poly-gcd-frac-list lst
)
409 (let ((gcd (car lst
)))
410 (while (and (setq lst
(cdr lst
)) (not (eq gcd
1)))
412 (setq gcd
(math-poly-gcd gcd
(car lst
)))))
413 (if lst
(setq lst
(math-poly-gcd-frac-list lst
)))
416 (defun math-poly-gcd-frac-list (lst)
417 (while (and lst
(not (eq (car-safe (car lst
)) 'frac
)))
418 (setq lst
(cdr lst
)))
420 (let ((denom (nth 2 (car lst
))))
421 (while (setq lst
(cdr lst
))
422 (if (eq (car-safe (car lst
)) 'frac
)
423 (setq denom
(calcFunc-lcm denom
(nth 2 (car lst
))))))
424 (list 'frac
1 denom
))
427 ;;; Compute the GCD of two monovariate polynomial lists.
428 ;;; Knuth section 4.6.1, algorithm C.
429 (defun math-poly-gcd-coefs (u v
)
430 (let ((d (math-poly-gcd (math-poly-gcd-list u
)
431 (math-poly-gcd-list v
)))
432 (g 1) (h 1) (z 0) hh r delta ghd
)
433 (while (and u v
(Math-zerop (car u
)) (Math-zerop (car v
)))
434 (setq u
(cdr u
) v
(cdr v
) z
(1+ z
)))
436 (setq u
(math-poly-div-list u d
)
437 v
(math-poly-div-list v d
)))
439 (setq delta
(- (length u
) (length v
)))
441 (setq r u u v v r delta
(- delta
)))
442 (setq r
(math-poly-pseudo-div u v
))
445 v
(math-poly-div-list r
(math-mul g
(math-pow h delta
)))
446 g
(nth (1- (length u
)) u
)
448 (math-mul (math-pow g delta
) (math-pow h
(- 1 delta
)))
449 (math-poly-div-exact (math-pow g delta
)
450 (math-pow h
(1- delta
))))))
453 (math-mul-list (math-poly-div-list v
(math-poly-gcd-list v
)) d
)))
454 (if (math-guess-if-neg (nth (1- (length v
)) v
))
455 (setq v
(math-mul-list v -
1)))
456 (while (>= (setq z
(1- z
)) 0)
461 ;;; Return true if is a factor containing no sums or quotients.
462 (defun math-atomic-factorp (expr)
463 (cond ((eq (car-safe expr
) '*)
464 (and (math-atomic-factorp (nth 1 expr
))
465 (math-atomic-factorp (nth 2 expr
))))
466 ((memq (car-safe expr
) '(+ -
/))
468 ((memq (car-safe expr
) '(^ neg
))
469 (math-atomic-factorp (nth 1 expr
)))
472 ;;; Find a suitable base for dividing a by b.
473 ;;; The base must exist in both expressions.
474 ;;; The degree in the numerator must be higher or equal than the
475 ;;; degree in the denominator.
476 ;;; If the above conditions are not met the quotient is just a remainder.
477 ;;; Return nil if this is the case.
479 (defun math-poly-div-base (a b
)
481 (and (setq a-base
(math-total-polynomial-base a
))
482 (setq b-base
(math-total-polynomial-base b
))
485 (let ((maybe (assoc (car (car a-base
)) b-base
)))
487 (if (>= (nth 1 (car a-base
)) (nth 1 maybe
))
488 (throw 'return
(car (car a-base
))))))
489 (setq a-base
(cdr a-base
)))))))
491 ;;; Same as above but for gcd algorithm.
492 ;;; Here there is no requirement that degree(a) > degree(b).
493 ;;; Take the base that has the highest degree considering both a and b.
494 ;;; ("a^20+b^21+x^3+a+b", "a+b^2+x^5+a^22+b^10") --> (a 22)
496 (defun math-poly-gcd-base (a b
)
498 (and (setq a-base
(math-total-polynomial-base a
))
499 (setq b-base
(math-total-polynomial-base b
))
501 (while (and a-base b-base
)
502 (if (> (nth 1 (car a-base
)) (nth 1 (car b-base
)))
503 (if (assoc (car (car a-base
)) b-base
)
504 (throw 'return
(car (car a-base
)))
505 (setq a-base
(cdr a-base
)))
506 (if (assoc (car (car b-base
)) a-base
)
507 (throw 'return
(car (car b-base
)))
508 (setq b-base
(cdr b-base
)))))))))
510 ;;; Sort a list of polynomial bases.
511 (defun math-sort-poly-base-list (lst)
512 (sort lst
(function (lambda (a b
)
513 (or (> (nth 1 a
) (nth 1 b
))
514 (and (= (nth 1 a
) (nth 1 b
))
515 (math-beforep (car a
) (car b
))))))))
517 ;;; Given an expression find all variables that are polynomial bases.
518 ;;; Return list in the form '( (var1 degree1) (var2 degree2) ... ).
519 ;;; Note dynamic scope of mpb-total-base.
520 (defun math-total-polynomial-base (expr)
521 (let ((mpb-total-base nil
))
522 (math-polynomial-base expr
'math-polynomial-p1
)
523 (math-sort-poly-base-list mpb-total-base
)))
525 (defun math-polynomial-p1 (subexpr)
526 (or (assoc subexpr mpb-total-base
)
527 (memq (car subexpr
) '(+ -
* / neg
))
528 (and (eq (car subexpr
) '^
) (natnump (nth 2 subexpr
)))
529 (let* ((math-poly-base-variable subexpr
)
530 (exponent (math-polynomial-p mpb-top-expr subexpr
)))
532 (setq mpb-total-base
(cons (list subexpr exponent
)
539 (defun calcFunc-factors (expr &optional var
)
540 (let ((math-factored-vars (if var t nil
))
542 (calc-prefer-frac t
))
544 (setq var
(math-polynomial-base expr
)))
545 (let ((res (math-factor-finish
546 (or (catch 'factor
(math-factor-expr-try var
))
548 (math-simplify (if (math-vectorp res
)
550 (list 'vec
(list 'vec res
1)))))))
552 (defun calcFunc-factor (expr &optional var
)
553 (let ((math-factored-vars nil
)
555 (calc-prefer-frac t
))
556 (math-simplify (math-factor-finish
558 (let ((math-factored-vars t
))
559 (or (catch 'factor
(math-factor-expr-try var
)) expr
))
560 (math-factor-expr expr
))))))
562 (defun math-factor-finish (x)
565 (if (eq (car x
) 'calcFunc-Fac-Prot
)
566 (math-factor-finish (nth 1 x
))
567 (cons (car x
) (mapcar 'math-factor-finish
(cdr x
))))))
569 (defun math-factor-protect (x)
570 (if (memq (car-safe x
) '(+ -
))
571 (list 'calcFunc-Fac-Prot x
)
574 (defun math-factor-expr (expr)
575 (cond ((eq math-factored-vars t
) expr
)
576 ((or (memq (car-safe expr
) '(* / ^ neg
))
577 (assq (car-safe expr
) calc-tweak-eqn-table
))
578 (cons (car expr
) (mapcar 'math-factor-expr
(cdr expr
))))
579 ((memq (car-safe expr
) '(+ -
))
580 (let* ((math-factored-vars math-factored-vars
)
581 (y (catch 'factor
(math-factor-expr-part expr
))))
587 (defun math-factor-expr-part (x) ; uses "expr"
588 (if (memq (car-safe x
) '(+ -
* / ^ neg
))
589 (while (setq x
(cdr x
))
590 (math-factor-expr-part (car x
)))
591 (and (not (Math-objvecp x
))
592 (not (assoc x math-factored-vars
))
593 (> (math-factor-contains expr x
) 1)
594 (setq math-factored-vars
(cons (list x
) math-factored-vars
))
595 (math-factor-expr-try x
))))
597 (defun math-factor-expr-try (x)
598 (if (eq (car-safe expr
) '*)
599 (let ((res1 (catch 'factor
(let ((expr (nth 1 expr
)))
600 (math-factor-expr-try x
))))
601 (res2 (catch 'factor
(let ((expr (nth 2 expr
)))
602 (math-factor-expr-try x
)))))
604 (throw 'factor
(math-accum-factors (or res1
(nth 1 expr
)) 1
605 (or res2
(nth 2 expr
))))))
606 (let* ((p (math-is-polynomial expr x
30 'gen
))
607 (math-poly-modulus (math-poly-modulus expr
))
610 (setq res
(math-factor-poly-coefs p
))
611 (throw 'factor res
)))))
613 (defun math-accum-factors (fac pow facs
)
615 (if (math-vectorp fac
)
617 (while (setq fac
(cdr fac
))
618 (setq facs
(math-accum-factors (nth 1 (car fac
))
619 (* pow
(nth 2 (car fac
)))
622 (if (and (eq (car-safe fac
) '^
) (natnump (nth 2 fac
)))
623 (setq pow
(* pow
(nth 2 fac
))
627 (or (math-vectorp facs
)
628 (setq facs
(if (eq facs
1) '(vec)
629 (list 'vec
(list 'vec facs
1)))))
631 (while (and (setq found
(cdr found
))
632 (not (equal fac
(nth 1 (car found
))))))
635 (setcar (cdr (cdr (car found
))) (+ pow
(nth 2 (car found
))))
637 ;; Put constant term first.
638 (if (and (cdr facs
) (Math-ratp (nth 1 (nth 1 facs
))))
639 (cons 'vec
(cons (nth 1 facs
) (cons (list 'vec fac pow
)
641 (cons 'vec
(cons (list 'vec fac pow
) (cdr facs
))))))))
642 (math-mul (math-pow fac pow
) facs
)))
644 (defun math-factor-poly-coefs (p &optional square-free
) ; uses "x"
649 ;; Strip off multiples of x.
650 ((Math-zerop (car p
))
652 (while (and p
(Math-zerop (car p
)))
653 (setq z
(1+ z
) p
(cdr p
)))
655 (setq p
(math-factor-poly-coefs p square-free
))
656 (setq p
(math-sort-terms (math-factor-expr (car p
)))))
657 (math-accum-factors x z
(math-factor-protect p
))))
659 ;; Factor out content.
660 ((and (not square-free
)
661 (not (eq 1 (setq t1
(math-mul (math-poly-gcd-list p
)
662 (if (math-guess-if-neg
663 (nth (1- (length p
)) p
))
665 (math-accum-factors t1
1 (math-factor-poly-coefs
666 (math-poly-div-list p t1
) 'cont
)))
668 ;; Check if linear in x.
670 (math-add (math-factor-protect
672 (math-factor-expr (car p
))))
673 (math-mul x
(math-factor-protect
675 (math-factor-expr (nth 1 p
)))))))
677 ;; If symbolic coefficients, use FactorRules.
679 (while (and pp
(or (Math-ratp (car pp
))
680 (and (eq (car (car pp
)) 'mod
)
681 (Math-integerp (nth 1 (car pp
)))
682 (Math-integerp (nth 2 (car pp
))))))
685 (let ((res (math-rewrite
686 (list 'calcFunc-thecoefs x
(cons 'vec p
))
687 '(var FactorRules var-FactorRules
))))
688 (or (and (eq (car-safe res
) 'calcFunc-thefactors
)
690 (math-vectorp (nth 2 res
))
693 (while (setq vec
(cdr vec
))
694 (setq facs
(math-accum-factors (car vec
) 1 facs
)))
696 (math-build-polynomial-expr p x
))))
698 ;; Check if rational coefficients (i.e., not modulo a prime).
699 ((eq math-poly-modulus
1)
701 ;; Check if there are any squared terms, or a content not = 1.
702 (if (or (eq square-free t
)
703 (equal (setq t1
(math-poly-gcd-coefs
704 p
(setq t2
(math-poly-deriv-coefs p
))))
707 ;; We now have a square-free polynomial with integer coefs.
708 ;; For now, we use a kludgey method that finds linear and
709 ;; quadratic terms using floating-point root-finding.
710 (if (setq t1
(let ((calc-symbolic-mode nil
))
711 (math-poly-all-roots nil p t
)))
712 (let ((roots (car t1
))
713 (csign (if (math-negp (nth (1- (length p
)) p
)) -
1 1))
718 (let ((coef0 (car (car roots
)))
719 (coef1 (cdr (car roots
))))
720 (setq expr
(math-accum-factors
722 (let ((den (math-lcm-denoms
724 (setq scale
(math-div scale den
))
727 (math-mul den
(math-pow x
2))
728 (math-mul (math-mul coef1 den
) x
))
729 (math-mul coef0 den
)))
730 (let ((den (math-lcm-denoms coef0
)))
731 (setq scale
(math-div scale den
))
732 (math-add (math-mul den x
)
733 (math-mul coef0 den
))))
736 (setq expr
(math-accum-factors
739 (math-build-polynomial-expr
740 (math-mul-list (nth 1 t1
) scale
)
742 (math-build-polynomial-expr p x
)) ; can't factor it.
744 ;; Separate out the squared terms (Knuth exercise 4.6.2-34).
745 ;; This step also divides out the content of the polynomial.
746 (let* ((cabs (math-poly-gcd-list p
))
747 (csign (if (math-negp (nth (1- (length p
)) p
)) -
1 1))
748 (t1s (math-mul-list t1 csign
))
750 (v (car (math-poly-div-coefs p t1s
)))
751 (w (car (math-poly-div-coefs t2 t1s
))))
753 (not (math-poly-zerop
754 (setq t2
(math-poly-simplify
756 w
1 (math-poly-deriv-coefs v
) -
1)))))
757 (setq t1
(math-poly-gcd-coefs v t2
)
759 v
(car (math-poly-div-coefs v t1
))
760 w
(car (math-poly-div-coefs t2 t1
))))
762 t2
(math-accum-factors (math-factor-poly-coefs v t
)
765 (setq t2
(math-accum-factors (math-factor-poly-coefs
770 (math-accum-factors (math-mul cabs csign
) 1 t2
))))
772 ;; Factoring modulo a prime.
773 ((and (= (length (setq temp
(math-poly-gcd-coefs
774 p
(math-poly-deriv-coefs p
))))
778 (setq temp
(nthcdr (nth 2 math-poly-modulus
) temp
)
779 p
(cons (car temp
) p
)))
780 (and (setq temp
(math-factor-poly-coefs p
))
781 (math-pow temp
(nth 2 math-poly-modulus
))))
783 (math-reject-arg nil
"*Modulo factorization not yet implemented")))))
785 (defun math-poly-deriv-coefs (p)
788 (while (setq p
(cdr p
))
789 (setq dp
(cons (math-mul (car p
) n
) dp
)
793 (defun math-factor-contains (x a
)
796 (if (memq (car-safe x
) '(+ -
* / neg
))
798 (while (setq x
(cdr x
))
799 (setq sum
(+ sum
(math-factor-contains (car x
) a
))))
801 (if (and (eq (car-safe x
) '^
)
803 (* (math-factor-contains (nth 1 x
) a
) (nth 2 x
))
810 ;;; Merge all quotients and expand/simplify the numerator
811 (defun calcFunc-nrat (expr)
812 (if (math-any-floats expr
)
813 (setq expr
(calcFunc-pfrac expr
)))
814 (if (or (math-vectorp expr
)
815 (assq (car-safe expr
) calc-tweak-eqn-table
))
816 (cons (car expr
) (mapcar 'calcFunc-nrat
(cdr expr
)))
817 (let* ((calc-prefer-frac t
)
818 (res (math-to-ratpoly expr
))
819 (num (math-simplify (math-sort-terms (calcFunc-expand (car res
)))))
820 (den (math-simplify (math-sort-terms (calcFunc-expand (cdr res
)))))
821 (g (math-poly-gcd num den
)))
823 (let ((num2 (math-poly-div num g
))
824 (den2 (math-poly-div den g
)))
825 (and (eq (cdr num2
) 0) (eq (cdr den2
) 0)
826 (setq num
(car num2
) den
(car den2
)))))
827 (math-simplify (math-div num den
)))))
829 ;;; Returns expressions (num . denom).
830 (defun math-to-ratpoly (expr)
831 (let ((res (math-to-ratpoly-rec expr
)))
832 (cons (math-simplify (car res
)) (math-simplify (cdr res
)))))
834 (defun math-to-ratpoly-rec (expr)
835 (cond ((Math-primp expr
)
837 ((memq (car expr
) '(+ -
))
838 (let ((r1 (math-to-ratpoly-rec (nth 1 expr
)))
839 (r2 (math-to-ratpoly-rec (nth 2 expr
))))
840 (if (equal (cdr r1
) (cdr r2
))
841 (cons (list (car expr
) (car r1
) (car r2
)) (cdr r1
))
843 (cons (list (car expr
)
844 (math-mul (car r1
) (cdr r2
))
848 (cons (list (car expr
)
850 (math-mul (car r2
) (cdr r1
)))
852 (let ((g (math-poly-gcd (cdr r1
) (cdr r2
))))
853 (let ((d1 (and (not (eq g
1)) (math-poly-div (cdr r1
) g
)))
854 (d2 (and (not (eq g
1)) (math-poly-div
855 (math-mul (car r1
) (cdr r2
))
857 (if (and (eq (cdr d1
) 0) (eq (cdr d2
) 0))
858 (cons (list (car expr
) (car d2
)
859 (math-mul (car r2
) (car d1
)))
860 (math-mul (car d1
) (cdr r2
)))
861 (cons (list (car expr
)
862 (math-mul (car r1
) (cdr r2
))
863 (math-mul (car r2
) (cdr r1
)))
864 (math-mul (cdr r1
) (cdr r2
)))))))))))
866 (let* ((r1 (math-to-ratpoly-rec (nth 1 expr
)))
867 (r2 (math-to-ratpoly-rec (nth 2 expr
)))
868 (g (math-mul (math-poly-gcd (car r1
) (cdr r2
))
869 (math-poly-gcd (cdr r1
) (car r2
)))))
871 (cons (math-mul (car r1
) (car r2
))
872 (math-mul (cdr r1
) (cdr r2
)))
873 (cons (math-poly-div-exact (math-mul (car r1
) (car r2
)) g
)
874 (math-poly-div-exact (math-mul (cdr r1
) (cdr r2
)) g
)))))
876 (let* ((r1 (math-to-ratpoly-rec (nth 1 expr
)))
877 (r2 (math-to-ratpoly-rec (nth 2 expr
))))
878 (if (and (eq (cdr r1
) 1) (eq (cdr r2
) 1))
879 (cons (car r1
) (car r2
))
880 (let ((g (math-mul (math-poly-gcd (car r1
) (car r2
))
881 (math-poly-gcd (cdr r1
) (cdr r2
)))))
883 (cons (math-mul (car r1
) (cdr r2
))
884 (math-mul (cdr r1
) (car r2
)))
885 (cons (math-poly-div-exact (math-mul (car r1
) (cdr r2
)) g
)
886 (math-poly-div-exact (math-mul (cdr r1
) (car r2
))
888 ((and (eq (car expr
) '^
) (integerp (nth 2 expr
)))
889 (let ((r1 (math-to-ratpoly-rec (nth 1 expr
))))
890 (if (> (nth 2 expr
) 0)
891 (cons (math-pow (car r1
) (nth 2 expr
))
892 (math-pow (cdr r1
) (nth 2 expr
)))
893 (cons (math-pow (cdr r1
) (- (nth 2 expr
)))
894 (math-pow (car r1
) (- (nth 2 expr
)))))))
895 ((eq (car expr
) 'neg
)
896 (let ((r1 (math-to-ratpoly-rec (nth 1 expr
))))
897 (cons (math-neg (car r1
)) (cdr r1
))))
901 (defun math-ratpoly-p (expr &optional var
)
902 (cond ((equal expr var
) 1)
903 ((Math-primp expr
) 0)
904 ((memq (car expr
) '(+ -
))
905 (let ((p1 (math-ratpoly-p (nth 1 expr
) var
))
907 (and p1
(setq p2
(math-ratpoly-p (nth 2 expr
) var
))
910 (let ((p1 (math-ratpoly-p (nth 1 expr
) var
))
912 (and p1
(setq p2
(math-ratpoly-p (nth 2 expr
) var
))
914 ((eq (car expr
) 'neg
)
915 (math-ratpoly-p (nth 1 expr
) var
))
917 (let ((p1 (math-ratpoly-p (nth 1 expr
) var
))
919 (and p1
(setq p2
(math-ratpoly-p (nth 2 expr
) var
))
921 ((and (eq (car expr
) '^
)
922 (integerp (nth 2 expr
)))
923 (let ((p1 (math-ratpoly-p (nth 1 expr
) var
)))
924 (and p1
(* p1
(nth 2 expr
)))))
926 ((math-poly-depends expr var
) nil
)
930 (defun calcFunc-apart (expr &optional var
)
931 (cond ((Math-primp expr
) expr
)
933 (math-add (calcFunc-apart (nth 1 expr
) var
)
934 (calcFunc-apart (nth 2 expr
) var
)))
936 (math-sub (calcFunc-apart (nth 1 expr
) var
)
937 (calcFunc-apart (nth 2 expr
) var
)))
938 ((not (math-ratpoly-p expr var
))
939 (math-reject-arg expr
"Expected a rational function"))
941 (let* ((calc-prefer-frac t
)
942 (rat (math-to-ratpoly expr
))
945 (qr (math-poly-div num den
))
949 (setq var
(math-polynomial-base den
)))
950 (math-add q
(or (and var
951 (math-expr-contains den var
)
952 (math-partial-fractions r den var
))
953 (math-div r den
)))))))
956 (defun math-padded-polynomial (expr var deg
)
957 (let ((p (math-is-polynomial expr var deg
)))
958 (append p
(make-list (- deg
(length p
)) 0))))
960 (defun math-partial-fractions (r den var
)
961 (let* ((fden (calcFunc-factors den var
))
962 (tdeg (math-polynomial-p den var
))
967 (tz (make-list (1- tdeg
) 0))
968 (calc-matrix-mode 'scalar
))
969 (and (not (and (= (length fden
) 2) (eq (nth 2 (nth 1 fden
)) 1)))
971 (while (setq fp
(cdr fp
))
972 (let ((rpt (nth 2 (car fp
)))
973 (deg (math-polynomial-p (nth 1 (car fp
)) var
))
979 (setq dvar
(append '(vec) lz
'(1) tz
)
983 dnum
(math-add dnum
(math-mul dvar
984 (math-pow var deg2
)))
985 dlist
(cons (and (= deg2
(1- deg
))
986 (math-pow (nth 1 (car fp
)) rpt
))
990 (while (setq fpp
(cdr fpp
))
992 (setq mult
(math-mul mult
993 (math-pow (nth 1 (car fpp
))
994 (nth 2 (car fpp
)))))))
995 (setq dnum
(math-mul dnum mult
)))
996 (setq eqns
(math-add eqns
(math-mul dnum
1002 (setq eqns
(math-div (cons 'vec
(math-padded-polynomial r var tdeg
))
1008 (cons 'vec
(math-padded-polynomial
1011 (and (math-vectorp eqns
)
1014 (setq eqns
(nreverse eqns
))
1016 (setq num
(cons (car eqns
) num
)
1019 (setq num
(math-build-polynomial-expr
1021 res
(math-add res
(math-div num
(car dlist
)))
1023 (setq dlist
(cdr dlist
)))
1024 (math-normalize res
)))))))
1028 (defun math-expand-term (expr)
1029 (cond ((and (eq (car-safe expr
) '*)
1030 (memq (car-safe (nth 1 expr
)) '(+ -
)))
1031 (math-add-or-sub (list '* (nth 1 (nth 1 expr
)) (nth 2 expr
))
1032 (list '* (nth 2 (nth 1 expr
)) (nth 2 expr
))
1033 nil
(eq (car (nth 1 expr
)) '-
)))
1034 ((and (eq (car-safe expr
) '*)
1035 (memq (car-safe (nth 2 expr
)) '(+ -
)))
1036 (math-add-or-sub (list '* (nth 1 expr
) (nth 1 (nth 2 expr
)))
1037 (list '* (nth 1 expr
) (nth 2 (nth 2 expr
)))
1038 nil
(eq (car (nth 2 expr
)) '-
)))
1039 ((and (eq (car-safe expr
) '/)
1040 (memq (car-safe (nth 1 expr
)) '(+ -
)))
1041 (math-add-or-sub (list '/ (nth 1 (nth 1 expr
)) (nth 2 expr
))
1042 (list '/ (nth 2 (nth 1 expr
)) (nth 2 expr
))
1043 nil
(eq (car (nth 1 expr
)) '-
)))
1044 ((and (eq (car-safe expr
) '^
)
1045 (memq (car-safe (nth 1 expr
)) '(+ -
))
1046 (integerp (nth 2 expr
))
1047 (if (> (nth 2 expr
) 0)
1048 (or (and (or (> math-mt-many
500000) (< math-mt-many -
500000))
1049 (math-expand-power (nth 1 expr
) (nth 2 expr
)
1053 (list '^
(nth 1 expr
) (1- (nth 2 expr
)))))
1054 (if (< (nth 2 expr
) 0)
1055 (list '/ 1 (list '^
(nth 1 expr
) (- (nth 2 expr
))))))))
1058 (defun calcFunc-expand (expr &optional many
)
1059 (math-normalize (math-map-tree 'math-expand-term expr many
)))
1061 (defun math-expand-power (x n
&optional var else-nil
)
1062 (or (and (natnump n
)
1063 (memq (car-safe x
) '(+ -
))
1066 (while (memq (car-safe x
) '(+ -
))
1067 (setq terms
(cons (if (eq (car x
) '-
)
1068 (math-neg (nth 2 x
))
1072 (setq terms
(cons x terms
))
1076 (or (math-expr-contains (car p
) var
)
1077 (setq terms
(delq (car p
) terms
)
1078 cterms
(cons (car p
) cterms
)))
1081 (setq terms
(cons (apply 'calcFunc-add cterms
)
1083 (if (= (length terms
) 2)
1087 (setq accum
(list '+ accum
1088 (list '* (calcFunc-choose n i
)
1090 (list '^
(nth 1 terms
) i
)
1091 (list '^
(car terms
)
1100 (setq accum
(list '+ accum
1101 (list '^
(car p1
) 2))
1103 (while (setq p2
(cdr p2
))
1104 (setq accum
(list '+ accum
1115 (setq accum
(list '+ accum
(list '^
(car p1
) 3))
1117 (while (setq p2
(cdr p2
))
1118 (setq accum
(list '+
1124 (list '^
(car p1
) 2)
1129 (list '^
(car p2
) 2))))
1131 (while (setq p3
(cdr p3
))
1132 (setq accum
(list '+ accum
1144 (defun calcFunc-expandpow (x n
)
1145 (math-normalize (math-expand-power x n
)))
1147 ;;; arch-tag: d2566c51-2ccc-45f1-8c50-f3462c2953ff
1148 ;;; calc-poly.el ends here