calc-poly.el

Revision 4220, 38.5 kB (checked in by miyoshi, 9 months ago)
Sync up with Emacs22.2.

Line
1	;;; calc-poly.el --- polynomial functions for Calc
2
3	;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
4	;; 2005, 2006, 2007, 2008 Free Software Foundation, Inc.
5
6	;; Author: David Gillespie <daveg@synaptics.com>
7	;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
8
9	;; This file is part of GNU Emacs.
10
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, or (at your option)
14	;; any later version.
15
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.
20
21	;; You should have received a copy of the GNU General Public License
22	;; along with GNU Emacs; see the file COPYING. If not, write to the
23	;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
24	;; Boston, MA 02110-1301, USA.
25
26	;;; Commentary:
27
28	;;; Code:
29
30	;; This file is autoloaded from calc-ext.el.
31
32	(require 'calc-ext)
33	(require 'calc-macs)
34
35	(defun calcFunc-pcont (expr &optional var)
36	(cond ((Math-primp expr)
37	(cond ((Math-zerop expr) 1)
38	((Math-messy-integerp expr) (math-trunc expr))
39	((Math-objectp expr) expr)
40	((or (equal expr var) (not var)) 1)
41	(t expr)))
42	((eq (car expr) '*)
43	(math-mul (calcFunc-pcont (nth 1 expr) var)
44	(calcFunc-pcont (nth 2 expr) var)))
45	((eq (car expr) '/)
46	(math-div (calcFunc-pcont (nth 1 expr) var)
47	(calcFunc-pcont (nth 2 expr) var)))
48	((and (eq (car expr) '^) (Math-natnump (nth 2 expr)))
49	(math-pow (calcFunc-pcont (nth 1 expr) var) (nth 2 expr)))
50	((memq (car expr) '(neg polar))
51	(calcFunc-pcont (nth 1 expr) var))
52	((consp var)
53	(let ((p (math-is-polynomial expr var)))
54	(if p
55	(let ((lead (nth (1- (length p)) p))
56	(cont (math-poly-gcd-list p)))
57	(if (math-guess-if-neg lead)
58	(math-neg cont)
59	cont))
60	1)))
61	((memq (car expr) '(+ - cplx sdev))
62	(let ((cont (calcFunc-pcont (nth 1 expr) var)))
63	(if (eq cont 1)
64	1
65	(let ((c2 (calcFunc-pcont (nth 2 expr) var)))
66	(if (and (math-negp cont)
67	(if (eq (car expr) '-) (math-posp c2) (math-negp c2)))
68	(math-neg (math-poly-gcd cont c2))
69	(math-poly-gcd cont c2))))))
70	(var expr)
71	(t 1)))
72
73	(defun calcFunc-pprim (expr &optional var)
74	(let ((cont (calcFunc-pcont expr var)))
75	(if (math-equal-int cont 1)
76	expr
77	(math-poly-div-exact expr cont var))))
78
79	(defun math-div-poly-const (expr c)
80	(cond ((memq (car-safe expr) '(+ -))
81	(list (car expr)
82	(math-div-poly-const (nth 1 expr) c)
83	(math-div-poly-const (nth 2 expr) c)))
84	(t (math-div expr c))))
85
86	(defun calcFunc-pdeg (expr &optional var)
87	(if (Math-zerop expr)
88	'(neg (var inf var-inf))
89	(if var
90	(or (math-polynomial-p expr var)
91	(math-reject-arg expr "Expected a polynomial"))
92	(math-poly-degree expr))))
93
94	(defun math-poly-degree (expr)
95	(cond ((Math-primp expr)
96	(if (eq (car-safe expr) 'var) 1 0))
97	((eq (car expr) 'neg)
98	(math-poly-degree (nth 1 expr)))
99	((eq (car expr) '*)
100	(+ (math-poly-degree (nth 1 expr))
101	(math-poly-degree (nth 2 expr))))
102	((eq (car expr) '/)
103	(- (math-poly-degree (nth 1 expr))
104	(math-poly-degree (nth 2 expr))))
105	((and (eq (car expr) '^) (natnump (nth 2 expr)))
106	(* (math-poly-degree (nth 1 expr)) (nth 2 expr)))
107	((memq (car expr) '(+ -))
108	(max (math-poly-degree (nth 1 expr))
109	(math-poly-degree (nth 2 expr))))
110	(t 1)))
111
112	(defun calcFunc-plead (expr var)
113	(cond ((eq (car-safe expr) '*)
114	(math-mul (calcFunc-plead (nth 1 expr) var)
115	(calcFunc-plead (nth 2 expr) var)))
116	((eq (car-safe expr) '/)
117	(math-div (calcFunc-plead (nth 1 expr) var)
118	(calcFunc-plead (nth 2 expr) var)))
119	((and (eq (car-safe expr) '^) (math-natnump (nth 2 expr)))
120	(math-pow (calcFunc-plead (nth 1 expr) var) (nth 2 expr)))
121	((Math-primp expr)
122	(if (equal expr var)
123	1
124	expr))
125	(t
126	(let ((p (math-is-polynomial expr var)))
127	(if (cdr p)
128	(nth (1- (length p)) p)
129	1)))))
130
131
132
133
134
135	;;; Polynomial quotient, remainder, and GCD.
136	;;; Originally by Ove Ewerlid (ewerlid@mizar.DoCS.UU.SE).
137	;;; Modifications and simplifications by daveg.
138
139	(defvar math-poly-modulus 1)
140
141	;;; Return gcd of two polynomials
142	(defun calcFunc-pgcd (pn pd)
143	(if (math-any-floats pn)
144	(math-reject-arg pn "Coefficients must be rational"))
145	(if (math-any-floats pd)
146	(math-reject-arg pd "Coefficients must be rational"))
147	(let ((calc-prefer-frac t)
148	(math-poly-modulus (math-poly-modulus pn pd)))
149	(math-poly-gcd pn pd)))
150
151	;;; Return only quotient to top of stack (nil if zero)
152
153	;; calc-poly-div-remainder is a local variable for
154	;; calc-poly-div (in calc-alg.el), but is used by
155	;; calcFunc-pdiv, which is called by calc-poly-div.
156	(defvar calc-poly-div-remainder)
157
158	(defun calcFunc-pdiv (pn pd &optional base)
159	(let* ((calc-prefer-frac t)
160	(math-poly-modulus (math-poly-modulus pn pd))
161	(res (math-poly-div pn pd base)))
162	(setq calc-poly-div-remainder (cdr res))
163	(car res)))
164
165	;;; Return only remainder to top of stack
166	(defun calcFunc-prem (pn pd &optional base)
167	(let ((calc-prefer-frac t)
168	(math-poly-modulus (math-poly-modulus pn pd)))
169	(cdr (math-poly-div pn pd base))))
170
171	(defun calcFunc-pdivrem (pn pd &optional base)
172	(let* ((calc-prefer-frac t)
173	(math-poly-modulus (math-poly-modulus pn pd))
174	(res (math-poly-div pn pd base)))
175	(list 'vec (car res) (cdr res))))
176
177	(defun calcFunc-pdivide (pn pd &optional base)
178	(let* ((calc-prefer-frac t)
179	(math-poly-modulus (math-poly-modulus pn pd))
180	(res (math-poly-div pn pd base)))
181	(math-add (car res) (math-div (cdr res) pd))))
182
183
184	;;; Multiply two terms, expanding out products of sums.
185	(defun math-mul-thru (lhs rhs)
186	(if (memq (car-safe lhs) '(+ -))
187	(list (car lhs)
188	(math-mul-thru (nth 1 lhs) rhs)
189	(math-mul-thru (nth 2 lhs) rhs))
190	(if (memq (car-safe rhs) '(+ -))
191	(list (car rhs)
192	(math-mul-thru lhs (nth 1 rhs))
193	(math-mul-thru lhs (nth 2 rhs)))
194	(math-mul lhs rhs))))
195
196	(defun math-div-thru (num den)
197	(if (memq (car-safe num) '(+ -))
198	(list (car num)
199	(math-div-thru (nth 1 num) den)
200	(math-div-thru (nth 2 num) den))
201	(math-div num den)))
202
203
204	;;; Sort the terms of a sum into canonical order.
205	(defun math-sort-terms (expr)
206	(if (memq (car-safe expr) '(+ -))
207	(math-list-to-sum
208	(sort (math-sum-to-list expr)
209	(function (lambda (a b) (math-beforep (car a) (car b))))))
210	expr))
211
212	(defun math-list-to-sum (lst)
213	(if (cdr lst)
214	(list (if (cdr (car lst)) '- '+)
215	(math-list-to-sum (cdr lst))
216	(car (car lst)))
217	(if (cdr (car lst))
218	(math-neg (car (car lst)))
219	(car (car lst)))))
220
221	(defun math-sum-to-list (tree &optional neg)
222	(cond ((eq (car-safe tree) '+)
223	(nconc (math-sum-to-list (nth 1 tree) neg)
224	(math-sum-to-list (nth 2 tree) neg)))
225	((eq (car-safe tree) '-)
226	(nconc (math-sum-to-list (nth 1 tree) neg)
227	(math-sum-to-list (nth 2 tree) (not neg))))
228	(t (list (cons tree neg)))))
229
230	;;; Check if the polynomial coefficients are modulo forms.
231	(defun math-poly-modulus (expr &optional expr2)
232	(or (math-poly-modulus-rec expr)
233	(and expr2 (math-poly-modulus-rec expr2))
234	1))
235
236	(defun math-poly-modulus-rec (expr)
237	(if (and (eq (car-safe expr) 'mod) (Math-natnump (nth 2 expr)))
238	(list 'mod 1 (nth 2 expr))
239	(and (memq (car-safe expr) '(+ - * /))
240	(or (math-poly-modulus-rec (nth 1 expr))
241	(math-poly-modulus-rec (nth 2 expr))))))
242
243
244	;;; Divide two polynomials. Return (quotient . remainder).
245	(defvar math-poly-div-base nil)
246	(defun math-poly-div (u v &optional math-poly-div-base)
247	(if math-poly-div-base
248	(math-do-poly-div u v)
249	(math-do-poly-div (calcFunc-expand u) (calcFunc-expand v))))
250
251	(defun math-poly-div-exact (u v &optional base)
252	(let ((res (math-poly-div u v base)))
253	(if (eq (cdr res) 0)
254	(car res)
255	(math-reject-arg (list 'vec u v) "Argument is not a polynomial"))))
256
257	(defun math-do-poly-div (u v)
258	(cond ((math-constp u)
259	(if (math-constp v)
260	(cons (math-div u v) 0)
261	(cons 0 u)))
262	((math-constp v)
263	(cons (if (eq v 1)
264	u
265	(if (memq (car-safe u) '(+ -))
266	(math-add-or-sub (math-poly-div-exact (nth 1 u) v)
267	(math-poly-div-exact (nth 2 u) v)
268	nil (eq (car u) '-))
269	(math-div u v)))
270	0))
271	((Math-equal u v)
272	(cons math-poly-modulus 0))
273	((and (math-atomic-factorp u) (math-atomic-factorp v))
274	(cons (math-simplify (math-div u v)) 0))
275	(t
276	(let ((base (or math-poly-div-base
277	(math-poly-div-base u v)))
278	vp up res)
279	(if (or (null base)
280	(null (setq vp (math-is-polynomial v base nil 'gen))))
281	(cons 0 u)
282	(setq up (math-is-polynomial u base nil 'gen)
283	res (math-poly-div-coefs up vp))
284	(cons (math-build-polynomial-expr (car res) base)
285	(math-build-polynomial-expr (cdr res) base)))))))
286
287	(defun math-poly-div-rec (u v)
288	(cond ((math-constp u)
289	(math-div u v))
290	((math-constp v)
291	(if (eq v 1)
292	u
293	(if (memq (car-safe u) '(+ -))
294	(math-add-or-sub (math-poly-div-rec (nth 1 u) v)
295	(math-poly-div-rec (nth 2 u) v)
296	nil (eq (car u) '-))
297	(math-div u v))))
298	((Math-equal u v) math-poly-modulus)
299	((and (math-atomic-factorp u) (math-atomic-factorp v))
300	(math-simplify (math-div u v)))
301	(math-poly-div-base
302	(math-div u v))
303	(t
304	(let ((base (math-poly-div-base u v))
305	vp up res)
306	(if (or (null base)
307	(null (setq vp (math-is-polynomial v base nil 'gen))))
308	(math-div u v)
309	(setq up (math-is-polynomial u base nil 'gen)
310	res (math-poly-div-coefs up vp))
311	(math-add (math-build-polynomial-expr (car res) base)
312	(math-div (math-build-polynomial-expr (cdr res) base)
313	v)))))))
314
315	;;; Divide two polynomials in coefficient-list form. Return (quot . rem).
316	(defun math-poly-div-coefs (u v)
317	(cond ((null v) (math-reject-arg nil "Division by zero"))
318	((< (length u) (length v)) (cons nil u))
319	((cdr u)
320	(let ((q nil)
321	(urev (reverse u))
322	(vrev (reverse v)))
323	(while
324	(let ((qk (math-poly-div-rec (math-simplify (car urev))
325	(car vrev)))
326	(up urev)
327	(vp vrev))
328	(if (or q (not (math-zerop qk)))
329	(setq q (cons qk q)))
330	(while (setq up (cdr up) vp (cdr vp))
331	(setcar up (math-sub (car up) (math-mul-thru qk (car vp)))))
332	(setq urev (cdr urev))
333	up))
334	(while (and urev (Math-zerop (car urev)))
335	(setq urev (cdr urev)))
336	(cons q (nreverse (mapcar 'math-simplify urev)))))
337	(t
338	(cons (list (math-poly-div-rec (car u) (car v)))
339	nil))))
340
341	;;; Perform a pseudo-division of polynomials. (See Knuth section 4.6.1.)
342	;;; This returns only the remainder from the pseudo-division.
343	(defun math-poly-pseudo-div (u v)
344	(cond ((null v) nil)
345	((< (length u) (length v)) u)
346	((or (cdr u) (cdr v))
347	(let ((urev (reverse u))
348	(vrev (reverse v))
349	up)
350	(while
351	(let ((vp vrev))
352	(setq up urev)
353	(while (setq up (cdr up) vp (cdr vp))
354	(setcar up (math-sub (math-mul-thru (car vrev) (car up))
355	(math-mul-thru (car urev) (car vp)))))
356	(setq urev (cdr urev))
357	up)
358	(while up
359	(setcar up (math-mul-thru (car vrev) (car up)))
360	(setq up (cdr up))))
361	(while (and urev (Math-zerop (car urev)))
362	(setq urev (cdr urev)))
363	(nreverse (mapcar 'math-simplify urev))))
364	(t nil)))
365
366	;;; Compute the GCD of two multivariate polynomials.
367	(defun math-poly-gcd (u v)
368	(cond ((Math-equal u v) u)
369	((math-constp u)
370	(if (Math-zerop u)
371	v
372	(calcFunc-gcd u (calcFunc-pcont v))))
373	((math-constp v)
374	(if (Math-zerop v)
375	v
376	(calcFunc-gcd v (calcFunc-pcont u))))
377	(t
378	(let ((base (math-poly-gcd-base u v)))
379	(if base
380	(math-simplify
381	(calcFunc-expand
382	(math-build-polynomial-expr
383	(math-poly-gcd-coefs (math-is-polynomial u base nil 'gen)
384	(math-is-polynomial v base nil 'gen))
385	base)))
386	(calcFunc-gcd (calcFunc-pcont u) (calcFunc-pcont u)))))))
387
388	(defun math-poly-div-list (lst a)
389	(if (eq a 1)
390	lst
391	(if (eq a -1)
392	(math-mul-list lst a)
393	(mapcar (function (lambda (x) (math-poly-div-exact x a))) lst))))
394
395	(defun math-mul-list (lst a)
396	(if (eq a 1)
397	lst
398	(if (eq a -1)
399	(mapcar 'math-neg lst)
400	(and (not (eq a 0))
401	(mapcar (function (lambda (x) (math-mul x a))) lst)))))
402
403	;;; Run GCD on all elements in a list.
404	(defun math-poly-gcd-list (lst)
405	(if (or (memq 1 lst) (memq -1 lst))
406	(math-poly-gcd-frac-list lst)
407	(let ((gcd (car lst)))
408	(while (and (setq lst (cdr lst)) (not (eq gcd 1)))
409	(or (eq (car lst) 0)
410	(setq gcd (math-poly-gcd gcd (car lst)))))
411	(if lst (setq lst (math-poly-gcd-frac-list lst)))
412	gcd)))
413
414	(defun math-poly-gcd-frac-list (lst)
415	(while (and lst (not (eq (car-safe (car lst)) 'frac)))
416	(setq lst (cdr lst)))
417	(if lst
418	(let ((denom (nth 2 (car lst))))
419	(while (setq lst (cdr lst))
420	(if (eq (car-safe (car lst)) 'frac)
421	(setq denom (calcFunc-lcm denom (nth 2 (car lst))))))
422	(list 'frac 1 denom))
423	1))
424
425	;;; Compute the GCD of two monovariate polynomial lists.
426	;;; Knuth section 4.6.1, algorithm C.
427	(defun math-poly-gcd-coefs (u v)
428	(let ((d (math-poly-gcd (math-poly-gcd-list u)
429	(math-poly-gcd-list v)))
430	(g 1) (h 1) (z 0) hh r delta ghd)
431	(while (and u v (Math-zerop (car u)) (Math-zerop (car v)))
432	(setq u (cdr u) v (cdr v) z (1+ z)))
433	(or (eq d 1)
434	(setq u (math-poly-div-list u d)
435	v (math-poly-div-list v d)))
436	(while (progn
437	(setq delta (- (length u) (length v)))
438	(if (< delta 0)
439	(setq r u u v v r delta (- delta)))
440	(setq r (math-poly-pseudo-div u v))
441	(cdr r))
442	(setq u v
443	v (math-poly-div-list r (math-mul g (math-pow h delta)))
444	g (nth (1- (length u)) u)
445	h (if (<= delta 1)
446	(math-mul (math-pow g delta) (math-pow h (- 1 delta)))
447	(math-poly-div-exact (math-pow g delta)
448	(math-pow h (1- delta))))))
449	(setq v (if r
450	(list d)
451	(math-mul-list (math-poly-div-list v (math-poly-gcd-list v)) d)))
452	(if (math-guess-if-neg (nth (1- (length v)) v))
453	(setq v (math-mul-list v -1)))
454	(while (>= (setq z (1- z)) 0)
455	(setq v (cons 0 v)))
456	v))
457
458
459	;;; Return true if is a factor containing no sums or quotients.
460	(defun math-atomic-factorp (expr)
461	(cond ((eq (car-safe expr) '*)
462	(and (math-atomic-factorp (nth 1 expr))
463	(math-atomic-factorp (nth 2 expr))))
464	((memq (car-safe expr) '(+ - /))
465	nil)
466	((memq (car-safe expr) '(^ neg))
467	(math-atomic-factorp (nth 1 expr)))
468	(t t)))
469
470	;;; Find a suitable base for dividing a by b.
471	;;; The base must exist in both expressions.
472	;;; The degree in the numerator must be higher or equal than the
473	;;; degree in the denominator.
474	;;; If the above conditions are not met the quotient is just a remainder.
475	;;; Return nil if this is the case.
476
477	(defun math-poly-div-base (a b)
478	(let (a-base b-base)
479	(and (setq a-base (math-total-polynomial-base a))
480	(setq b-base (math-total-polynomial-base b))
481	(catch 'return
482	(while a-base
483	(let ((maybe (assoc (car (car a-base)) b-base)))
484	(if maybe
485	(if (>= (nth 1 (car a-base)) (nth 1 maybe))
486	(throw 'return (car (car a-base))))))
487	(setq a-base (cdr a-base)))))))
488
489	;;; Same as above but for gcd algorithm.
490	;;; Here there is no requirement that degree(a) > degree(b).
491	;;; Take the base that has the highest degree considering both a and b.
492	;;; ("a^20+b^21+x^3+a+b", "a+b^2+x^5+a^22+b^10") --> (a 22)
493
494	(defun math-poly-gcd-base (a b)
495	(let (a-base b-base)
496	(and (setq a-base (math-total-polynomial-base a))
497	(setq b-base (math-total-polynomial-base b))
498	(catch 'return
499	(while (and a-base b-base)
500	(if (> (nth 1 (car a-base)) (nth 1 (car b-base)))
501	(if (assoc (car (car a-base)) b-base)
502	(throw 'return (car (car a-base)))
503	(setq a-base (cdr a-base)))
504	(if (assoc (car (car b-base)) a-base)
505	(throw 'return (car (car b-base)))
506	(setq b-base (cdr b-base)))))))))
507
508	;;; Sort a list of polynomial bases.
509	(defun math-sort-poly-base-list (lst)
510	(sort lst (function (lambda (a b)
511	(or (> (nth 1 a) (nth 1 b))
512	(and (= (nth 1 a) (nth 1 b))
513	(math-beforep (car a) (car b))))))))
514
515	;;; Given an expression find all variables that are polynomial bases.
516	;;; Return list in the form '( (var1 degree1) (var2 degree2) ... ).
517
518	;; The variable math-poly-base-total-base is local to
519	;; math-total-polynomial-base, but is used by math-polynomial-p1,
520	;; which is called by math-total-polynomial-base.
521	(defvar math-poly-base-total-base)
522
523	(defun math-total-polynomial-base (expr)
524	(let ((math-poly-base-total-base nil))
525	(math-polynomial-base expr 'math-polynomial-p1)
526	(math-sort-poly-base-list math-poly-base-total-base)))
527
528	;; The variable math-poly-base-top-expr is local to math-polynomial-base
529	;; in calc-alg.el, but is used by math-polynomial-p1 which is called
530	;; by math-polynomial-base.
531	(defvar math-poly-base-top-expr)
532
533	(defun math-polynomial-p1 (subexpr)
534	(or (assoc subexpr math-poly-base-total-base)
535	(memq (car subexpr) '(+ - * / neg))
536	(and (eq (car subexpr) '^) (natnump (nth 2 subexpr)))
537	(let* ((math-poly-base-variable subexpr)
538	(exponent (math-polynomial-p math-poly-base-top-expr subexpr)))
539	(if exponent
540	(setq math-poly-base-total-base (cons (list subexpr exponent)
541	math-poly-base-total-base)))))
542	nil)
543
544	;; The variable math-factored-vars is local to calcFunc-factors and
545	;; calcFunc-factor, but is used by math-factor-expr and
546	;; math-factor-expr-part, which are called (directly and indirectly) by
547	;; calcFunc-factor and calcFunc-factors.
548	(defvar math-factored-vars)
549
550	;; The variable math-fact-expr is local to calcFunc-factors,
551	;; calcFunc-factor and math-factor-expr, but is used by math-factor-expr-try
552	;; and math-factor-expr-part, which are called (directly and indirectly) by
553	;; calcFunc-factor, calcFunc-factors and math-factor-expr.
554	(defvar math-fact-expr)
555
556	;; The variable math-to-list is local to calcFunc-factors and
557	;; calcFunc-factor, but is used by math-accum-factors, which is
558	;; called (indirectly) by calcFunc-factors and calcFunc-factor.
559	(defvar math-to-list)
560
561	(defun calcFunc-factors (math-fact-expr &optional var)
562	(let ((math-factored-vars (if var t nil))
563	(math-to-list t)
564	(calc-prefer-frac t))
565	(or var
566	(setq var (math-polynomial-base math-fact-expr)))
567	(let ((res (math-factor-finish
568	(or (catch 'factor (math-factor-expr-try var))
569	math-fact-expr))))
570	(math-simplify (if (math-vectorp res)
571	res
572	(list 'vec (list 'vec res 1)))))))
573
574	(defun calcFunc-factor (math-fact-expr &optional var)
575	(let ((math-factored-vars nil)
576	(math-to-list nil)
577	(calc-prefer-frac t))
578	(math-simplify (math-factor-finish
579	(if var
580	(let ((math-factored-vars t))
581	(or (catch 'factor (math-factor-expr-try var)) math-fact-expr))
582	(math-factor-expr math-fact-expr))))))
583
584	(defun math-factor-finish (x)
585	(if (Math-primp x)
586	x
587	(if (eq (car x) 'calcFunc-Fac-Prot)
588	(math-factor-finish (nth 1 x))
589	(cons (car x) (mapcar 'math-factor-finish (cdr x))))))
590
591	(defun math-factor-protect (x)
592	(if (memq (car-safe x) '(+ -))
593	(list 'calcFunc-Fac-Prot x)
594	x))
595
596	(defun math-factor-expr (math-fact-expr)
597	(cond ((eq math-factored-vars t) math-fact-expr)
598	((or (memq (car-safe math-fact-expr) '(* / ^ neg))
599	(assq (car-safe math-fact-expr) calc-tweak-eqn-table))
600	(cons (car math-fact-expr) (mapcar 'math-factor-expr (cdr math-fact-expr))))
601	((memq (car-safe math-fact-expr) '(+ -))
602	(let* ((math-factored-vars math-factored-vars)
603	(y (catch 'factor (math-factor-expr-part math-fact-expr))))
604	(if y
605	(math-factor-expr y)
606	math-fact-expr)))
607	(t math-fact-expr)))
608
609	(defun math-factor-expr-part (x) ; uses "expr"
610	(if (memq (car-safe x) '(+ - * / ^ neg))
611	(while (setq x (cdr x))
612	(math-factor-expr-part (car x)))
613	(and (not (Math-objvecp x))
614	(not (assoc x math-factored-vars))
615	(> (math-factor-contains math-fact-expr x) 1)
616	(setq math-factored-vars (cons (list x) math-factored-vars))
617	(math-factor-expr-try x))))
618
619	;; The variable math-fet-x is local to math-factor-expr-try, but is
620	;; used by math-factor-poly-coefs, which is called by math-factor-expr-try.
621	(defvar math-fet-x)
622
623	(defun math-factor-expr-try (math-fet-x)
624	(if (eq (car-safe math-fact-expr) '*)
625	(let ((res1 (catch 'factor (let ((math-fact-expr (nth 1 math-fact-expr)))
626	(math-factor-expr-try math-fet-x))))
627	(res2 (catch 'factor (let ((math-fact-expr (nth 2 math-fact-expr)))
628	(math-factor-expr-try math-fet-x)))))
629	(and (or res1 res2)
630	(throw 'factor (math-accum-factors (or res1 (nth 1 math-fact-expr)) 1
631	(or res2 (nth 2 math-fact-expr))))))
632	(let* ((p (math-is-polynomial math-fact-expr math-fet-x 30 'gen))
633	(math-poly-modulus (math-poly-modulus math-fact-expr))
634	res)
635	(and (cdr p)
636	(setq res (math-factor-poly-coefs p))
637	(throw 'factor res)))))
638
639	(defun math-accum-factors (fac pow facs)
640	(if math-to-list
641	(if (math-vectorp fac)
642	(progn
643	(while (setq fac (cdr fac))
644	(setq facs (math-accum-factors (nth 1 (car fac))
645	(* pow (nth 2 (car fac)))
646	facs)))
647	facs)
648	(if (and (eq (car-safe fac) '^) (natnump (nth 2 fac)))
649	(setq pow (* pow (nth 2 fac))
650	fac (nth 1 fac)))
651	(if (eq fac 1)
652	facs
653	(or (math-vectorp facs)
654	(setq facs (if (eq facs 1) '(vec)
655	(list 'vec (list 'vec facs 1)))))
656	(let ((found facs))
657	(while (and (setq found (cdr found))
658	(not (equal fac (nth 1 (car found))))))
659	(if found
660	(progn
661	(setcar (cdr (cdr (car found))) (+ pow (nth 2 (car found))))
662	facs)
663	;; Put constant term first.
664	(if (and (cdr facs) (Math-ratp (nth 1 (nth 1 facs))))
665	(cons 'vec (cons (nth 1 facs) (cons (list 'vec fac pow)
666	(cdr (cdr facs)))))
667	(cons 'vec (cons (list 'vec fac pow) (cdr facs))))))))
668	(math-mul (math-pow fac pow) facs)))
669
670	(defun math-factor-poly-coefs (p &optional square-free) ; uses "x"
671	(let (t1 t2 temp)
672	(cond ((not (cdr p))
673	(or (car p) 0))
674
675	;; Strip off multiples of math-fet-x.
676	((Math-zerop (car p))
677	(let ((z 0))
678	(while (and p (Math-zerop (car p)))
679	(setq z (1+ z) p (cdr p)))
680	(if (cdr p)
681	(setq p (math-factor-poly-coefs p square-free))
682	(setq p (math-sort-terms (math-factor-expr (car p)))))
683	(math-accum-factors math-fet-x z (math-factor-protect p))))
684
685	;; Factor out content.
686	((and (not square-free)
687	(not (eq 1 (setq t1 (math-mul (math-poly-gcd-list p)
688	(if (math-guess-if-neg
689	(nth (1- (length p)) p))
690	-1 1))))))
691	(math-accum-factors t1 1 (math-factor-poly-coefs
692	(math-poly-div-list p t1) 'cont)))
693
694	;; Check if linear in math-fet-x.
695	((not (cdr (cdr p)))
696	(math-sort-terms
697	(math-add (math-factor-protect
698	(math-sort-terms
699	(math-factor-expr (car p))))
700	(math-mul math-fet-x (math-factor-protect
701	(math-sort-terms
702	(math-factor-expr (nth 1 p))))))))
703
704	;; If symbolic coefficients, use FactorRules.
705	((let ((pp p))
706	(while (and pp (or (Math-ratp (car pp))
707	(and (eq (car (car pp)) 'mod)
708	(Math-integerp (nth 1 (car pp)))
709	(Math-integerp (nth 2 (car pp))))))
710	(setq pp (cdr pp)))
711	pp)
712	(let ((res (math-rewrite
713	(list 'calcFunc-thecoefs math-fet-x (cons 'vec p))
714	'(var FactorRules var-FactorRules))))
715	(or (and (eq (car-safe res) 'calcFunc-thefactors)
716	(= (length res) 3)
717	(math-vectorp (nth 2 res))
718	(let ((facs 1)
719	(vec (nth 2 res)))
720	(while (setq vec (cdr vec))
721	(setq facs (math-accum-factors (car vec) 1 facs)))
722	facs))
723	(math-build-polynomial-expr p math-fet-x))))
724
725	;; Check if rational coefficients (i.e., not modulo a prime).
726	((eq math-poly-modulus 1)
727
728	;; Check if there are any squared terms, or a content not = 1.
729	(if (or (eq square-free t)
730	(equal (setq t1 (math-poly-gcd-coefs
731	p (setq t2 (math-poly-deriv-coefs p))))
732	'(1)))
733
734	;; We now have a square-free polynomial with integer coefs.
735	;; For now, we use a kludgey method that finds linear and
736	;; quadratic terms using floating-point root-finding.
737	(if (setq t1 (let ((calc-symbolic-mode nil))
738	(math-poly-all-roots nil p t)))
739	(let ((roots (car t1))
740	(csign (if (math-negp (nth (1- (length p)) p)) -1 1))
741	(expr 1)
742	(unfac (nth 1 t1))
743	(scale (nth 2 t1)))
744	(while roots
745	(let ((coef0 (car (car roots)))
746	(coef1 (cdr (car roots))))
747	(setq expr (math-accum-factors
748	(if coef1
749	(let ((den (math-lcm-denoms
750	coef0 coef1)))
751	(setq scale (math-div scale den))
752	(math-add
753	(math-add
754	(math-mul den (math-pow math-fet-x 2))
755	(math-mul (math-mul coef1 den)
756	math-fet-x))
757	(math-mul coef0 den)))
758	(let ((den (math-lcm-denoms coef0)))
759	(setq scale (math-div scale den))
760	(math-add (math-mul den math-fet-x)
761	(math-mul coef0 den))))
762	1 expr)
763	roots (cdr roots))))
764	(setq expr (math-accum-factors
765	expr 1
766	(math-mul csign
767	(math-build-polynomial-expr
768	(math-mul-list (nth 1 t1) scale)
769	math-fet-x)))))
770	(math-build-polynomial-expr p math-fet-x)) ; can't factor it.
771
772	;; Separate out the squared terms (Knuth exercise 4.6.2-34).
773	;; This step also divides out the content of the polynomial.
774	(let* ((cabs (math-poly-gcd-list p))
775	(csign (if (math-negp (nth (1- (length p)) p)) -1 1))
776	(t1s (math-mul-list t1 csign))
777	(uu nil)
778	(v (car (math-poly-div-coefs p t1s)))
779	(w (car (math-poly-div-coefs t2 t1s))))
780	(while
781	(not (math-poly-zerop
782	(setq t2 (math-poly-simplify
783	(math-poly-mix
784	w 1 (math-poly-deriv-coefs v) -1)))))
785	(setq t1 (math-poly-gcd-coefs v t2)
786	uu (cons t1 uu)
787	v (car (math-poly-div-coefs v t1))
788	w (car (math-poly-div-coefs t2 t1))))
789	(setq t1 (length uu)
790	t2 (math-accum-factors (math-factor-poly-coefs v t)
791	(1+ t1) 1))
792	(while uu
793	(setq t2 (math-accum-factors (math-factor-poly-coefs
794	(car uu) t)
795	t1 t2)
796	t1 (1- t1)
797	uu (cdr uu)))
798	(math-accum-factors (math-mul cabs csign) 1 t2))))
799
800	;; Factoring modulo a prime.
801	((and (= (length (setq temp (math-poly-gcd-coefs
802	p (math-poly-deriv-coefs p))))
803	(length p)))
804	(setq p (car temp))
805	(while (cdr temp)
806	(setq temp (nthcdr (nth 2 math-poly-modulus) temp)
807	p (cons (car temp) p)))
808	(and (setq temp (math-factor-poly-coefs p))
809	(math-pow temp (nth 2 math-poly-modulus))))
810	(t
811	(math-reject-arg nil *"Modulo factorization not yet implemented"**)))))
812
813	(defun math-poly-deriv-coefs (p)