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

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