1;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;; 2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3;;; The data in this file contains enhancments. ;;;;; 4;;; ;;;;; 5;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;; 6;;; All rights reserved ;;;;; 7;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 8;;; (c) Copyright 1982 Massachusetts Institute of Technology ;;; 9;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 10 11(in-package :maxima) 12 13(macsyma-module asum) 14 15(load-macsyma-macros rzmac) 16 17(declare-top (special opers *a *n $factlim sum msump *i *opers-list opers-list $ratsimpexpons makef $factorial_expand)) 18 19(loop for (x y) on '(%cot %tan %csc %sin %sec %cos %coth %tanh %csch %sinh %sech %cosh) 20 by #'cddr do (putprop x y 'recip) (putprop y x 'recip)) 21 22(defmvar $zeta%pi t) 23 24;; polynomial predicates and other such things 25 26(defun poly? (exp var) 27 (cond ((or (atom exp) (free exp var))) 28 ((member (caar exp) '(mtimes mplus) :test #'eq) 29 (do ((exp (cdr exp) (cdr exp))) 30 ((null exp) t) 31 (and (null (poly? (car exp) var)) (return nil)))) 32 ((and (eq (caar exp) 'mexpt) 33 (integerp (caddr exp)) 34 (> (caddr exp) 0)) 35 (poly? (cadr exp) var)))) 36 37(defun smono (x var) 38 (smonogen x var t)) 39 40(defun smonop (x var) 41 (smonogen x var nil)) 42 43(defun smonogen (x var fl) ; fl indicates whether to return *a *n 44 (cond ((free x var) (and fl (setq *n 0 *a x)) t) 45 ((atom x) (and fl (setq *n (setq *a 1))) t) 46 ((and (listp (car x)) 47 (eq (caar x) 'mtimes)) 48 (do ((x (cdr x) (cdr x)) 49 (a '(1)) (n '(0))) 50 ((null x) 51 (and fl (setq *n (addn n nil) *a (muln a nil))) t) 52 (let (*a *n) 53 (if (smonogen (car x) var fl) 54 (and fl (setq a (cons *a a) n (cons *n n))) 55 (return nil))))) 56 ((and (listp (car x)) 57 (eq (caar x) 'mexpt)) 58 (cond ((and (free (caddr x) var) (eq (cadr x) var)) 59 (and fl (setq *n (caddr x) *a 1)) t))))) 60 61;; factorial stuff 62 63(defmvar $factlim 100000) ; set to a big integer which will work (not -1) 64(defvar makef nil) 65 66(defmfun $genfact (&rest l) 67 (cons '(%genfact) l)) 68 69(defun gfact (n %m i) 70 (cond ((minusp %m) (improper-arg-err %m '$genfact)) 71 ((= %m 0) 1) 72 (t (prog (ans) 73 (setq ans n) 74 a (if (= %m 1) (return ans)) 75 (setq n (m- n i) %m (1- %m) ans (m* ans n)) 76 (go a))))) 77 78;; From Richard Fateman's paper, "Comments on Factorial Programs", 79;; http://www.cs.berkeley.edu/~fateman/papers/factorial.pdf 80;; 81;; k(n,m) = n*(n-m)*(n-2*m)*... 82;; 83;; (k n 1) is n! 84;; 85;; This is much faster (3-4 times) than the original factorial 86;; function. 87 88(defun k (n m) 89 (if (<= n m) 90 n 91 (* (k n (* 2 m)) 92 (k (- n m) (* 2 m))))) 93 94(defun factorial (n) 95 (if (zerop n) 96 1 97 (k n 1))) 98 99;;; Factorial has mirror symmetry 100 101(defprop mfactorial t commutes-with-conjugate) 102 103(defun simpfact (x y z) 104 (oneargcheck x) 105 (setq y (simpcheck (cadr x) z)) 106 (cond ((and (mnump y) 107 (eq ($sign y) '$neg) 108 (zerop1 (sub (simplify (list '(%truncate) y)) y))) 109 ;; Negative integer or a real representation of a negative integer. 110 (merror (intl:gettext "factorial: factorial of negative integer ~:M not defined.") y)) 111 ((or (floatp y) 112 ($bfloatp y) 113 (and (not (integerp y)) 114 (not (ratnump y)) 115 (or (and (complex-number-p y 'float-or-rational-p) 116 (or $numer 117 (floatp ($realpart y)) 118 (floatp ($imagpart y)))) 119 (and (complex-number-p y 'bigfloat-or-number-p) 120 (or $numer 121 ($bfloatp ($realpart y)) 122 ($bfloatp ($imagpart y)))))) 123 (and (not makef) (ratnump y) (equal (caddr y) 2))) 124 ;; Numerically evaluate for real or complex argument in float or 125 ;; bigfloat precision using the Gamma function 126 (simplify (list '(%gamma) (add 1 y)))) 127 ((eq y '$inf) '$inf) 128 ((and $factorial_expand 129 (mplusp y) 130 (integerp (cadr y))) 131 ;; factorial(n+m) and m integer. Expand. 132 (let ((m (cadr y)) 133 (n (simplify (cons '(mplus) (cddr y))))) 134 (cond ((>= m 0) 135 (mul 136 (simplify (list '($pochhammer) (add n 1) m)) 137 (simplify (list '(mfactorial) n)))) 138 ((< m 0) 139 (setq m (- m)) 140 (div 141 (mul (power -1 m) (simplify (list '(mfactorial) n))) 142 ;; We factor to get the ordering (n-1)*(n-2)*... 143 ($factor 144 (simplify (list '($pochhammer) (mul -1 n) m)))))))) 145 ((or (not (fixnump y)) (not (> y -1))) 146 (eqtest (list '(mfactorial) y) x)) 147 ((or (minusp $factlim) (not (> y $factlim))) 148 (factorial y)) 149 (t (eqtest (list '(mfactorial) y) x)))) 150 151(defun makegamma1 (e) 152 (cond ((atom e) e) 153 ((eq (caar e) 'mfactorial) 154 (list '(%gamma) (list '(mplus) 1 (makegamma1 (cadr e))))) 155 156 ;; Begin code copied from orthopoly/orthopoly-init.lisp 157 ;; Do pochhammer(x,n) ==> gamma(x+n)/gamma(x). 158 159 ((eq (caar e) '$pochhammer) 160 (let ((x (makegamma1 (nth 1 e))) 161 (n (makegamma1 (nth 2 e)))) 162 (div (take '(%gamma) (add x n)) (take '(%gamma) x)))) 163 164 ;; (gamma(x/z+1)*z^floor(y))/gamma(x/z-floor(y)+1) 165 166 ((eq (caar e) '%genfact) 167 (let ((x (makegamma1 (nth 1 e))) 168 (y (makegamma1 (nth 2 e))) 169 (z (makegamma1 (nth 3 e)))) 170 (setq y (take '($floor) y)) 171 (div 172 (mul 173 (take '(%gamma) (add (div x z) 1)) 174 (power z y)) 175 (take '(%gamma) (sub (add (div x z) 1) y))))) 176 ;; End code copied from orthopoly/orthopoly-init.lisp 177 178 ;; Double factorial 179 180 ((eq (caar e) '%double_factorial) 181 (let ((x (makegamma1 (nth 1 e)))) 182 (mul 183 (power 184 (div 2 '$%pi) 185 (mul 186 (div 1 4) 187 (sub 1 (simplify (list '(%cos) (mul '$%pi x)))))) 188 (power 2 (div x 2)) 189 (simplify (list '(%gamma) (add 1 (div x 2))))))) 190 191 ((eq (caar e) '%elliptic_kc) 192 ;; Complete elliptic integral of the first kind 193 (cond ((alike1 (cadr e) '((rat simp) 1 2)) 194 ;; K(1/2) = gamma(1/4)/4/sqrt(pi) 195 '((mtimes simp) ((rat simp) 1 4) 196 ((mexpt simp) $%pi ((rat simp) -1 2)) 197 ((mexpt simp) ((%gamma simp) ((rat simp) 1 4)) 2))) 198 ((or (alike1 (cadr e) 199 '((mtimes simp) ((rat simp) 1 4) 200 ((mplus simp) 2 201 ((mexpt simp) 3 ((rat simp) 1 2))))) 202 (alike1 (cadr e) 203 '((mplus simp) ((rat simp) 1 2) 204 ((mtimes simp) ((rat simp) 1 4) 205 ((mexpt simp) 3 ((rat simp) 1 2))))) 206 (alike1 (cadr e) 207 ;; 1/(8-4*sqrt(3)) 208 '((mexpt simp) 209 ((mplus simp) 8 210 ((mtimes simp) -4 211 ((mexpt simp) 3 ((rat simp) 1 2)))) 212 -1))) 213 ;; K((2+sqrt(3)/4)) 214 '((mtimes simp) ((rat simp) 1 4) 215 ((mexpt simp) 3 ((rat simp) 1 4)) 216 ((mexpt simp) $%pi ((rat simp) -1 2)) 217 ((%gamma simp) ((rat simp) 1 6)) 218 ((%gamma simp) ((rat simp) 1 3)))) 219 ((or (alike1 (cadr e) 220 ;; (2-sqrt(3))/4 221 '((mtimes simp) ((rat simp) 1 4) 222 ((mplus simp) 2 223 ((mtimes simp) -1 224 ((mexpt simp) 3 ((rat simp) 1 2)))))) 225 (alike1 (cadr e) 226 ;; 1/2-sqrt(3)/4 227 '((mplus simp) ((rat simp) 1 2) 228 ((mtimes simp) ((rat simp) -1 4) 229 ((mexpt simp) 3 ((rat simp) 1 2))))) 230 (alike (cadr e) 231 ;; 1/(4*sqrt(3)+8) 232 '((mexpt simp) 233 ((mplus simp) 8 234 ((mtimes simp) 4 235 ((mexpt simp) 3 ((rat simp) 1 2)))) 236 -1))) 237 ;; K((2-sqrt(3))/4) 238 '((mtimes simp) ((rat simp) 1 4) 239 ((mexpt simp) 3 ((rat simp) -1 4)) 240 ((mexpt simp) $%pi ((rat simp) -1 2)) 241 ((%gamma simp) ((rat simp) 1 6)) 242 ((%gamma simp) ((rat simp) 1 3)))) 243 ((or 244 ;; (3-2*sqrt(2))/(3+2*sqrt(2)) 245 (alike1 (cadr e) 246 '((mtimes simp) 247 ((mplus simp) 3 248 ((mtimes simp) -2 249 ((mexpt simp) 2 ((rat simp) 1 2)))) 250 ((mexpt simp) 251 ((mplus simp) 3 252 ((mtimes simp) 2 253 ((mexpt simp) 2 ((rat simp) 1 2)))) -1))) 254 ;; 17 - 12*sqrt(2) 255 (alike1 (cadr e) 256 '((mplus simp) 17 257 ((mtimes simp) -12 258 ((mexpt simp) 2 ((rat simp) 1 2))))) 259 ;; (2*SQRT(2) - 3)/(2*SQRT(2) + 3) 260 (alike1 (cadr e) 261 '((mtimes simp) -1 262 ((mplus simp) -3 263 ((mtimes simp) 2 264 ((mexpt simp) 2 ((rat simp) 1 2)))) 265 ((mexpt simp) 266 ((mplus simp) 3 267 ((mtimes simp) 2 268 ((mexpt simp) 2 ((rat simp) 1 2)))) 269 -1)))) 270 '((mtimes simp) ((rat simp) 1 8) 271 ((mexpt simp) 2 ((rat simp) -1 2)) 272 ((mplus simp) 1 ((mexpt simp) 2 ((rat simp) 1 2))) 273 ((mexpt simp) $%pi ((rat simp) -1 2)) 274 ((mexpt simp) ((%gamma simp) ((rat simp) 1 4)) 2))) 275 (t 276 ;; Give up 277 e))) 278 ((eq (caar e) '%elliptic_ec) 279 ;; Complete elliptic integral of the second kind 280 (cond ((alike1 (cadr e) '((rat simp) 1 2)) 281 ;; 2*E(1/2) - K(1/2) = 2*%pi^(3/2)*gamma(1/4)^(-2) 282 '((mplus simp) 283 ((mtimes simp) ((mexpt simp) $%pi ((rat simp) 3 2)) 284 ((mexpt simp) 285 ((%gamma simp irreducible) ((rat simp) 1 4)) -2)) 286 ((mtimes simp) ((rat simp) 1 8) 287 ((mexpt simp) $%pi ((rat simp) -1 2)) 288 ((mexpt simp) ((%gamma simp) ((rat simp) 1 4)) 2)))) 289 ((or (alike1 (cadr e) 290 '((mtimes simp) ((rat simp) 1 4) 291 ((mplus simp) 2 292 ((mtimes simp) -1 293 ((mexpt simp) 3 ((rat simp) 1 2)))))) 294 (alike1 (cadr e) 295 '((mplus simp) ((rat simp) 1 2) 296 ((mtimes simp) ((rat simp) -1 4) 297 ((mexpt simp) 3 ((rat simp) 1 2)))))) 298 ;; E((2-sqrt(3))/4) 299 ;; 300 ;; %pi/4/sqrt(3) = K*(E-(sqrt(3)+1)/2/sqrt(3)*K) 301 '((mplus simp) 302 ((mtimes simp) ((mexpt simp) 3 ((rat simp) -1 4)) 303 ((mexpt simp) $%pi ((rat simp) 3 2)) 304 ((mexpt simp) ((%gamma simp) ((rat simp) 1 6)) -1) 305 ((mexpt simp) ((%gamma simp) ((rat simp) 1 3)) -1)) 306 ((mtimes simp) ((rat simp) 1 8) 307 ((mexpt simp) 3 ((rat simp) -3 4)) 308 ((mexpt simp) $%pi ((rat simp) -1 2)) 309 ((%gamma simp) ((rat simp) 1 6)) 310 ((%gamma simp) ((rat simp) 1 3))) 311 ((mtimes simp) ((rat simp) 1 8) 312 ((mexpt simp) 3 ((rat simp) -1 4)) 313 ((mexpt simp) $%pi ((rat simp) -1 2)) 314 ((%gamma simp) ((rat simp) 1 6)) 315 ((%gamma simp) ((rat simp) 1 3))))) 316 ((or (alike1 (cadr e) 317 '((mtimes simp) ((rat simp) 1 4) 318 ((mplus simp) 2 319 ((mexpt simp) 3 ((rat simp) 1 2))))) 320 (alike1 (cadr e) 321 '((mplus simp) ((rat simp) 1 2) 322 ((mtimes simp) ((rat simp) 1 4) 323 ((mexpt simp) 3 ((rat simp) 1 2)))))) 324 ;; E((2+sqrt(3))/4) 325 ;; 326 ;; %pi*sqrt(3)/4 = K1*(E1-(sqrt(3)-1)/2/sqrt(3)*K1) 327 '((mplus simp) 328 ((mtimes simp) 3 ((mexpt simp) 3 ((rat simp) -3 4)) 329 ((mexpt simp) $%pi ((rat simp) 3 2)) 330 ((mexpt simp) ((%gamma simp) ((rat simp) 1 6)) -1) 331 ((mexpt simp) ((%gamma simp) ((rat simp) 1 3)) -1)) 332 ((mtimes simp) ((rat simp) 3 8) 333 ((mexpt simp) 3 ((rat simp) -3 4)) 334 ((mexpt simp) $%pi ((rat simp) -1 2)) 335 ((%gamma simp) ((rat simp) 1 6)) 336 ((%gamma simp) ((rat simp) 1 3))) 337 ((mtimes simp) ((rat simp) -1 8) 338 ((mexpt simp) 3 ((rat simp) -1 4)) 339 ((mexpt simp) $%pi ((rat simp) -1 2)) 340 ((%gamma simp) ((rat simp) 1 6)) 341 ((%gamma simp) ((rat simp) 1 3))))) 342 (t 343 e))) 344 (t (recur-apply #'makegamma1 e)))) 345 346(defun simpgfact (x vestigial z) 347 (declare (ignore vestigial)) 348 (arg-count-check 3 x) 349 (setq z (mapcar #'(lambda (q) (simpcheck q z)) (cdr x))) 350 (let ((a (car z)) (b (take '($floor) (cadr z))) (c (caddr z))) 351 (cond ((and (fixnump a) 352 (fixnump b) 353 (fixnump c)) 354 (if (and (> a -1) 355 (> b -1) 356 (or (<= c a) (= b 0)) 357 (<= b (/ a c))) 358 (gfact a b c) 359 (merror (intl:gettext "genfact: generalized factorial not defined for given arguments.")))) 360 (t (eqtest (list '(%genfact) a 361 (if (and (not (atom b)) 362 (eq (caar b) '$floor)) 363 (cadr b) 364 b) 365 c) 366 x))))) 367 368;; sum begins 369 370(defmvar $cauchysum nil 371 "When multiplying together sums with INF as their upper limit, 372causes the Cauchy product to be used rather than the usual product. 373In the Cauchy product the index of the inner summation is a function of 374the index of the outer one rather than varying independently." 375 modified-commands '$sum) 376 377(defmvar $gensumnum 0 378 "The numeric suffix used to generate the next variable of 379summation. If it is set to FALSE then the index will consist only of 380GENINDEX with no numeric suffix." 381 modified-commands '$sum 382 setting-predicate #'(lambda (x) (or (null x) (integerp x)))) 383 384(defmvar $genindex '$i 385 "The alphabetic prefix used to generate the next variable of 386summation when necessary." 387 modified-commands '$sum 388 setting-predicate #'symbolp) 389 390(defmvar $zerobern t) 391(defmvar $simpsum nil) 392(defmvar $simpproduct nil) 393 394(defvar *infsumsimp t) 395 396;; These variables should be initialized where they belong. 397 398(defmvar $cflength 1) 399(defmvar $taylordepth 3) 400(defmvar $maxtaydiff 4) 401(defmvar $verbose nil) 402(defvar *trunclist nil) 403(defvar ps-bmt-disrep t) 404(defvar silent-taylor-flag nil) 405 406(defmacro sum-arg (sum) 407 `(cadr ,sum)) 408 409(defmacro sum-index (sum) 410 `(caddr ,sum)) 411 412(defmacro sum-lower (sum) 413 `(cadddr ,sum)) 414 415(defmacro sum-upper (sum) 416 `(cadr (cdddr ,sum))) 417 418(defmspec $sum (l) 419 (arg-count-check 4 l) 420 (setq l (cdr l)) 421 (dosum (car l) (cadr l) (meval (caddr l)) (meval (cadddr l)) t :evaluate-summand t)) 422 423 424(defmspec $lsum (l) 425 (arg-count-check 3 l) 426 (setq l (cdr l)) 427 ;;(or (= (length l) 3) (wna-err '$lsum)) 428 (let ((form (car l)) 429 (ind (cadr l)) 430 (lis (meval (caddr l))) 431 (ans 0)) 432 (or (symbolp ind) (merror (intl:gettext "lsum: second argument must be a variable; found ~M") ind)) 433 (cond (($listp lis) 434 (loop for v in (cdr lis) 435 with lind = (cons ind nil) 436 for w = (cons v nil) 437 do 438 (setq ans (add* ans (mbinding (lind w) (meval form))))) 439 ans) 440 (t `((%lsum) ,form ,ind ,lis))))) 441 442(defun simpsum (x y z) 443 (let (($ratsimpexpons t)) 444 (setq y (simplifya (sum-arg x) z))) 445 (simpsum1 y (sum-index x) (simplifya (sum-lower x) z) 446 (simplifya (sum-upper x) z))) 447 448; This function was SIMPSUM1 until the sum/product code was revised Nov 2005. 449; The revised code punts back to this function since this code knows 450; some simplifications not handled by the revised code. -- Robert Dodier 451 452(defun simpsum1-save (exp i lo hi) 453 (cond ((not (symbolp i)) (merror (intl:gettext "sum: index must be a symbol; found ~M") i)) 454 ((equal lo hi) (mbinding ((list i) (list hi)) (meval exp))) 455 ((and (atom exp) 456 (not (eq exp i)) 457 (getl '%sum '($outative $linear))) 458 (freesum exp lo hi 1)) 459 ((null $simpsum) (list (get '%sum 'msimpind) exp i lo hi)) 460 ((and (or (eq lo '$minf) 461 (alike1 lo '((mtimes simp) -1 $inf))) 462 (equal hi '$inf)) 463 (let ((pos-part (simpsum2 exp i 0 '$inf)) 464 (neg-part (simpsum2 (maxima-substitute (m- i) i exp) i 1 '$inf))) 465 (cond 466 ((or (eq neg-part '$und) 467 (eq pos-part '$und)) 468 '$und) 469 ((eq pos-part '$inf) 470 (if (eq neg-part '$minf) '$und '$inf)) 471 ((eq pos-part '$minf) 472 (if (eq neg-part '$inf) '$und '$minf)) 473 ((or (eq neg-part '$inf) (eq neg-part '$minf)) 474 neg-part) 475 (t (m+ neg-part pos-part))))) 476 ((or (eq lo '$minf) 477 (alike1 lo '((mtimes simp) -1 '$inf))) 478 (simpsum2 (maxima-substitute (m- i) i exp) i (m- hi) '$inf)) 479 (t (simpsum2 exp i lo hi)))) 480 481;; DOSUM, MEVALSUMARG, DO%SUM -- general principles 482 483;; - evaluate the summand/productand 484;; - substitute a gensym for the index variable and make assertions (via assume) about the gensym index 485;; - return 0/1 for empty sum/product. sumhack/prodhack are ignored 486;; - distribute sum/product over mbags when listarith = true 487 488(defun dosum (expr ind low hi sump &key (evaluate-summand t)) 489 (setq low (ratdisrep low) hi (ratdisrep hi)) ;; UGH, GAG WITH ME A SPOON 490 (if (not (symbolp ind)) 491 (merror (intl:gettext "~:M: index must be a symbol; found ~M") (if sump '$sum '$product) ind)) 492 (unwind-protect 493 (prog (u *i lind l*i *hl) 494 (setq lind (cons ind nil)) 495 (cond 496 ((not (fixnump (setq *hl (mfuncall '$floor (m- hi low))))) 497 (if evaluate-summand (setq expr (mevalsumarg expr ind low hi))) 498 (return (cons (if sump '(%sum) '(%product)) 499 (list expr ind low hi)))) 500 ((signp l *hl) 501 (return (if sump 0 1)))) 502 (setq *i low l*i (list *i) u (if sump 0 1)) 503 lo (setq u 504 (if sump 505 (add u (resimplify (let* ((foo (mbinding (lind l*i) (meval expr))) 506 (bar (subst-if-not-freeof *i ind foo))) 507 bar))) 508 (mul u (resimplify (let* ((foo (mbinding (lind l*i) (meval expr))) 509 (bar (subst-if-not-freeof *i ind foo))) 510 bar))))) 511 (when (zerop *hl) (return u)) 512 (setq *hl (1- *hl)) 513 (setq *i (car (rplaca l*i (m+ *i 1)))) 514 (go lo)))) 515 516(defun subst-if-not-freeof (x y expr) 517 (if ($freeof y expr) 518 expr 519 (if (atom expr) 520 x 521 (let* ((args (cdr expr)) 522 (L (eval `(mapcar (lambda (a) (subst-if-not-freeof ',x ',y a)) ',args)))) 523 (cons (car expr) L))))) 524 525(defun mevalsumarg (expr ind low hi) 526 (if (let (($prederror nil)) 527 (eq (mevalp `((mlessp) ,hi ,low)) t)) 528 0) 529 530 (let ((gensym-ind (gensym))) 531 (if (apparently-integer low) 532 (meval `(($declare) ,gensym-ind $integer))) 533 (assume (list '(mgeqp) gensym-ind low)) 534 (if (not (eq hi '$inf)) 535 (assume (list '(mgeqp) hi gensym-ind))) 536 (let ((msump t) (foo) (summand)) 537 (setq summand 538 (if (and (not (atom expr)) (get (caar expr) 'mevalsumarg-macro)) 539 (funcall (get (caar expr) 'mevalsumarg-macro) expr) 540 expr)) 541 (let (($simp nil)) 542 (setq summand ($substitute gensym-ind ind summand))) 543 (setq foo (mbinding ((list gensym-ind) (list gensym-ind)) 544 (resimplify (meval summand)))) 545 ;; At this point we do not switch off simplification to preserve 546 ;; the achieved simplification of the summand (DK 02/2010). 547 (let (($simp t)) 548 (setq foo ($substitute ind gensym-ind foo))) 549 (if (not (eq hi '$inf)) 550 (forget (list '(mgeqp) hi gensym-ind))) 551 (forget (list '(mgeqp) gensym-ind low)) 552 (if (apparently-integer low) 553 (meval `(($remove) ,gensym-ind $integer))) 554 foo))) 555 556(defun apparently-integer (x) 557 (or ($integerp x) ($featurep x '$integer))) 558 559(defun do%sum (l op) 560 (if (not (= (length l) 4)) (wna-err op)) 561 (let ((ind (cadr l))) 562 (if (mquotep ind) (setq ind (cadr ind))) 563 (if (not (symbolp ind)) 564 (merror (intl:gettext "~:M: index must be a symbol; found ~M") op ind)) 565 (let ((low (caddr l)) 566 (hi (cadddr l))) 567 (list (mevalsumarg (car l) ind low hi) 568 ind (meval (caddr l)) (meval (cadddr l)))))) 569 570(defun simpsum1 (e k lo hi) 571 (with-new-context (context) 572 (let ((acc 0) (n) (sgn) ($prederror nil) (i (gensym)) (ex)) 573 (setq lo ($ratdisrep lo)) 574 (setq hi ($ratdisrep hi)) 575 576 (setq n ($limit (add 1 (sub hi lo)))) 577 (setq sgn ($sign n)) 578 579 (if (not (eq t (csign lo))) (mfuncall '$assume `((mgeqp) ,i ,lo))) 580 (if (not (eq t (csign hi))) (mfuncall '$assume `((mgeqp) ,hi ,i))) 581 582 (setq ex (subst i k e)) 583 (setq ex (subst i k ex)) 584 585 (setq acc 586 (cond ((and (eq n '$inf) ($freeof i ex)) 587 (setq sgn (csign ex)) 588 (cond ((eq sgn '$pos) '$inf) 589 ((eq sgn '$neg) '$minf) 590 ((eq sgn '$zero) 0) 591 (t `((%sum simp) ,ex ,i ,lo ,hi)))) 592 593 ((and (mbagp e) $listarith) 594 (simplifya 595 `((,(caar e)) ,@(mapcar #'(lambda (s) (mfuncall '$sum s k lo hi)) (margs e))) t)) 596 597 ((or (eq sgn '$neg) (eq sgn '$zero) (eq sgn '$nz)) 0) 598 599 ((like ex 0) 0) 600 601 (($freeof i ex) (mult n ex)) 602 603 ((and (integerp n) (eq sgn '$pos) $simpsum) 604 (dotimes (j n acc) 605 (setq acc (add acc (resimplify (subst (add j lo) i ex)))))) 606 607 (t 608 (setq ex (subst '%sum '$sum ex)) 609 `((%sum simp) ,(subst k i ex) ,k ,lo ,hi)))) 610 611 (setq acc (subst k i acc)) 612 613 ;; If expression is still a summation, 614 ;; punt to previous simplification code. 615 616 (if (and $simpsum (op-equalp acc '$sum '%sum)) 617 (let* ((args (cdr acc)) (e (first args)) (i (second args)) (i0 (third args)) (i1 (fourth args))) 618 (setq acc (simpsum1-save e i i0 i1)))) 619 620 acc))) 621 622(defun simpprod1 (e k lo hi) 623 (with-new-context (context) 624 (let ((acc 1) (n) (sgn) ($prederror nil) (i (gensym)) (ex) (ex-mag) (realp)) 625 626 (setq lo ($ratdisrep lo)) 627 (setq hi ($ratdisrep hi)) 628 (setq n ($limit (add 1 (sub hi lo)))) 629 (setq sgn ($sign n)) 630 631 (if (not (eq t (csign lo))) (mfuncall '$assume `((mgeqp) ,i ,lo))) 632 (if (not (eq t (csign hi))) (mfuncall '$assume `((mgeqp) ,hi ,i))) 633 634 (setq ex (subst i k e)) 635 (setq ex (subst i k ex)) 636 637 (setq acc 638 (cond 639 ((like ex 1) 1) 640 641 ((and (eq n '$inf) ($freeof i ex)) 642 (setq ex-mag (mfuncall '$cabs ex)) 643 (setq realp (mfuncall '$imagpart ex)) 644 (setq realp (mevalp `((mequal) 0 ,realp))) 645 646 (cond ((eq t (mevalp `((mlessp) ,ex-mag 1))) 0) 647 ((and (eq realp t) (eq t (mevalp `((mgreaterp) ,ex 1)))) '$inf) 648 ((eq t (mevalp `((mgreaterp) ,ex-mag 1))) '$infinity) 649 ((eq t (mevalp `((mequal) 1 ,ex-mag))) '$und) 650 (t `((%product) ,e ,k ,lo ,hi)))) 651 652 ((or (eq sgn '$neg) (eq sgn '$zero) (eq sgn '$nz)) 653 1) 654 655 ((and (mbagp e) $listarith) 656 (simplifya 657 `((,(caar e)) ,@(mapcar #'(lambda (s) (mfuncall '$product s k lo hi)) (margs e))) t)) 658 659 (($freeof i ex) (power ex n)) 660 661 ((and (integerp n) (eq sgn '$pos) $simpproduct) 662 (dotimes (j n acc) 663 (setq acc (mult acc (resimplify (subst (add j lo) i ex)))))) 664 665 (t 666 (setq ex (subst '%product '$product ex)) 667 `((%product simp) ,(subst k i ex) ,k ,lo ,hi)))) 668 669 ;; Hmm, this is curious... don't call existing product simplifications 670 ;; if index range is infinite -- what's up with that?? 671 672 (if (and $simpproduct (op-equalp acc '$product '%product) (not (like n '$inf))) 673 (let* ((args (cdr acc)) (e (first args)) (i (second args)) (i0 (third args)) (i1 (fourth args))) 674 (setq acc (simpprod1-save e i i0 i1)))) 675 676 (setq acc (subst k i acc)) 677 (setq acc (subst '%product '$product acc)) 678 679 acc))) 680 681; This function was SIMPPROD1 until the sum/product code was revised Nov 2005. 682; The revised code punts back to this function since this code knows 683; some simplifications not handled by the revised code. -- Robert Dodier 684 685(defun simpprod1-save (exp i lo hi) 686 (let (u) 687 (cond ((not (symbolp i)) (merror (intl:gettext "product: index must be a symbol; found ~M") i)) 688 ((alike1 lo hi) 689 (let ((valist (list i))) 690 (mbinding (valist (list hi)) 691 (meval exp)))) 692 ((eq ($sign (setq u (m- hi lo))) '$neg) 693 (cond ((eq ($sign (add2 u 1)) '$zero) 1) 694 (t (merror (intl:gettext "product: lower bound ~M greater than upper bound ~M") lo hi)))) 695 ((atom exp) 696 (cond ((null (eq exp i)) 697 (power* exp (list '(mplus) hi 1 (list '(mtimes) -1 lo)))) 698 ((let ((lot (asksign lo))) 699 (cond ((equal lot '$zero) 0) 700 ((eq lot '$positive) 701 (m// (list '(mfactorial) hi) 702 (list '(mfactorial) (list '(mplus) lo -1)))) 703 ((m* (list '(mfactorial) 704 (list '(mabs) lo)) 705 (cond ((member (asksign hi) '($zero $positive) :test #'eq) 706 0) 707 (t (prog1 708 (m^ -1 (m+ hi lo 1)) 709 (setq hi (list '(mabs) hi))))) 710 (list '(mfactorial) hi)))))))) 711 ((list '(%product simp) exp i lo hi))))) 712 713 714;; multiplication of sums 715 716(defun gensumindex () 717 (intern (format nil "~S~D" $genindex (incf $gensumnum)))) 718 719(defun sumtimes (x y) 720 (cond ((null x) y) 721 ((null y) x) 722 ((or (safe-zerop x) (safe-zerop y)) 0) 723 ((or (atom x) (not (eq (caar x) '%sum))) (sumultin x y)) 724 ((or (atom y) (not (eq (caar y) '%sum))) (sumultin y x)) 725 (t (let (u v i j) 726 (if (great (sum-arg x) (sum-arg y)) (setq u y v x) (setq u x v y)) 727 (setq i (let ((ind (gensumindex))) 728 (setq u (subst ind (sum-index u) u)) ind)) 729 (setq j (let ((ind (gensumindex))) 730 (setq v (subst ind (sum-index v) v)) ind)) 731 (if (and $cauchysum (eq (sum-upper u) '$inf) 732 (eq (sum-upper v) '$inf)) 733 (list '(%sum) 734 (list '(%sum) 735 (sumtimes (maxima-substitute j i (sum-arg u)) 736 (maxima-substitute (m- i j) j (sum-arg v))) 737 j (sum-lower u) (m- i (sum-lower v))) 738 i (m+ (sum-lower u) (sum-lower v)) '$inf) 739 (list '(%sum) 740 (list '(%sum) (sumtimes (sum-arg u) (sum-arg v)) 741 j (sum-lower v) (sum-upper v)) 742 i (sum-lower u) (sum-upper u))))))) 743 744(defun sumultin (x s) ; Multiplies x into a sum adjusting indices. 745 (cond ((or (atom s) (not (eq (caar s) '%sum))) (m* x s)) 746 ((free x (sum-index s)) 747 (list* (car s) (sumultin x (sum-arg s)) (cddr s))) 748 (t (let ((ind (gensumindex))) 749 (list* (car s) 750 (sumultin x (subst ind (sum-index s) (sum-arg s))) 751 ind 752 (cdddr s)))))) 753 754;; addition of sums 755 756(defun sumpls (sum out) 757 (prog (l) 758 (if (null out) (return (cons sum nil))) 759 (setq out (setq l (cons nil out))) 760 a (if (null (cdr out)) (return (cons sum (cdr l)))) 761 (and (not (atom (cadr out))) 762 (consp (caadr out)) 763 (eq (caar (cadr out)) '%sum) 764 (alike1 (sum-arg (cadr out)) (sum-arg sum)) 765 (alike1 (sum-index (cadr out)) (sum-index sum)) 766 (cond ((onediff (sum-upper (cadr out)) (sum-lower sum)) 767 (setq sum (list (car sum) 768 (sum-arg sum) 769 (sum-index sum) 770 (sum-lower (cadr out)) 771 (sum-upper sum))) 772 (rplacd out (cddr out)) 773 (go a)) 774 ((onediff (sum-upper sum) (sum-lower (cadr out))) 775 (setq sum (list (car sum) 776 (sum-arg sum) 777 (sum-index sum) 778 (sum-lower sum) 779 (sum-upper (cadr out)))) 780 (rplacd out (cddr out)) 781 (go a)))) 782 (setq out (cdr out)) 783 (go a))) 784 785(defun onediff (x y) 786 (equal 1 (m- y x))) 787 788(defun freesum (e b a q) 789 (m* e q (m- (m+ a 1) b))) 790 791;; linear operator stuff 792 793(defparameter *opers-list '(($linear . linearize1))) 794(defparameter opers (list '$linear)) 795 796(defun oper-apply (e z) 797 (cond ((null opers-list) 798 (let ((w (get (caar e) 'operators))) 799 (if w (funcall w e 1 z) (simpargs e z)))) 800 ((get (caar e) (caar opers-list)) 801 (let ((opers-list (cdr opers-list)) 802 (fun (cdar opers-list))) 803 (funcall fun e z))) 804 (t (let ((opers-list (cdr opers-list))) 805 (oper-apply e z))))) 806 807;; Define an operator simplification, the same as antisymmetric, commutative, linear, etc. 808;; Here OP = operator name, FN = function of 1 argument to carry out operator-specific simplification. 809;; 1. push operator name onto OPERS 810;; 2. update $OPPROPERTIES 811;; 3. push operator name and glue code onto *OPERS-LIST 812;; 4. declare operator name as a feature, so declare(..., <op>) is recognized 813 814(defmfun $define_opproperty (op fn) 815 (unless (symbolp op) 816 (merror "define_opproperty: first argument must be a symbol; found: ~M" op)) 817 (unless (or (symbolp fn) (and (consp fn) (eq (caar fn) 'lambda))) 818 (merror "define_opproperty: second argument must be a symbol or lambda expression; found: ~M" fn)) 819 (push op opers) 820 (setq $opproperties (cons '(mlist simp) (reverse opers))) 821 (let 822 ((fn-glue (coerce (if (symbolp fn) 823 `(lambda (e z) 824 (declare (ignorable z)) 825 (if (or (fboundp ',fn) (mget ',fn 'mexpr)) 826 (let ((e1 (let ($simp) (mfuncall ',fn e)))) 827 (if ($mapatom e1) e1 (oper-apply e1 nil))) 828 (list '(,fn) (let ((*opers-list (cdr *opers-list))) (oper-apply e z))))) 829 `(lambda (e z) 830 (declare (ignore z)) 831 (let ((e1 (let ($simp) (mfuncall ',fn e)))) 832 (if ($mapatom e1) e1 (oper-apply e1 nil))))) 833 'function))) 834 (push `(,op . ,fn-glue) *opers-list)) 835 (mfuncall '$declare op '$feature)) 836 837(defun linearize1 (e z) ; z = t means args already simplified. 838 (linearize2 (cons (car e) (mapcar #'(lambda (q) (simpcheck q z)) (cdr e))) 839 nil)) 840 841(defun opident (op) 842 (cond ((eq op 'mplus) 0) 843 ((eq op 'mtimes) 1))) 844 845(defun rem-const (e) ;removes constantp stuff 846 (do ((l (cdr e) (cdr l)) 847 (a (list (opident (caar e)))) 848 (b (list (opident (caar e))))) 849 ((null l) 850 (cons (simplifya (cons (list (caar e)) a) nil) 851 (simplifya (cons (list (caar e)) b) nil))) 852 (if ($constantp (car l)) 853 (setq a (cons (car l) a)) 854 (setq b (cons (car l) b))))) 855 856(defun linearize2 (e times) 857 (cond ((linearconst e)) 858 ((atom (cadr e)) (oper-apply e t)) 859 ((eq (caar (cadr e)) 'mplus) 860 (addn (mapcar #'(lambda (q) 861 (linearize2 (list* (car e) q (cddr e)) nil)) 862 (cdr (cadr e))) 863 t)) 864 ((and (eq (caar (cadr e)) 'mtimes) (null times)) 865 (let ((z (if (and (cddr e) 866 (or (atom (caddr e)) 867 ($subvarp (caddr e)))) 868 (partition (cadr e) (caddr e) 1) 869 (rem-const (cadr e)))) 870 (w)) 871 (setq w (linearize2 (list* (car e) 872 (simplifya (cdr z) t) 873 (cddr e)) 874 t)) 875 (linearize3 w e (car z)))) 876 (t (oper-apply e t)))) 877 878(defun linearconst (e) 879 (if (or (mnump (cadr e)) 880 (constant (cadr e)) 881 (and (cddr e) 882 (or (atom (caddr e)) (member 'array (cdar (caddr e)) :test #'eq)) 883 (free (cadr e) (caddr e)))) 884 (if (or (zerop1 (cadr e)) 885 (and (member (caar e) '(%sum %integrate) :test #'eq) 886 (= (length e) 5) 887 (or (eq (cadddr e) '$minf) 888 (member (car (cddddr e)) '($inf $infinity) :test #'eq)) 889 (eq ($asksign (cadr e)) '$zero))) 890 0 891 (let ((w (oper-apply (list* (car e) 1 (cddr e)) t))) 892 (linearize3 w e (cadr e)))))) 893 894(defun linearize3 (w e x) 895 (let (w1) 896 (if (and (member w '($inf $minf $infinity) :test #'eq) (safe-zerop x)) 897 (merror (intl:gettext "LINEARIZE3: undefined form 0*inf: ~M") e)) 898 (setq w (mul2 (simplifya x t) w)) 899 (cond ((or (atom w) (getl (caar w) '($outative $linear))) (setq w1 1)) 900 ((eq (caar w) 'mtimes) 901 (setq w1 (cons '(mtimes) nil)) 902 (do ((w2 (cdr w) (cdr w2))) 903 ((null w2) (setq w1 (nreverse w1))) 904 (if (or (atom (car w2)) 905 (not (getl (caaar w2) '($outative $linear)))) 906 (setq w1 (cons (car w2) w1))))) 907 (t (setq w1 w))) 908 (if (and (not (atom w1)) (or (among '$inf w1) (among '$minf w1))) 909 (infsimp w) 910 w))) 911 912(setq opers (cons '$additive opers) 913 *opers-list (cons '($additive . additive) *opers-list)) 914 915(defun rem-opers-p (p) 916 (cond ((eq (caar opers-list) p) 917 (setq opers-list (cdr p))) 918 ((do ((l opers-list (cdr l))) 919 ((null l)) 920 (if (eq (caar (cdr l)) p) 921 (return (rplacd l (cddr l)))))))) 922 923(defun additive (e z) 924 (cond ((get (caar e) '$outative) ; Really a linearize! 925 (setq opers-list (copy-list opers-list)) 926 (rem-opers-p '$outative) 927 (linearize1 e z)) 928 ((mplusp (cadr e)) 929 (addn (mapcar #'(lambda (q) 930 (let ((opers-list *opers-list)) 931 (oper-apply (list* (car e) q (cddr e)) z))) 932 (cdr (cadr e))) 933 z)) 934 (t (oper-apply e z)))) 935 936(setq opers (cons '$multiplicative opers) 937 *opers-list (cons '($multiplicative . multiplicative) *opers-list)) 938 939(defun multiplicative (e z) 940 (cond ((mtimesp (cadr e)) 941 (muln (mapcar #'(lambda (q) 942 (let ((opers-list *opers-list)) 943 (oper-apply (list* (car e) q (cddr e)) z))) 944 (cdr (cadr e))) 945 z)) 946 (t (oper-apply e z)))) 947 948(setq opers (cons '$outative opers) 949 *opers-list (cons '($outative . outative) *opers-list)) 950 951(defun outative (e z) 952 (setq e (cons (car e) (mapcar #'(lambda (q) (simpcheck q z)) (cdr e)))) 953 (cond ((get (caar e) '$additive) 954 (setq opers-list (copy-list opers-list )) 955 (rem-opers-p '$additive) 956 (linearize1 e t)) 957 ((linearconst e)) 958 ((mtimesp (cadr e)) 959 (let ((u (if (and (cddr e) 960 (or (atom (caddr e)) 961 ($subvarp (caddr e)))) 962 (partition (cadr e) (caddr e) 1) 963 (rem-const (cadr e)))) 964 (w)) 965 (setq w (oper-apply (list* (car e) 966 (simplifya (cdr u) t) 967 (cddr e)) 968 t)) 969 (linearize3 w e (car u)))) 970 (t (oper-apply e t)))) 971 972(defprop %sum t $outative) 973(defprop %sum t opers) 974(defprop %integrate t $outative) 975(defprop %integrate t opers) 976(defprop %limit t $outative) 977(defprop %limit t opers) 978 979(setq opers (cons '$evenfun opers) 980 *opers-list (cons '($evenfun . evenfun) *opers-list)) 981 982(setq opers (cons '$oddfun opers) 983 *opers-list (cons '($oddfun . oddfun) *opers-list)) 984 985(defun evenfun (e z) 986 (if (or (null (cdr e)) (cddr e)) 987 (merror (intl:gettext "Function declared 'even' takes exactly one argument; found ~M") e)) 988 (let ((arg (simpcheck (cadr e) z))) 989 (oper-apply (list (car e) (if (mminusp arg) (neg arg) arg)) t))) 990 991(defun oddfun (e z) 992 (if (or (null (cdr e)) (cddr e)) 993 (merror (intl:gettext "Function declared 'odd' takes exactly one argument; found ~M") e)) 994 (let ((arg (simpcheck (cadr e) z))) 995 (if (mminusp arg) (neg (oper-apply (list (car e) (neg arg)) t)) 996 (oper-apply (list (car e) arg) t)))) 997 998(setq opers (cons '$commutative opers) 999 *opers-list (cons '($commutative . commutative1) *opers-list)) 1000 1001(setq opers (cons '$symmetric opers) 1002 *opers-list (cons '($symmetric . commutative1) *opers-list)) 1003 1004(defun commutative1 (e z) 1005 (oper-apply (cons (car e) 1006 (reverse 1007 (sort (mapcar #'(lambda (q) (simpcheck q z)) 1008 (cdr e)) 1009 'great))) 1010 t)) 1011 1012(setq opers (cons '$antisymmetric opers) 1013 *opers-list (cons '($antisymmetric . antisym) *opers-list)) 1014 1015(defun antisym (e z) 1016 (when (and $dotscrules (mnctimesp e)) 1017 (let ($dotexptsimp) 1018 (setq e (simpnct e 1 nil)))) 1019 (if ($atom e) e (antisym1 e z))) 1020 1021(defun antisym1 (e z) 1022 (let ((antisym-sign nil) 1023 (l (mapcar #'(lambda (q) (simpcheck q z)) (cdr e)))) 1024 (when (or (not (eq (caar e) 'mnctimes)) (freel l 'mnctimes)) 1025 (multiple-value-setq (l antisym-sign) (bbsort1 l))) 1026 (cond ((equal l 0) 0) 1027 ((prog1 1028 (null antisym-sign) 1029 (setq e (oper-apply (cons (car e) l) t))) 1030 e) 1031 (t (neg e))))) 1032 1033(defun bbsort1 (l) 1034 (prog (sl sl1 antisym-sign) 1035 (if (or (null l) (null (cdr l))) (return (values l antisym-sign)) 1036 (setq sl (list nil (car l)))) 1037 loop (setq l (cdr l)) 1038 (if (null l) (return (values (nreverse (cdr sl)) antisym-sign))) 1039 (setq sl1 sl) 1040 loop1(cond ((null (cdr sl1)) (rplacd sl1 (cons (car l) nil))) 1041 ((alike1 (car l) (cadr sl1)) (return (values 0 nil))) 1042 ((great (car l) (cadr sl1)) (rplacd sl1 (cons (car l) (cdr sl1)))) 1043 (t (setq antisym-sign (not antisym-sign) sl1 (cdr sl1)) (go loop1))) 1044 (go loop))) 1045 1046(setq opers (cons '$nary opers) 1047 *opers-list (cons '($nary . nary1) *opers-list)) 1048 1049(defun nary1 (e z) 1050 (oper-apply (nary2 e z) z)) 1051 1052(defun nary2 (e z) 1053 (do 1054 ((l (cdr e) (cdr l)) (ans) (some-change)) 1055 1056 ((null l) 1057 (if some-change 1058 (nary2 (cons (car e) (nreverse ans)) z) 1059 (simpargs e z))) 1060 1061 (setq 1062 ans (if (and (not (atom (car l))) (eq (caaar l) (caar e))) 1063 (progn 1064 (setq some-change t) 1065 (nconc (reverse (cdar l)) ans)) 1066 (cons (car l) ans))))) 1067 1068(setq opers (cons '$lassociative opers) 1069 *opers-list (cons '($lassociative . lassociative) *opers-list)) 1070 1071(defun lassociative (e z) 1072 (let* 1073 ((ans0 (oper-apply (cons (car e) (total-nary e)) z)) 1074 (ans (if (consp ans0) (cdr ans0)))) 1075 (cond ((or (null (cddr ans)) (not (eq (caar ans0) (caar e)))) ans0) 1076 ((do ((newans (list (car e) (car ans) (cadr ans)) 1077 (list (car e) newans (car ans))) 1078 (ans (cddr ans) (cdr ans))) 1079 ((null ans) newans)))))) 1080 1081(setq opers (cons '$rassociative opers) 1082 *opers-list (cons '($rassociative . rassociative) *opers-list)) 1083 1084(defun rassociative (e z) 1085 (let* 1086 ((ans0 (oper-apply (cons (car e) (total-nary e)) z)) 1087 (ans (if (consp ans0) (cdr ans0)))) 1088 (cond ((or (null (cddr ans)) (not (eq (caar ans0) (caar e)))) ans0) 1089 (t (setq ans (nreverse ans)) 1090 (do ((newans (list (car e) (cadr ans) (car ans)) 1091 (list (car e) (car ans) newans)) 1092 (ans (cddr ans) (cdr ans))) 1093 ((null ans) newans)))))) 1094 1095(defun total-nary (e) 1096 (do ((l (cdr e) (cdr l)) (ans)) 1097 ((null l) (nreverse ans)) 1098 (setq ans (if (and (not (atom (car l))) (eq (caaar l) (caar e))) 1099 (nconc (reverse (total-nary (car l))) ans) 1100 (cons (car l) ans))))) 1101 1102(defparameter $opproperties (cons '(mlist simp) (reverse opers))) 1103