code.delx.au - gnu-emacs/blob - lisp/calc/calc-poly.el

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