cgit.sxemacs.org Git - packages/blob - xemacs-packages/calc/calc-poly.el

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