Remove unnecessary stack overflow dependency
[emacs.git] / lisp / calc / calc-poly.el
blob1dab3c474aa0f097ce8175352924cad71f933f2f
1 ;;; calc-poly.el --- polynomial functions for Calc
3 ;; Copyright (C) 1990-1993, 2001-2015 Free Software Foundation, Inc.
5 ;; Author: David Gillespie <daveg@synaptics.com>
6 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
8 ;; This file is part of GNU Emacs.
10 ;; GNU Emacs is free software: you can redistribute it and/or modify
11 ;; it under the terms of the GNU General Public License as published by
12 ;; the Free Software Foundation, either version 3 of the License, or
13 ;; (at your option) any later version.
15 ;; GNU Emacs is distributed in the hope that it will be useful,
16 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
17 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 ;; GNU General Public License for more details.
20 ;; You should have received a copy of the GNU General Public License
21 ;; along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>.
23 ;;; Commentary:
25 ;;; Code:
27 ;; This file is autoloaded from calc-ext.el.
29 (require 'calc-ext)
30 (require 'calc-macs)
32 (defun calcFunc-pcont (expr &optional var)
33 (cond ((Math-primp expr)
34 (cond ((Math-zerop expr) 1)
35 ((Math-messy-integerp expr) (math-trunc expr))
36 ((Math-objectp expr) expr)
37 ((or (equal expr var) (not var)) 1)
38 (t expr)))
39 ((eq (car expr) '*)
40 (math-mul (calcFunc-pcont (nth 1 expr) var)
41 (calcFunc-pcont (nth 2 expr) var)))
42 ((eq (car expr) '/)
43 (math-div (calcFunc-pcont (nth 1 expr) var)
44 (calcFunc-pcont (nth 2 expr) var)))
45 ((and (eq (car expr) '^) (Math-natnump (nth 2 expr)))
46 (math-pow (calcFunc-pcont (nth 1 expr) var) (nth 2 expr)))
47 ((memq (car expr) '(neg polar))
48 (calcFunc-pcont (nth 1 expr) var))
49 ((consp var)
50 (let ((p (math-is-polynomial expr var)))
51 (if p
52 (let ((lead (nth (1- (length p)) p))
53 (cont (math-poly-gcd-list p)))
54 (if (math-guess-if-neg lead)
55 (math-neg cont)
56 cont))
57 1)))
58 ((memq (car expr) '(+ - cplx sdev))
59 (let ((cont (calcFunc-pcont (nth 1 expr) var)))
60 (if (eq cont 1)
62 (let ((c2 (calcFunc-pcont (nth 2 expr) var)))
63 (if (and (math-negp cont)
64 (if (eq (car expr) '-) (math-posp c2) (math-negp c2)))
65 (math-neg (math-poly-gcd cont c2))
66 (math-poly-gcd cont c2))))))
67 (var expr)
68 (t 1)))
70 (defun calcFunc-pprim (expr &optional var)
71 (let ((cont (calcFunc-pcont expr var)))
72 (if (math-equal-int cont 1)
73 expr
74 (math-poly-div-exact expr cont var))))
76 (defun math-div-poly-const (expr c)
77 (cond ((memq (car-safe expr) '(+ -))
78 (list (car expr)
79 (math-div-poly-const (nth 1 expr) c)
80 (math-div-poly-const (nth 2 expr) c)))
81 (t (math-div expr c))))
83 (defun calcFunc-pdeg (expr &optional var)
84 (if (Math-zerop expr)
85 '(neg (var inf var-inf))
86 (if var
87 (or (math-polynomial-p expr var)
88 (math-reject-arg expr "Expected a polynomial"))
89 (math-poly-degree expr))))
91 (defun math-poly-degree (expr)
92 (cond ((Math-primp expr)
93 (if (eq (car-safe expr) 'var) 1 0))
94 ((eq (car expr) 'neg)
95 (math-poly-degree (nth 1 expr)))
96 ((eq (car expr) '*)
97 (+ (math-poly-degree (nth 1 expr))
98 (math-poly-degree (nth 2 expr))))
99 ((eq (car expr) '/)
100 (- (math-poly-degree (nth 1 expr))
101 (math-poly-degree (nth 2 expr))))
102 ((and (eq (car expr) '^) (natnump (nth 2 expr)))
103 (* (math-poly-degree (nth 1 expr)) (nth 2 expr)))
104 ((memq (car expr) '(+ -))
105 (max (math-poly-degree (nth 1 expr))
106 (math-poly-degree (nth 2 expr))))
107 (t 1)))
109 (defun calcFunc-plead (expr var)
110 (cond ((eq (car-safe expr) '*)
111 (math-mul (calcFunc-plead (nth 1 expr) var)
112 (calcFunc-plead (nth 2 expr) var)))
113 ((eq (car-safe expr) '/)
114 (math-div (calcFunc-plead (nth 1 expr) var)
115 (calcFunc-plead (nth 2 expr) var)))
116 ((and (eq (car-safe expr) '^) (math-natnump (nth 2 expr)))
117 (math-pow (calcFunc-plead (nth 1 expr) var) (nth 2 expr)))
118 ((Math-primp expr)
119 (if (equal expr var)
121 expr))
123 (let ((p (math-is-polynomial expr var)))
124 (if (cdr p)
125 (nth (1- (length p)) p)
126 1)))))
132 ;;; Polynomial quotient, remainder, and GCD.
133 ;;; Originally by Ove Ewerlid (ewerlid@mizar.DoCS.UU.SE).
134 ;;; Modifications and simplifications by daveg.
136 (defvar math-poly-modulus 1)
138 ;;; Return gcd of two polynomials
139 (defun calcFunc-pgcd (pn pd)
140 (if (math-any-floats pn)
141 (math-reject-arg pn "Coefficients must be rational"))
142 (if (math-any-floats pd)
143 (math-reject-arg pd "Coefficients must be rational"))
144 (let ((calc-prefer-frac t)
145 (math-poly-modulus (math-poly-modulus pn pd)))
146 (math-poly-gcd pn pd)))
148 ;;; Return only quotient to top of stack (nil if zero)
150 ;; calc-poly-div-remainder is a local variable for
151 ;; calc-poly-div (in calc-alg.el), but is used by
152 ;; calcFunc-pdiv, which is called by calc-poly-div.
153 (defvar calc-poly-div-remainder)
155 (defun calcFunc-pdiv (pn pd &optional base)
156 (let* ((calc-prefer-frac t)
157 (math-poly-modulus (math-poly-modulus pn pd))
158 (res (math-poly-div pn pd base)))
159 (setq calc-poly-div-remainder (cdr res))
160 (car res)))
162 ;;; Return only remainder to top of stack
163 (defun calcFunc-prem (pn pd &optional base)
164 (let ((calc-prefer-frac t)
165 (math-poly-modulus (math-poly-modulus pn pd)))
166 (cdr (math-poly-div pn pd base))))
168 (defun calcFunc-pdivrem (pn pd &optional base)
169 (let* ((calc-prefer-frac t)
170 (math-poly-modulus (math-poly-modulus pn pd))
171 (res (math-poly-div pn pd base)))
172 (list 'vec (car res) (cdr res))))
174 (defun calcFunc-pdivide (pn pd &optional base)
175 (let* ((calc-prefer-frac t)
176 (math-poly-modulus (math-poly-modulus pn pd))
177 (res (math-poly-div pn pd base)))
178 (math-add (car res) (math-div (cdr res) pd))))
181 ;;; Multiply two terms, expanding out products of sums.
182 (defun math-mul-thru (lhs rhs)
183 (if (memq (car-safe lhs) '(+ -))
184 (list (car lhs)
185 (math-mul-thru (nth 1 lhs) rhs)
186 (math-mul-thru (nth 2 lhs) rhs))
187 (if (memq (car-safe rhs) '(+ -))
188 (list (car rhs)
189 (math-mul-thru lhs (nth 1 rhs))
190 (math-mul-thru lhs (nth 2 rhs)))
191 (math-mul lhs rhs))))
193 (defun math-div-thru (num den)
194 (if (memq (car-safe num) '(+ -))
195 (list (car num)
196 (math-div-thru (nth 1 num) den)
197 (math-div-thru (nth 2 num) den))
198 (math-div num den)))
201 ;;; Sort the terms of a sum into canonical order.
202 (defun math-sort-terms (expr)
203 (if (memq (car-safe expr) '(+ -))
204 (math-list-to-sum
205 (sort (math-sum-to-list expr)
206 (function (lambda (a b) (math-beforep (car a) (car b))))))
207 expr))
209 (defun math-list-to-sum (lst)
210 (if (cdr lst)
211 (list (if (cdr (car lst)) '- '+)
212 (math-list-to-sum (cdr lst))
213 (car (car lst)))
214 (if (cdr (car lst))
215 (math-neg (car (car lst)))
216 (car (car lst)))))
218 (defun math-sum-to-list (tree &optional neg)
219 (cond ((eq (car-safe tree) '+)
220 (nconc (math-sum-to-list (nth 1 tree) neg)
221 (math-sum-to-list (nth 2 tree) neg)))
222 ((eq (car-safe tree) '-)
223 (nconc (math-sum-to-list (nth 1 tree) neg)
224 (math-sum-to-list (nth 2 tree) (not neg))))
225 (t (list (cons tree neg)))))
227 ;;; Check if the polynomial coefficients are modulo forms.
228 (defun math-poly-modulus (expr &optional expr2)
229 (or (math-poly-modulus-rec expr)
230 (and expr2 (math-poly-modulus-rec expr2))
233 (defun math-poly-modulus-rec (expr)
234 (if (and (eq (car-safe expr) 'mod) (Math-natnump (nth 2 expr)))
235 (list 'mod 1 (nth 2 expr))
236 (and (memq (car-safe expr) '(+ - * /))
237 (or (math-poly-modulus-rec (nth 1 expr))
238 (math-poly-modulus-rec (nth 2 expr))))))
241 ;;; Divide two polynomials. Return (quotient . remainder).
242 (defvar math-poly-div-base nil)
243 (defun math-poly-div (u v &optional math-poly-div-base)
244 (if math-poly-div-base
245 (math-do-poly-div u v)
246 (math-do-poly-div (calcFunc-expand u) (calcFunc-expand v))))
248 (defun math-poly-div-exact (u v &optional base)
249 (let ((res (math-poly-div u v base)))
250 (if (eq (cdr res) 0)
251 (car res)
252 (math-reject-arg (list 'vec u v) "Argument is not a polynomial"))))
254 (defun math-do-poly-div (u v)
255 (cond ((math-constp u)
256 (if (math-constp v)
257 (cons (math-div u v) 0)
258 (cons 0 u)))
259 ((math-constp v)
260 (cons (if (eq v 1)
262 (if (memq (car-safe u) '(+ -))
263 (math-add-or-sub (math-poly-div-exact (nth 1 u) v)
264 (math-poly-div-exact (nth 2 u) v)
265 nil (eq (car u) '-))
266 (math-div u v)))
268 ((Math-equal u v)
269 (cons math-poly-modulus 0))
270 ((and (math-atomic-factorp u) (math-atomic-factorp v))
271 (cons (math-simplify (math-div u v)) 0))
273 (let ((base (or math-poly-div-base
274 (math-poly-div-base u v)))
275 vp up res)
276 (if (or (null base)
277 (null (setq vp (math-is-polynomial v base nil 'gen))))
278 (cons 0 u)
279 (setq up (math-is-polynomial u base nil 'gen)
280 res (math-poly-div-coefs up vp))
281 (cons (math-build-polynomial-expr (car res) base)
282 (math-build-polynomial-expr (cdr res) base)))))))
284 (defun math-poly-div-rec (u v)
285 (cond ((math-constp u)
286 (math-div u v))
287 ((math-constp v)
288 (if (eq v 1)
290 (if (memq (car-safe u) '(+ -))
291 (math-add-or-sub (math-poly-div-rec (nth 1 u) v)
292 (math-poly-div-rec (nth 2 u) v)
293 nil (eq (car u) '-))
294 (math-div u v))))
295 ((Math-equal u v) math-poly-modulus)
296 ((and (math-atomic-factorp u) (math-atomic-factorp v))
297 (math-simplify (math-div u v)))
298 (math-poly-div-base
299 (math-div u v))
301 (let ((base (math-poly-div-base u v))
302 vp up res)
303 (if (or (null base)
304 (null (setq vp (math-is-polynomial v base nil 'gen))))
305 (math-div u v)
306 (setq up (math-is-polynomial u base nil 'gen)
307 res (math-poly-div-coefs up vp))
308 (math-add (math-build-polynomial-expr (car res) base)
309 (math-div (math-build-polynomial-expr (cdr res) base)
310 v)))))))
312 ;;; Divide two polynomials in coefficient-list form. Return (quot . rem).
313 (defun math-poly-div-coefs (u v)
314 (cond ((null v) (math-reject-arg nil "Division by zero"))
315 ((< (length u) (length v)) (cons nil u))
316 ((cdr u)
317 (let ((q nil)
318 (urev (reverse u))
319 (vrev (reverse v)))
320 (while
321 (let ((qk (math-poly-div-rec (math-simplify (car urev))
322 (car vrev)))
323 (up urev)
324 (vp vrev))
325 (if (or q (not (math-zerop qk)))
326 (setq q (cons qk q)))
327 (while (setq up (cdr up) vp (cdr vp))
328 (setcar up (math-sub (car up) (math-mul-thru qk (car vp)))))
329 (setq urev (cdr urev))
330 up))
331 (while (and urev (Math-zerop (car urev)))
332 (setq urev (cdr urev)))
333 (cons q (nreverse (mapcar 'math-simplify urev)))))
335 (cons (list (math-poly-div-rec (car u) (car v)))
336 nil))))
338 ;;; Perform a pseudo-division of polynomials. (See Knuth section 4.6.1.)
339 ;;; This returns only the remainder from the pseudo-division.
340 (defun math-poly-pseudo-div (u v)
341 (cond ((null v) nil)
342 ((< (length u) (length v)) u)
343 ((or (cdr u) (cdr v))
344 (let ((urev (reverse u))
345 (vrev (reverse v))
347 (while
348 (let ((vp vrev))
349 (setq up urev)
350 (while (setq up (cdr up) vp (cdr vp))
351 (setcar up (math-sub (math-mul-thru (car vrev) (car up))
352 (math-mul-thru (car urev) (car vp)))))
353 (setq urev (cdr urev))
355 (while up
356 (setcar up (math-mul-thru (car vrev) (car up)))
357 (setq up (cdr up))))
358 (while (and urev (Math-zerop (car urev)))
359 (setq urev (cdr urev)))
360 (nreverse (mapcar 'math-simplify urev))))
361 (t nil)))
363 ;;; Compute the GCD of two multivariate polynomials.
364 (defun math-poly-gcd (u v)
365 (cond ((Math-equal u v) u)
366 ((math-constp u)
367 (if (Math-zerop u)
369 (calcFunc-gcd u (calcFunc-pcont v))))
370 ((math-constp v)
371 (if (Math-zerop v)
373 (calcFunc-gcd v (calcFunc-pcont u))))
375 (let ((base (math-poly-gcd-base u v)))
376 (if base
377 (math-simplify
378 (calcFunc-expand
379 (math-build-polynomial-expr
380 (math-poly-gcd-coefs (math-is-polynomial u base nil 'gen)
381 (math-is-polynomial v base nil 'gen))
382 base)))
383 (calcFunc-gcd (calcFunc-pcont u) (calcFunc-pcont u)))))))
385 (defun math-poly-div-list (lst a)
386 (if (eq a 1)
388 (if (eq a -1)
389 (math-mul-list lst a)
390 (mapcar (function (lambda (x) (math-poly-div-exact x a))) lst))))
392 (defun math-mul-list (lst a)
393 (if (eq a 1)
395 (if (eq a -1)
396 (mapcar 'math-neg lst)
397 (and (not (eq a 0))
398 (mapcar (function (lambda (x) (math-mul x a))) lst)))))
400 ;;; Run GCD on all elements in a list.
401 (defun math-poly-gcd-list (lst)
402 (if (or (memq 1 lst) (memq -1 lst))
403 (math-poly-gcd-frac-list lst)
404 (let ((gcd (car lst)))
405 (while (and (setq lst (cdr lst)) (not (eq gcd 1)))
406 (or (eq (car lst) 0)
407 (setq gcd (math-poly-gcd gcd (car lst)))))
408 (if lst (setq lst (math-poly-gcd-frac-list lst)))
409 gcd)))
411 (defun math-poly-gcd-frac-list (lst)
412 (while (and lst (not (eq (car-safe (car lst)) 'frac)))
413 (setq lst (cdr lst)))
414 (if lst
415 (let ((denom (nth 2 (car lst))))
416 (while (setq lst (cdr lst))
417 (if (eq (car-safe (car lst)) 'frac)
418 (setq denom (calcFunc-lcm denom (nth 2 (car lst))))))
419 (list 'frac 1 denom))
422 ;;; Compute the GCD of two univariate polynomial lists.
423 ;;; Knuth section 4.6.1, algorithm C.
424 (defun math-poly-gcd-coefs (u v)
425 (let ((d (math-poly-gcd (math-poly-gcd-list u)
426 (math-poly-gcd-list v)))
427 (g 1) (h 1) (z 0) hh r delta ghd)
428 (while (and u v (Math-zerop (car u)) (Math-zerop (car v)))
429 (setq u (cdr u) v (cdr v) z (1+ z)))
430 (or (eq d 1)
431 (setq u (math-poly-div-list u d)
432 v (math-poly-div-list v d)))
433 (while (progn
434 (setq delta (- (length u) (length v)))
435 (if (< delta 0)
436 (setq r u u v v r delta (- delta)))
437 (setq r (math-poly-pseudo-div u v))
438 (cdr r))
439 (setq u v
440 v (math-poly-div-list r (math-mul g (math-pow h delta)))
441 g (nth (1- (length u)) u)
442 h (if (<= delta 1)
443 (math-mul (math-pow g delta) (math-pow h (- 1 delta)))
444 (math-poly-div-exact (math-pow g delta)
445 (math-pow h (1- delta))))))
446 (setq v (if r
447 (list d)
448 (math-mul-list (math-poly-div-list v (math-poly-gcd-list v)) d)))
449 (if (math-guess-if-neg (nth (1- (length v)) v))
450 (setq v (math-mul-list v -1)))
451 (while (>= (setq z (1- z)) 0)
452 (setq v (cons 0 v)))
456 ;;; Return true if is a factor containing no sums or quotients.
457 (defun math-atomic-factorp (expr)
458 (cond ((eq (car-safe expr) '*)
459 (and (math-atomic-factorp (nth 1 expr))
460 (math-atomic-factorp (nth 2 expr))))
461 ((memq (car-safe expr) '(+ - /))
462 nil)
463 ((memq (car-safe expr) '(^ neg))
464 (math-atomic-factorp (nth 1 expr)))
465 (t t)))
467 ;;; Find a suitable base for dividing a by b.
468 ;;; The base must exist in both expressions.
469 ;;; The degree in the numerator must be higher or equal than the
470 ;;; degree in the denominator.
471 ;;; If the above conditions are not met the quotient is just a remainder.
472 ;;; Return nil if this is the case.
474 (defun math-poly-div-base (a b)
475 (let (a-base b-base)
476 (and (setq a-base (math-total-polynomial-base a))
477 (setq b-base (math-total-polynomial-base b))
478 (catch 'return
479 (while a-base
480 (let ((maybe (assoc (car (car a-base)) b-base)))
481 (if maybe
482 (if (>= (nth 1 (car a-base)) (nth 1 maybe))
483 (throw 'return (car (car a-base))))))
484 (setq a-base (cdr a-base)))))))
486 ;;; Same as above but for gcd algorithm.
487 ;;; Here there is no requirement that degree(a) > degree(b).
488 ;;; Take the base that has the highest degree considering both a and b.
489 ;;; ("a^20+b^21+x^3+a+b", "a+b^2+x^5+a^22+b^10") --> (a 22)
491 (defun math-poly-gcd-base (a b)
492 (let (a-base b-base)
493 (and (setq a-base (math-total-polynomial-base a))
494 (setq b-base (math-total-polynomial-base b))
495 (catch 'return
496 (while (and a-base b-base)
497 (if (> (nth 1 (car a-base)) (nth 1 (car b-base)))
498 (if (assoc (car (car a-base)) b-base)
499 (throw 'return (car (car a-base)))
500 (setq a-base (cdr a-base)))
501 (if (assoc (car (car b-base)) a-base)
502 (throw 'return (car (car b-base)))
503 (setq b-base (cdr b-base)))))))))
505 ;;; Sort a list of polynomial bases.
506 (defun math-sort-poly-base-list (lst)
507 (sort lst (function (lambda (a b)
508 (or (> (nth 1 a) (nth 1 b))
509 (and (= (nth 1 a) (nth 1 b))
510 (math-beforep (car a) (car b))))))))
512 ;;; Given an expression find all variables that are polynomial bases.
513 ;;; Return list in the form '( (var1 degree1) (var2 degree2) ... ).
515 ;; The variable math-poly-base-total-base is local to
516 ;; math-total-polynomial-base, but is used by math-polynomial-p1,
517 ;; which is called by math-total-polynomial-base.
518 (defvar math-poly-base-total-base)
520 (defun math-total-polynomial-base (expr)
521 (let ((math-poly-base-total-base nil))
522 (math-polynomial-base expr 'math-polynomial-p1)
523 (math-sort-poly-base-list math-poly-base-total-base)))
525 ;; The variable math-poly-base-top-expr is local to math-polynomial-base
526 ;; in calc-alg.el, but is used by math-polynomial-p1 which is called
527 ;; by math-polynomial-base.
528 (defvar math-poly-base-top-expr)
530 (defun math-polynomial-p1 (subexpr)
531 (or (assoc subexpr math-poly-base-total-base)
532 (memq (car subexpr) '(+ - * / neg))
533 (and (eq (car subexpr) '^) (natnump (nth 2 subexpr)))
534 (let* ((math-poly-base-variable subexpr)
535 (exponent (math-polynomial-p math-poly-base-top-expr subexpr)))
536 (if exponent
537 (setq math-poly-base-total-base (cons (list subexpr exponent)
538 math-poly-base-total-base)))))
539 nil)
541 ;; The variable math-factored-vars is local to calcFunc-factors and
542 ;; calcFunc-factor, but is used by math-factor-expr and
543 ;; math-factor-expr-part, which are called (directly and indirectly) by
544 ;; calcFunc-factor and calcFunc-factors.
545 (defvar math-factored-vars)
547 ;; The variable math-fact-expr is local to calcFunc-factors,
548 ;; calcFunc-factor and math-factor-expr, but is used by math-factor-expr-try
549 ;; and math-factor-expr-part, which are called (directly and indirectly) by
550 ;; calcFunc-factor, calcFunc-factors and math-factor-expr.
551 (defvar math-fact-expr)
553 ;; The variable math-to-list is local to calcFunc-factors and
554 ;; calcFunc-factor, but is used by math-accum-factors, which is
555 ;; called (indirectly) by calcFunc-factors and calcFunc-factor.
556 (defvar math-to-list)
558 (defun calcFunc-factors (math-fact-expr &optional var)
559 (let ((math-factored-vars (if var t nil))
560 (math-to-list t)
561 (calc-prefer-frac t))
562 (or var
563 (setq var (math-polynomial-base math-fact-expr)))
564 (let ((res (math-factor-finish
565 (or (catch 'factor (math-factor-expr-try var))
566 math-fact-expr))))
567 (math-simplify (if (math-vectorp res)
569 (list 'vec (list 'vec res 1)))))))
571 (defun calcFunc-factor (math-fact-expr &optional var)
572 (let ((math-factored-vars nil)
573 (math-to-list nil)
574 (calc-prefer-frac t))
575 (math-simplify (math-factor-finish
576 (if var
577 (let ((math-factored-vars t))
578 (or (catch 'factor (math-factor-expr-try var)) math-fact-expr))
579 (math-factor-expr math-fact-expr))))))
581 (defun math-factor-finish (x)
582 (if (Math-primp x)
584 (if (eq (car x) 'calcFunc-Fac-Prot)
585 (math-factor-finish (nth 1 x))
586 (cons (car x) (mapcar 'math-factor-finish (cdr x))))))
588 (defun math-factor-protect (x)
589 (if (memq (car-safe x) '(+ -))
590 (list 'calcFunc-Fac-Prot x)
593 (defun math-factor-expr (math-fact-expr)
594 (cond ((eq math-factored-vars t) math-fact-expr)
595 ((or (memq (car-safe math-fact-expr) '(* / ^ neg))
596 (assq (car-safe math-fact-expr) calc-tweak-eqn-table))
597 (cons (car math-fact-expr) (mapcar 'math-factor-expr (cdr math-fact-expr))))
598 ((memq (car-safe math-fact-expr) '(+ -))
599 (let* ((math-factored-vars math-factored-vars)
600 (y (catch 'factor (math-factor-expr-part math-fact-expr))))
601 (if y
602 (math-factor-expr y)
603 math-fact-expr)))
604 (t math-fact-expr)))
606 (defun math-factor-expr-part (x) ; uses "expr"
607 (if (memq (car-safe x) '(+ - * / ^ neg))
608 (while (setq x (cdr x))
609 (math-factor-expr-part (car x)))
610 (and (not (Math-objvecp x))
611 (not (assoc x math-factored-vars))
612 (> (math-factor-contains math-fact-expr x) 1)
613 (setq math-factored-vars (cons (list x) math-factored-vars))
614 (math-factor-expr-try x))))
616 ;; The variable math-fet-x is local to math-factor-expr-try, but is
617 ;; used by math-factor-poly-coefs, which is called by math-factor-expr-try.
618 (defvar math-fet-x)
620 (defun math-factor-expr-try (math-fet-x)
621 (if (eq (car-safe math-fact-expr) '*)
622 (let ((res1 (catch 'factor (let ((math-fact-expr (nth 1 math-fact-expr)))
623 (math-factor-expr-try math-fet-x))))
624 (res2 (catch 'factor (let ((math-fact-expr (nth 2 math-fact-expr)))
625 (math-factor-expr-try math-fet-x)))))
626 (and (or res1 res2)
627 (throw 'factor (math-accum-factors (or res1 (nth 1 math-fact-expr)) 1
628 (or res2 (nth 2 math-fact-expr))))))
629 (let* ((p (math-is-polynomial math-fact-expr math-fet-x 30 'gen))
630 (math-poly-modulus (math-poly-modulus math-fact-expr))
631 res)
632 (and (cdr p)
633 (setq res (math-factor-poly-coefs p))
634 (throw 'factor res)))))
636 (defun math-accum-factors (fac pow facs)
637 (if math-to-list
638 (if (math-vectorp fac)
639 (progn
640 (while (setq fac (cdr fac))
641 (setq facs (math-accum-factors (nth 1 (car fac))
642 (* pow (nth 2 (car fac)))
643 facs)))
644 facs)
645 (if (and (eq (car-safe fac) '^) (natnump (nth 2 fac)))
646 (setq pow (* pow (nth 2 fac))
647 fac (nth 1 fac)))
648 (if (eq fac 1)
649 facs
650 (or (math-vectorp facs)
651 (setq facs (if (eq facs 1) '(vec)
652 (list 'vec (list 'vec facs 1)))))
653 (let ((found facs))
654 (while (and (setq found (cdr found))
655 (not (equal fac (nth 1 (car found))))))
656 (if found
657 (progn
658 (setcar (cdr (cdr (car found))) (+ pow (nth 2 (car found))))
659 facs)
660 ;; Put constant term first.
661 (if (and (cdr facs) (Math-ratp (nth 1 (nth 1 facs))))
662 (cons 'vec (cons (nth 1 facs) (cons (list 'vec fac pow)
663 (cdr (cdr facs)))))
664 (cons 'vec (cons (list 'vec fac pow) (cdr facs))))))))
665 (math-mul (math-pow fac pow) (math-factor-protect facs))))
667 (defun math-factor-poly-coefs (p &optional square-free) ; uses "x"
668 (let (t1 t2 temp)
669 (cond ((not (cdr p))
670 (or (car p) 0))
672 ;; Strip off multiples of math-fet-x.
673 ((Math-zerop (car p))
674 (let ((z 0))
675 (while (and p (Math-zerop (car p)))
676 (setq z (1+ z) p (cdr p)))
677 (if (cdr p)
678 (setq p (math-factor-poly-coefs p square-free))
679 (setq p (math-sort-terms (math-factor-expr (car p)))))
680 (math-accum-factors math-fet-x z (math-factor-protect p))))
682 ;; Factor out content.
683 ((and (not square-free)
684 (not (eq 1 (setq t1 (math-mul (math-poly-gcd-list p)
685 (if (math-guess-if-neg
686 (nth (1- (length p)) p))
687 -1 1))))))
688 (math-accum-factors t1 1 (math-factor-poly-coefs
689 (math-poly-div-list p t1) 'cont)))
691 ;; Check if linear in math-fet-x.
692 ((not (cdr (cdr p)))
693 (math-sort-terms
694 (math-add (math-factor-protect
695 (math-sort-terms
696 (math-factor-expr (car p))))
697 (math-mul math-fet-x (math-factor-protect
698 (math-sort-terms
699 (math-factor-expr (nth 1 p))))))))
701 ;; If symbolic coefficients, use FactorRules.
702 ((let ((pp p))
703 (while (and pp (or (Math-ratp (car pp))
704 (and (eq (car (car pp)) 'mod)
705 (Math-integerp (nth 1 (car pp)))
706 (Math-integerp (nth 2 (car pp))))))
707 (setq pp (cdr pp)))
709 (let ((res (math-rewrite
710 (list 'calcFunc-thecoefs math-fet-x (cons 'vec p))
711 '(var FactorRules var-FactorRules))))
712 (or (and (eq (car-safe res) 'calcFunc-thefactors)
713 (= (length res) 3)
714 (math-vectorp (nth 2 res))
715 (let ((facs 1)
716 (vec (nth 2 res)))
717 (while (setq vec (cdr vec))
718 (setq facs (math-accum-factors (car vec) 1 facs)))
719 facs))
720 (math-build-polynomial-expr p math-fet-x))))
722 ;; Check if rational coefficients (i.e., not modulo a prime).
723 ((eq math-poly-modulus 1)
725 ;; Check if there are any squared terms, or a content not = 1.
726 (if (or (eq square-free t)
727 (equal (setq t1 (math-poly-gcd-coefs
728 p (setq t2 (math-poly-deriv-coefs p))))
729 '(1)))
731 ;; We now have a square-free polynomial with integer coefs.
732 ;; For now, we use a kludgy method that finds linear and
733 ;; quadratic terms using floating-point root-finding.
734 (if (setq t1 (let ((calc-symbolic-mode nil))
735 (math-poly-all-roots nil p t)))
736 (let ((roots (car t1))
737 (csign (if (math-negp (nth (1- (length p)) p)) -1 1))
738 (expr 1)
739 (unfac (nth 1 t1))
740 (scale (nth 2 t1)))
741 (while roots
742 (let ((coef0 (car (car roots)))
743 (coef1 (cdr (car roots))))
744 (setq expr (math-accum-factors
745 (if coef1
746 (let ((den (math-lcm-denoms
747 coef0 coef1)))
748 (setq scale (math-div scale den))
749 (math-add
750 (math-add
751 (math-mul den (math-pow math-fet-x 2))
752 (math-mul (math-mul coef1 den)
753 math-fet-x))
754 (math-mul coef0 den)))
755 (let ((den (math-lcm-denoms coef0)))
756 (setq scale (math-div scale den))
757 (math-add (math-mul den math-fet-x)
758 (math-mul coef0 den))))
759 1 expr)
760 roots (cdr roots))))
761 (setq expr (math-accum-factors
762 expr 1
763 (math-mul csign
764 (math-build-polynomial-expr
765 (math-mul-list (nth 1 t1) scale)
766 math-fet-x)))))
767 (math-build-polynomial-expr p math-fet-x)) ; can't factor it.
769 ;; Separate out the squared terms (Knuth exercise 4.6.2-34).
770 ;; This step also divides out the content of the polynomial.
771 (let* ((cabs (math-poly-gcd-list p))
772 (csign (if (math-negp (nth (1- (length p)) p)) -1 1))
773 (t1s (math-mul-list t1 csign))
774 (uu nil)
775 (v (car (math-poly-div-coefs p t1s)))
776 (w (car (math-poly-div-coefs t2 t1s))))
777 (while
778 (not (math-poly-zerop
779 (setq t2 (math-poly-simplify
780 (math-poly-mix
781 w 1 (math-poly-deriv-coefs v) -1)))))
782 (setq t1 (math-poly-gcd-coefs v t2)
783 uu (cons t1 uu)
784 v (car (math-poly-div-coefs v t1))
785 w (car (math-poly-div-coefs t2 t1))))
786 (setq t1 (length uu)
787 t2 (math-accum-factors (math-factor-poly-coefs v t)
788 (1+ t1) 1))
789 (while uu
790 (setq t2 (math-accum-factors (math-factor-poly-coefs
791 (car uu) t)
792 t1 t2)
793 t1 (1- t1)
794 uu (cdr uu)))
795 (math-accum-factors (math-mul cabs csign) 1 t2))))
797 ;; Factoring modulo a prime.
798 ((and (= (length (setq temp (math-poly-gcd-coefs
799 p (math-poly-deriv-coefs p))))
800 (length p)))
801 (setq p (car temp))
802 (while (cdr temp)
803 (setq temp (nthcdr (nth 2 math-poly-modulus) temp)
804 p (cons (car temp) p)))
805 (and (setq temp (math-factor-poly-coefs p))
806 (math-pow temp (nth 2 math-poly-modulus))))
808 (math-reject-arg nil "*Modulo factorization not yet implemented")))))
810 (defun math-poly-deriv-coefs (p)
811 (let ((n 1)
812 (dp nil))
813 (while (setq p (cdr p))
814 (setq dp (cons (math-mul (car p) n) dp)
815 n (1+ n)))
816 (nreverse dp)))
818 (defun math-factor-contains (x a)
819 (if (equal x a)
821 (if (memq (car-safe x) '(+ - * / neg))
822 (let ((sum 0))
823 (while (setq x (cdr x))
824 (setq sum (+ sum (math-factor-contains (car x) a))))
825 sum)
826 (if (and (eq (car-safe x) '^)
827 (natnump (nth 2 x)))
828 (* (math-factor-contains (nth 1 x) a) (nth 2 x))
829 0))))
835 ;;; Merge all quotients and expand/simplify the numerator
836 (defun calcFunc-nrat (expr)
837 (if (math-any-floats expr)
838 (setq expr (calcFunc-pfrac expr)))
839 (if (or (math-vectorp expr)
840 (assq (car-safe expr) calc-tweak-eqn-table))
841 (cons (car expr) (mapcar 'calcFunc-nrat (cdr expr)))
842 (let* ((calc-prefer-frac t)
843 (res (math-to-ratpoly expr))
844 (num (math-simplify (math-sort-terms (calcFunc-expand (car res)))))
845 (den (math-simplify (math-sort-terms (calcFunc-expand (cdr res)))))
846 (g (math-poly-gcd num den)))
847 (or (eq g 1)
848 (let ((num2 (math-poly-div num g))
849 (den2 (math-poly-div den g)))
850 (and (eq (cdr num2) 0) (eq (cdr den2) 0)
851 (setq num (car num2) den (car den2)))))
852 (math-simplify (math-div num den)))))
854 ;;; Returns expressions (num . denom).
855 (defun math-to-ratpoly (expr)
856 (let ((res (math-to-ratpoly-rec expr)))
857 (cons (math-simplify (car res)) (math-simplify (cdr res)))))
859 (defun math-to-ratpoly-rec (expr)
860 (cond ((Math-primp expr)
861 (cons expr 1))
862 ((memq (car expr) '(+ -))
863 (let ((r1 (math-to-ratpoly-rec (nth 1 expr)))
864 (r2 (math-to-ratpoly-rec (nth 2 expr))))
865 (if (equal (cdr r1) (cdr r2))
866 (cons (list (car expr) (car r1) (car r2)) (cdr r1))
867 (if (eq (cdr r1) 1)
868 (cons (list (car expr)
869 (math-mul (car r1) (cdr r2))
870 (car r2))
871 (cdr r2))
872 (if (eq (cdr r2) 1)
873 (cons (list (car expr)
874 (car r1)
875 (math-mul (car r2) (cdr r1)))
876 (cdr r1))
877 (let ((g (math-poly-gcd (cdr r1) (cdr r2))))
878 (let ((d1 (and (not (eq g 1)) (math-poly-div (cdr r1) g)))
879 (d2 (and (not (eq g 1)) (math-poly-div
880 (math-mul (car r1) (cdr r2))
881 g))))
882 (if (and (eq (cdr d1) 0) (eq (cdr d2) 0))
883 (cons (list (car expr) (car d2)
884 (math-mul (car r2) (car d1)))
885 (math-mul (car d1) (cdr r2)))
886 (cons (list (car expr)
887 (math-mul (car r1) (cdr r2))
888 (math-mul (car r2) (cdr r1)))
889 (math-mul (cdr r1) (cdr r2)))))))))))
890 ((eq (car expr) '*)
891 (let* ((r1 (math-to-ratpoly-rec (nth 1 expr)))
892 (r2 (math-to-ratpoly-rec (nth 2 expr)))
893 (g (math-mul (math-poly-gcd (car r1) (cdr r2))
894 (math-poly-gcd (cdr r1) (car r2)))))
895 (if (eq g 1)
896 (cons (math-mul (car r1) (car r2))
897 (math-mul (cdr r1) (cdr r2)))
898 (cons (math-poly-div-exact (math-mul (car r1) (car r2)) g)
899 (math-poly-div-exact (math-mul (cdr r1) (cdr r2)) g)))))
900 ((eq (car expr) '/)
901 (let* ((r1 (math-to-ratpoly-rec (nth 1 expr)))
902 (r2 (math-to-ratpoly-rec (nth 2 expr))))
903 (if (and (eq (cdr r1) 1) (eq (cdr r2) 1))
904 (cons (car r1) (car r2))
905 (let ((g (math-mul (math-poly-gcd (car r1) (car r2))
906 (math-poly-gcd (cdr r1) (cdr r2)))))
907 (if (eq g 1)
908 (cons (math-mul (car r1) (cdr r2))
909 (math-mul (cdr r1) (car r2)))
910 (cons (math-poly-div-exact (math-mul (car r1) (cdr r2)) g)
911 (math-poly-div-exact (math-mul (cdr r1) (car r2))
912 g)))))))
913 ((and (eq (car expr) '^) (integerp (nth 2 expr)))
914 (let ((r1 (math-to-ratpoly-rec (nth 1 expr))))
915 (if (> (nth 2 expr) 0)
916 (cons (math-pow (car r1) (nth 2 expr))
917 (math-pow (cdr r1) (nth 2 expr)))
918 (cons (math-pow (cdr r1) (- (nth 2 expr)))
919 (math-pow (car r1) (- (nth 2 expr)))))))
920 ((eq (car expr) 'neg)
921 (let ((r1 (math-to-ratpoly-rec (nth 1 expr))))
922 (cons (math-neg (car r1)) (cdr r1))))
923 (t (cons expr 1))))
926 (defun math-ratpoly-p (expr &optional var)
927 (cond ((equal expr var) 1)
928 ((Math-primp expr) 0)
929 ((memq (car expr) '(+ -))
930 (let ((p1 (math-ratpoly-p (nth 1 expr) var))
932 (and p1 (setq p2 (math-ratpoly-p (nth 2 expr) var))
933 (max p1 p2))))
934 ((eq (car expr) '*)
935 (let ((p1 (math-ratpoly-p (nth 1 expr) var))
937 (and p1 (setq p2 (math-ratpoly-p (nth 2 expr) var))
938 (+ p1 p2))))
939 ((eq (car expr) 'neg)
940 (math-ratpoly-p (nth 1 expr) var))
941 ((eq (car expr) '/)
942 (let ((p1 (math-ratpoly-p (nth 1 expr) var))
944 (and p1 (setq p2 (math-ratpoly-p (nth 2 expr) var))
945 (- p1 p2))))
946 ((and (eq (car expr) '^)
947 (integerp (nth 2 expr)))
948 (let ((p1 (math-ratpoly-p (nth 1 expr) var)))
949 (and p1 (* p1 (nth 2 expr)))))
950 ((not var) 1)
951 ((math-poly-depends expr var) nil)
952 (t 0)))
955 (defun calcFunc-apart (expr &optional var)
956 (cond ((Math-primp expr) expr)
957 ((eq (car expr) '+)
958 (math-add (calcFunc-apart (nth 1 expr) var)
959 (calcFunc-apart (nth 2 expr) var)))
960 ((eq (car expr) '-)
961 (math-sub (calcFunc-apart (nth 1 expr) var)
962 (calcFunc-apart (nth 2 expr) var)))
963 ((and var (not (math-ratpoly-p expr var)))
964 (math-reject-arg expr "Expected a rational function"))
966 (let* ((calc-prefer-frac t)
967 (rat (math-to-ratpoly expr))
968 (num (car rat))
969 (den (cdr rat)))
970 (or var
971 (setq var (math-polynomial-base den)))
972 (if (not (math-ratpoly-p expr var))
973 (math-reject-arg expr "Expected a rational function")
974 (let* ((qr (math-poly-div num den))
975 (q (car qr))
976 (r (cdr qr)))
977 (math-add q (or (and var
978 (math-expr-contains den var)
979 (math-partial-fractions r den var))
980 (math-div r den)))))))))
983 (defun math-padded-polynomial (expr var deg)
984 "Return a polynomial as list of coefficients.
985 If EXPR is of the form \"a + bx + cx^2 + ...\" in the variable VAR, return
986 the list (a b c ...) with at least DEG elements, else return NIL."
987 (let ((p (math-is-polynomial expr var deg)))
988 (append p (make-list (- deg (length p)) 0))))
990 (defun math-partial-fractions (r den var)
991 "Return R divided by DEN expressed in partial fractions of VAR.
992 All whole factors of DEN have already been split off from R.
993 If no partial fraction representation can be found, return nil."
994 (let* ((fden (calcFunc-factors den var))
995 (tdeg (math-polynomial-p den var))
996 (fp fden)
997 (dlist nil)
998 (eqns 0)
999 (lz nil)
1000 (tz (make-list (1- tdeg) 0))
1001 (calc-matrix-mode 'scalar))
1002 (and (not (and (= (length fden) 2) (eq (nth 2 (nth 1 fden)) 1)))
1003 (progn
1004 (while (setq fp (cdr fp))
1005 (let ((rpt (nth 2 (car fp)))
1006 (deg (math-polynomial-p (nth 1 (car fp)) var))
1007 dnum dvar deg2)
1008 (while (> rpt 0)
1009 (setq deg2 deg
1010 dnum 0)
1011 (while (> deg2 0)
1012 (setq dvar (append '(vec) lz '(1) tz)
1013 lz (cons 0 lz)
1014 tz (cdr tz)
1015 deg2 (1- deg2)
1016 dnum (math-add dnum (math-mul dvar
1017 (math-pow var deg2)))
1018 dlist (cons (and (= deg2 (1- deg))
1019 (math-pow (nth 1 (car fp)) rpt))
1020 dlist)))
1021 (let ((fpp fden)
1022 (mult 1))
1023 (while (setq fpp (cdr fpp))
1024 (or (eq fpp fp)
1025 (setq mult (math-mul mult
1026 (math-pow (nth 1 (car fpp))
1027 (nth 2 (car fpp)))))))
1028 (setq dnum (math-mul dnum mult)))
1029 (setq eqns (math-add eqns (math-mul dnum
1030 (math-pow
1031 (nth 1 (car fp))
1032 (- (nth 2 (car fp))
1033 rpt))))
1034 rpt (1- rpt)))))
1035 (setq eqns (math-div (cons 'vec (math-padded-polynomial r var tdeg))
1036 (math-transpose
1037 (cons 'vec
1038 (mapcar
1039 (function
1040 (lambda (x)
1041 (cons 'vec (math-padded-polynomial
1042 x var tdeg))))
1043 (cdr eqns))))))
1044 (and (math-vectorp eqns)
1045 (let ((res 0)
1046 (num nil))
1047 (setq eqns (nreverse eqns))
1048 (while eqns
1049 (setq num (cons (car eqns) num)
1050 eqns (cdr eqns))
1051 (if (car dlist)
1052 (setq num (math-build-polynomial-expr
1053 (nreverse num) var)
1054 res (math-add res (math-div num (car dlist)))
1055 num nil))
1056 (setq dlist (cdr dlist)))
1057 (math-normalize res)))))))
1061 (defun math-expand-term (expr)
1062 (cond ((and (eq (car-safe expr) '*)
1063 (memq (car-safe (nth 1 expr)) '(+ -)))
1064 (math-add-or-sub (list '* (nth 1 (nth 1 expr)) (nth 2 expr))
1065 (list '* (nth 2 (nth 1 expr)) (nth 2 expr))
1066 nil (eq (car (nth 1 expr)) '-)))
1067 ((and (eq (car-safe expr) '*)
1068 (memq (car-safe (nth 2 expr)) '(+ -)))
1069 (math-add-or-sub (list '* (nth 1 expr) (nth 1 (nth 2 expr)))
1070 (list '* (nth 1 expr) (nth 2 (nth 2 expr)))
1071 nil (eq (car (nth 2 expr)) '-)))
1072 ((and (eq (car-safe expr) '/)
1073 (memq (car-safe (nth 1 expr)) '(+ -)))
1074 (math-add-or-sub (list '/ (nth 1 (nth 1 expr)) (nth 2 expr))
1075 (list '/ (nth 2 (nth 1 expr)) (nth 2 expr))
1076 nil (eq (car (nth 1 expr)) '-)))
1077 ((and (eq (car-safe expr) '^)
1078 (memq (car-safe (nth 1 expr)) '(+ -))
1079 (integerp (nth 2 expr))
1080 (if (and
1081 (or (math-known-matrixp (nth 1 (nth 1 expr)))
1082 (math-known-matrixp (nth 2 (nth 1 expr)))
1083 (and
1084 calc-matrix-mode
1085 (not (eq calc-matrix-mode 'scalar))
1086 (not (and (math-known-scalarp (nth 1 (nth 1 expr)))
1087 (math-known-scalarp (nth 2 (nth 1 expr)))))))
1088 (> (nth 2 expr) 1))
1089 (if (= (nth 2 expr) 2)
1090 (math-add-or-sub (list '* (nth 1 (nth 1 expr)) (nth 1 expr))
1091 (list '* (nth 2 (nth 1 expr)) (nth 1 expr))
1092 nil (eq (car (nth 1 expr)) '-))
1093 (math-add-or-sub (list '* (nth 1 (nth 1 expr))
1094 (list '^ (nth 1 expr)
1095 (1- (nth 2 expr))))
1096 (list '* (nth 2 (nth 1 expr))
1097 (list '^ (nth 1 expr)
1098 (1- (nth 2 expr))))
1099 nil (eq (car (nth 1 expr)) '-)))
1100 (if (> (nth 2 expr) 0)
1101 (or (and (or (> math-mt-many 500000) (< math-mt-many -500000))
1102 (math-expand-power (nth 1 expr) (nth 2 expr)
1103 nil t))
1104 (list '*
1105 (nth 1 expr)
1106 (list '^ (nth 1 expr) (1- (nth 2 expr)))))
1107 (if (< (nth 2 expr) 0)
1108 (list '/ 1 (list '^ (nth 1 expr) (- (nth 2 expr)))))))))
1109 (t expr)))
1111 (defun calcFunc-expand (expr &optional many)
1112 (math-normalize (math-map-tree 'math-expand-term expr many)))
1114 (defun math-expand-power (x n &optional var else-nil)
1115 (or (and (natnump n)
1116 (memq (car-safe x) '(+ -))
1117 (let ((terms nil)
1118 (cterms nil))
1119 (while (memq (car-safe x) '(+ -))
1120 (setq terms (cons (if (eq (car x) '-)
1121 (math-neg (nth 2 x))
1122 (nth 2 x))
1123 terms)
1124 x (nth 1 x)))
1125 (setq terms (cons x terms))
1126 (if var
1127 (let ((p terms))
1128 (while p
1129 (or (math-expr-contains (car p) var)
1130 (setq terms (delq (car p) terms)
1131 cterms (cons (car p) cterms)))
1132 (setq p (cdr p)))
1133 (if cterms
1134 (setq terms (cons (apply 'calcFunc-add cterms)
1135 terms)))))
1136 (if (= (length terms) 2)
1137 (let ((i 0)
1138 (accum 0))
1139 (while (<= i n)
1140 (setq accum (list '+ accum
1141 (list '* (calcFunc-choose n i)
1142 (list '*
1143 (list '^ (nth 1 terms) i)
1144 (list '^ (car terms)
1145 (- n i)))))
1146 i (1+ i)))
1147 accum)
1148 (if (= n 2)
1149 (let ((accum 0)
1150 (p1 terms)
1152 (while p1
1153 (setq accum (list '+ accum
1154 (list '^ (car p1) 2))
1155 p2 p1)
1156 (while (setq p2 (cdr p2))
1157 (setq accum (list '+ accum
1158 (list '* 2 (list '*
1159 (car p1)
1160 (car p2))))))
1161 (setq p1 (cdr p1)))
1162 accum)
1163 (if (= n 3)
1164 (let ((accum 0)
1165 (p1 terms)
1166 p2 p3)
1167 (while p1
1168 (setq accum (list '+ accum (list '^ (car p1) 3))
1169 p2 p1)
1170 (while (setq p2 (cdr p2))
1171 (setq accum (list '+
1172 (list '+
1173 accum
1174 (list '* 3
1175 (list
1177 (list '^ (car p1) 2)
1178 (car p2))))
1179 (list '* 3
1180 (list
1181 '* (car p1)
1182 (list '^ (car p2) 2))))
1183 p3 p2)
1184 (while (setq p3 (cdr p3))
1185 (setq accum (list '+ accum
1186 (list '* 6
1187 (list '*
1188 (car p1)
1189 (list
1190 '* (car p2)
1191 (car p3))))))))
1192 (setq p1 (cdr p1)))
1193 accum))))))
1194 (and (not else-nil)
1195 (list '^ x n))))
1197 (defun calcFunc-expandpow (x n)
1198 (math-normalize (math-expand-power x n)))
1200 (provide 'calc-poly)
1202 ;;; calc-poly.el ends here