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

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