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