1;;; calc-poly.el --- polynomial functions for Calc
2
3;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
4;;   2005, 2006, 2007 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 2, 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)
814  (let ((n 1)
815	(dp nil))
816    (while (setq p (cdr p))
817      (setq dp (cons (math-mul (car p) n) dp)
818	    n (1+ n)))
819    (nreverse dp)))
820
821(defun math-factor-contains (x a)
822  (if (equal x a)
823      1
824    (if (memq (car-safe x) '(+ - * / neg))
825	(let ((sum 0))
826	  (while (setq x (cdr x))
827	    (setq sum (+ sum (math-factor-contains (car x) a))))
828	  sum)
829      (if (and (eq (car-safe x) '^)
830	       (natnump (nth 2 x)))
831	  (* (math-factor-contains (nth 1 x) a) (nth 2 x))
832	0))))
833
834
835
836
837
838;;; Merge all quotients and expand/simplify the numerator
839(defun calcFunc-nrat (expr)
840  (if (math-any-floats expr)
841      (setq expr (calcFunc-pfrac expr)))
842  (if (or (math-vectorp expr)
843	  (assq (car-safe expr) calc-tweak-eqn-table))
844      (cons (car expr) (mapcar 'calcFunc-nrat (cdr expr)))
845    (let* ((calc-prefer-frac t)
846	   (res (math-to-ratpoly expr))
847	   (num (math-simplify (math-sort-terms (calcFunc-expand (car res)))))
848	   (den (math-simplify (math-sort-terms (calcFunc-expand (cdr res)))))
849	   (g (math-poly-gcd num den)))
850      (or (eq g 1)
851	  (let ((num2 (math-poly-div num g))
852		(den2 (math-poly-div den g)))
853	    (and (eq (cdr num2) 0) (eq (cdr den2) 0)
854		 (setq num (car num2) den (car den2)))))
855      (math-simplify (math-div num den)))))
856
857;;; Returns expressions (num . denom).
858(defun math-to-ratpoly (expr)
859  (let ((res (math-to-ratpoly-rec expr)))
860    (cons (math-simplify (car res)) (math-simplify (cdr res)))))
861
862(defun math-to-ratpoly-rec (expr)
863  (cond ((Math-primp expr)
864	 (cons expr 1))
865	((memq (car expr) '(+ -))
866	 (let ((r1 (math-to-ratpoly-rec (nth 1 expr)))
867	       (r2 (math-to-ratpoly-rec (nth 2 expr))))
868	   (if (equal (cdr r1) (cdr r2))
869	       (cons (list (car expr) (car r1) (car r2)) (cdr r1))
870	     (if (eq (cdr r1) 1)
871		 (cons (list (car expr)
872			     (math-mul (car r1) (cdr r2))
873			     (car r2))
874		       (cdr r2))
875	       (if (eq (cdr r2) 1)
876		   (cons (list (car expr)
877			       (car r1)
878			       (math-mul (car r2) (cdr r1)))
879			 (cdr r1))
880		 (let ((g (math-poly-gcd (cdr r1) (cdr r2))))
881		   (let ((d1 (and (not (eq g 1)) (math-poly-div (cdr r1) g)))
882			 (d2 (and (not (eq g 1)) (math-poly-div
883						  (math-mul (car r1) (cdr r2))
884						  g))))
885		     (if (and (eq (cdr d1) 0) (eq (cdr d2) 0))
886			 (cons (list (car expr) (car d2)
887				     (math-mul (car r2) (car d1)))
888			       (math-mul (car d1) (cdr r2)))
889		       (cons (list (car expr)
890				   (math-mul (car r1) (cdr r2))
891				   (math-mul (car r2) (cdr r1)))
892			     (math-mul (cdr r1) (cdr r2)))))))))))
893	((eq (car expr) '*)
894	 (let* ((r1 (math-to-ratpoly-rec (nth 1 expr)))
895		(r2 (math-to-ratpoly-rec (nth 2 expr)))
896		(g (math-mul (math-poly-gcd (car r1) (cdr r2))
897			     (math-poly-gcd (cdr r1) (car r2)))))
898	   (if (eq g 1)
899	       (cons (math-mul (car r1) (car r2))
900		     (math-mul (cdr r1) (cdr r2)))
901	     (cons (math-poly-div-exact (math-mul (car r1) (car r2)) g)
902		   (math-poly-div-exact (math-mul (cdr r1) (cdr r2)) g)))))
903	((eq (car expr) '/)
904	 (let* ((r1 (math-to-ratpoly-rec (nth 1 expr)))
905		(r2 (math-to-ratpoly-rec (nth 2 expr))))
906	   (if (and (eq (cdr r1) 1) (eq (cdr r2) 1))
907	       (cons (car r1) (car r2))
908	     (let ((g (math-mul (math-poly-gcd (car r1) (car r2))
909				(math-poly-gcd (cdr r1) (cdr r2)))))
910	       (if (eq g 1)
911		   (cons (math-mul (car r1) (cdr r2))
912			 (math-mul (cdr r1) (car r2)))
913		 (cons (math-poly-div-exact (math-mul (car r1) (cdr r2)) g)
914		       (math-poly-div-exact (math-mul (cdr r1) (car r2))
915					    g)))))))
916	((and (eq (car expr) '^) (integerp (nth 2 expr)))
917	 (let ((r1 (math-to-ratpoly-rec (nth 1 expr))))
918	   (if (> (nth 2 expr) 0)
919	       (cons (math-pow (car r1) (nth 2 expr))
920		     (math-pow (cdr r1) (nth 2 expr)))
921	     (cons (math-pow (cdr r1) (- (nth 2 expr)))
922		   (math-pow (car r1) (- (nth 2 expr)))))))
923	((eq (car expr) 'neg)
924	 (let ((r1 (math-to-ratpoly-rec (nth 1 expr))))
925	   (cons (math-neg (car r1)) (cdr r1))))
926	(t (cons expr 1))))
927
928
929(defun math-ratpoly-p (expr &optional var)
930  (cond ((equal expr var) 1)
931	((Math-primp expr) 0)
932	((memq (car expr) '(+ -))
933	 (let ((p1 (math-ratpoly-p (nth 1 expr) var))
934	       p2)
935	   (and p1 (setq p2 (math-ratpoly-p (nth 2 expr) var))
936		(max p1 p2))))
937	((eq (car expr) '*)
938	 (let ((p1 (math-ratpoly-p (nth 1 expr) var))
939	       p2)
940	   (and p1 (setq p2 (math-ratpoly-p (nth 2 expr) var))
941		(+ p1 p2))))
942	((eq (car expr) 'neg)
943	 (math-ratpoly-p (nth 1 expr) var))
944	((eq (car expr) '/)
945	 (let ((p1 (math-ratpoly-p (nth 1 expr) var))
946	       p2)
947	   (and p1 (setq p2 (math-ratpoly-p (nth 2 expr) var))
948		(- p1 p2))))
949	((and (eq (car expr) '^)
950	      (integerp (nth 2 expr)))
951	 (let ((p1 (math-ratpoly-p (nth 1 expr) var)))
952	   (and p1 (* p1 (nth 2 expr)))))
953	((not var) 1)
954	((math-poly-depends expr var) nil)
955	(t 0)))
956
957
958(defun calcFunc-apart (expr &optional var)
959  (cond ((Math-primp expr) expr)
960	((eq (car expr) '+)
961	 (math-add (calcFunc-apart (nth 1 expr) var)
962		   (calcFunc-apart (nth 2 expr) var)))
963	((eq (car expr) '-)
964	 (math-sub (calcFunc-apart (nth 1 expr) var)
965		   (calcFunc-apart (nth 2 expr) var)))
966	((not (math-ratpoly-p expr var))
967	 (math-reject-arg expr "Expected a rational function"))
968	(t
969	 (let* ((calc-prefer-frac t)
970		(rat (math-to-ratpoly expr))
971		(num (car rat))
972		(den (cdr rat))
973		(qr (math-poly-div num den))
974		(q (car qr))
975		(r (cdr qr)))
976	   (or var
977	       (setq var (math-polynomial-base 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)))))))
982
983
984(defun math-padded-polynomial (expr var deg)
985  (let ((p (math-is-polynomial expr var deg)))
986    (append p (make-list (- deg (length p)) 0))))
987
988(defun math-partial-fractions (r den var)
989  (let* ((fden (calcFunc-factors den var))
990	 (tdeg (math-polynomial-p den var))
991	 (fp fden)
992	 (dlist nil)
993	 (eqns 0)
994	 (lz nil)
995	 (tz (make-list (1- tdeg) 0))
996	 (calc-matrix-mode 'scalar))
997    (and (not (and (= (length fden) 2) (eq (nth 2 (nth 1 fden)) 1)))
998	 (progn
999	   (while (setq fp (cdr fp))
1000	     (let ((rpt (nth 2 (car fp)))
1001		   (deg (math-polynomial-p (nth 1 (car fp)) var))
1002		   dnum dvar deg2)
1003	       (while (> rpt 0)
1004		 (setq deg2 deg
1005		       dnum 0)
1006		 (while (> deg2 0)
1007		   (setq dvar (append '(vec) lz '(1) tz)
1008			 lz (cons 0 lz)
1009			 tz (cdr tz)
1010			 deg2 (1- deg2)
1011			 dnum (math-add dnum (math-mul dvar
1012						       (math-pow var deg2)))
1013			 dlist (cons (and (= deg2 (1- deg))
1014					  (math-pow (nth 1 (car fp)) rpt))
1015				     dlist)))
1016		 (let ((fpp fden)
1017		       (mult 1))
1018		   (while (setq fpp (cdr fpp))
1019		     (or (eq fpp fp)
1020			 (setq mult (math-mul mult
1021					      (math-pow (nth 1 (car fpp))
1022							(nth 2 (car fpp)))))))
1023		   (setq dnum (math-mul dnum mult)))
1024		 (setq eqns (math-add eqns (math-mul dnum
1025						     (math-pow
1026						      (nth 1 (car fp))
1027						      (- (nth 2 (car fp))
1028							 rpt))))
1029		       rpt (1- rpt)))))
1030	   (setq eqns (math-div (cons 'vec (math-padded-polynomial r var tdeg))
1031				(math-transpose
1032				 (cons 'vec
1033				       (mapcar
1034					(function
1035					 (lambda (x)
1036					   (cons 'vec (math-padded-polynomial
1037						       x var tdeg))))
1038					(cdr eqns))))))
1039	   (and (math-vectorp eqns)
1040		(let ((res 0)
1041		      (num nil))
1042		  (setq eqns (nreverse eqns))
1043		  (while eqns
1044		    (setq num (cons (car eqns) num)
1045			  eqns (cdr eqns))
1046		    (if (car dlist)
1047			(setq num (math-build-polynomial-expr
1048				   (nreverse num) var)
1049			      res (math-add res (math-div num (car dlist)))
1050			      num nil))
1051		    (setq dlist (cdr dlist)))
1052		  (math-normalize res)))))))
1053
1054
1055
1056(defun math-expand-term (expr)
1057  (cond ((and (eq (car-safe expr) '*)
1058	      (memq (car-safe (nth 1 expr)) '(+ -)))
1059	 (math-add-or-sub (list '* (nth 1 (nth 1 expr)) (nth 2 expr))
1060			  (list '* (nth 2 (nth 1 expr)) (nth 2 expr))
1061			  nil (eq (car (nth 1 expr)) '-)))
1062	((and (eq (car-safe expr) '*)
1063	      (memq (car-safe (nth 2 expr)) '(+ -)))
1064	 (math-add-or-sub (list '* (nth 1 expr) (nth 1 (nth 2 expr)))
1065			  (list '* (nth 1 expr) (nth 2 (nth 2 expr)))
1066			  nil (eq (car (nth 2 expr)) '-)))
1067	((and (eq (car-safe expr) '/)
1068	      (memq (car-safe (nth 1 expr)) '(+ -)))
1069	 (math-add-or-sub (list '/ (nth 1 (nth 1 expr)) (nth 2 expr))
1070			  (list '/ (nth 2 (nth 1 expr)) (nth 2 expr))
1071			  nil (eq (car (nth 1 expr)) '-)))
1072	((and (eq (car-safe expr) '^)
1073	      (memq (car-safe (nth 1 expr)) '(+ -))
1074	      (integerp (nth 2 expr))
1075              (if (and
1076                   (or (math-known-matrixp (nth 1 (nth 1 expr)))
1077                       (math-known-matrixp (nth 2 (nth 1 expr)))
1078                       (and
1079                        calc-matrix-mode
1080                        (not (eq calc-matrix-mode 'scalar))
1081                        (not (and (math-known-scalarp (nth 1 (nth 1 expr)))
1082                                  (math-known-scalarp (nth 2 (nth 1 expr)))))))
1083                   (> (nth 2 expr) 1))
1084                  (if (= (nth 2 expr) 2)
1085                      (math-add-or-sub (list '* (nth 1 (nth 1 expr)) (nth 1 expr))
1086                                       (list '* (nth 2 (nth 1 expr)) (nth 1 expr))
1087                                       nil (eq (car (nth 1 expr)) '-))
1088                    (math-add-or-sub (list '* (nth 1 (nth 1 expr))
1089                                           (list '^ (nth 1 expr)
1090                                                 (1- (nth 2 expr))))
1091                                     (list '* (nth 2 (nth 1 expr))
1092                                           (list '^ (nth 1 expr)
1093                                                 (1- (nth 2 expr))))
1094                                     nil (eq (car (nth 1 expr)) '-)))
1095                (if (> (nth 2 expr) 0)
1096                    (or (and (or (> math-mt-many 500000) (< math-mt-many -500000))
1097                             (math-expand-power (nth 1 expr) (nth 2 expr)
1098                                                nil t))
1099                        (list '*
1100                              (nth 1 expr)
1101                              (list '^ (nth 1 expr) (1- (nth 2 expr)))))
1102                  (if (< (nth 2 expr) 0)
1103                      (list '/ 1 (list '^ (nth 1 expr) (- (nth 2 expr)))))))))
1104	(t expr)))
1105
1106(defun calcFunc-expand (expr &optional many)
1107  (math-normalize (math-map-tree 'math-expand-term expr many)))
1108
1109(defun math-expand-power (x n &optional var else-nil)
1110  (or (and (natnump n)
1111	   (memq (car-safe x) '(+ -))
1112	   (let ((terms nil)
1113		 (cterms nil))
1114	     (while (memq (car-safe x) '(+ -))
1115	       (setq terms (cons (if (eq (car x) '-)
1116				     (math-neg (nth 2 x))
1117				   (nth 2 x))
1118				 terms)
1119		     x (nth 1 x)))
1120	     (setq terms (cons x terms))
1121	     (if var
1122		 (let ((p terms))
1123		   (while p
1124		     (or (math-expr-contains (car p) var)
1125			 (setq terms (delq (car p) terms)
1126			       cterms (cons (car p) cterms)))
1127		     (setq p (cdr p)))
1128		   (if cterms
1129		       (setq terms (cons (apply 'calcFunc-add cterms)
1130					 terms)))))
1131	     (if (= (length terms) 2)
1132		 (let ((i 0)
1133		       (accum 0))
1134		   (while (<= i n)
1135		     (setq accum (list '+ accum
1136				       (list '* (calcFunc-choose n i)
1137					     (list '*
1138						   (list '^ (nth 1 terms) i)
1139						   (list '^ (car terms)
1140							 (- n i)))))
1141			   i (1+ i)))
1142		   accum)
1143	       (if (= n 2)
1144		   (let ((accum 0)
1145			 (p1 terms)
1146			 p2)
1147		     (while p1
1148		       (setq accum (list '+ accum
1149					 (list '^ (car p1) 2))
1150			     p2 p1)
1151		       (while (setq p2 (cdr p2))
1152			 (setq accum (list '+ accum
1153					   (list '* 2 (list '*
1154							    (car p1)
1155							    (car p2))))))
1156		       (setq p1 (cdr p1)))
1157		     accum)
1158		 (if (= n 3)
1159		     (let ((accum 0)
1160			   (p1 terms)
1161			   p2 p3)
1162		       (while p1
1163			 (setq accum (list '+ accum (list '^ (car p1) 3))
1164			       p2 p1)
1165			 (while (setq p2 (cdr p2))
1166			   (setq accum (list '+
1167					     (list '+
1168						   accum
1169						   (list '* 3
1170							 (list
1171							  '*
1172							  (list '^ (car p1) 2)
1173							  (car p2))))
1174					     (list '* 3
1175						   (list
1176						    '* (car p1)
1177						    (list '^ (car p2) 2))))
1178				 p3 p2)
1179			   (while (setq p3 (cdr p3))
1180			     (setq accum (list '+ accum
1181					       (list '* 6
1182						     (list '*
1183							   (car p1)
1184							   (list
1185							    '* (car p2)
1186							    (car p3))))))))
1187			 (setq p1 (cdr p1)))
1188		       accum))))))
1189      (and (not else-nil)
1190	   (list '^ x n))))
1191
1192(defun calcFunc-expandpow (x n)
1193  (math-normalize (math-expand-power x n)))
1194
1195(provide 'calc-poly)
1196
1197;;; arch-tag: d2566c51-2ccc-45f1-8c50-f3462c2953ff
1198;;; calc-poly.el ends here
1199