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

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