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 combin) 14 15(declare-top (special *mfactl *factlist donel nn* dn* *ans* *var* 16 $zerobern *n $cflength *a* *a $prevfib $next_lucas 17 *infsumsimp *times *plus sum usum makef 18 varlist genvar $sumsplitfact $ratfac $simpsum 19 $prederror $listarith 20 $ratprint $zeta%pi $bftorat)) 21 22(load-macsyma-macros mhayat rzmac ratmac) 23 24;; minfactorial and factcomb stuff 25 26(defmfun $makefact (e) 27 (let ((makef t)) (if (atom e) e (simplify (makefact1 e))))) 28 29(defun makefact1 (e) 30 (cond ((atom e) e) 31 ((eq (caar e) '%binomial) 32 (subst (makefact1 (cadr e)) 'x 33 (subst (makefact1 (caddr e)) 'y 34 '((mtimes) ((mfactorial) x) 35 ((mexpt) ((mfactorial) y) -1) 36 ((mexpt) ((mfactorial) ((mplus) x ((mtimes) -1 y))) 37 -1))))) 38 ((eq (caar e) '%gamma) 39 (list '(mfactorial) (list '(mplus) -1 (makefact1 (cadr e))))) 40 ((eq (caar e) '$beta) 41 (makefact1 (subst (cadr e) 'x 42 (subst (caddr e) 'y 43 '((mtimes) ((%gamma) x) 44 ((%gamma) y) 45 ((mexpt) ((%gamma) ((mplus) x y)) -1)))))) 46 (t (recur-apply #'makefact1 e)))) 47 48(defmfun $makegamma (e) 49 (if (atom e) e (simplify (makegamma1 ($makefact e))))) 50 51(defmfun $minfactorial (e) 52 (let (*mfactl *factlist) 53 (if (specrepp e) (setq e (specdisrep e))) 54 (getfact e) 55 (mapl #'evfac1 *factlist) 56 (setq e (evfact e)))) 57 58(defun evfact (e) 59 (cond ((atom e) e) 60 ((eq (caar e) 'mfactorial) 61 ;; Replace factorial with simplified expression from *factlist. 62 (simplifya (cdr (assoc (cadr e) *factlist :test #'equal)) nil)) 63 ((member (caar e) '(%sum %derivative %integrate %product) :test #'eq) 64 (cons (list (caar e)) (cons (evfact (cadr e)) (cddr e)))) 65 (t (recur-apply #'evfact e)))) 66 67(defun adfactl (e l) 68 (let (n) 69 (cond ((null l) (push (list e) *mfactl)) 70 ((numberp (setq n ($ratsimp `((mplus) ,e ((mtimes) -1 ,(caar l)))))) 71 (cond ((plusp n) 72 (rplacd (car l) (cons e (cdar l)))) 73 ((rplaca l (cons e (car l)))))) 74 ((adfactl e (cdr l)))))) 75 76(defun getfact (e) 77 (cond ((atom e) nil) 78 ((eq (caar e) 'mfactorial) 79 (and (null (member (cadr e) *factlist :test #'equal)) 80 (prog2 81 (push (cadr e) *factlist) 82 (adfactl (cadr e) *mfactl)))) 83 ((member (caar e) '(%sum %derivative %integrate %product) :test #'eq) 84 (getfact (cadr e))) 85 ((mapc #'getfact (cdr e))))) 86 87(defun evfac1 (e) 88 (do ((al *mfactl (cdr al))) 89 ((member (car e) (car al) :test #'equal) 90 (rplaca e 91 (cons (car e) 92 (list '(mtimes) 93 (gfact (car e) 94 ($ratsimp (list '(mplus) (car e) 95 (list '(mtimes) -1 (caar al)))) 1) 96 (list '(mfactorial) (caar al)))))))) 97 98(defmfun $factcomb (e) 99 (let ((varlist varlist ) genvar $ratfac (ratrep (and (not (atom e)) (eq (caar e) 'mrat)))) 100 (and ratrep (setq e (ratdisrep e))) 101 (setq e (factcomb e) 102 e (cond ((atom e) e) 103 (t (simplify (cons (list (caar e)) 104 (mapcar #'factcomb1 (cdr e))))))) 105 (or $sumsplitfact (setq e ($minfactorial e))) 106 (if ratrep (ratf e) e))) 107 108(defun factcomb1 (e) 109 (cond ((free e 'mfactorial) e) 110 ((member (caar e) '(mplus mtimes mexpt) :test #'eq) 111 (cons (list (caar e)) (mapcar #'factcomb1 (cdr e)))) 112 (t (setq e (factcomb e)) 113 (if (atom e) 114 e 115 (cons (list (caar e)) (mapcar #'factcomb1 (cdr e))))))) 116 117(defun factcomb (e) 118 (cond ((atom e) e) 119 ((free e 'mfactorial) e) 120 ((member (caar e) '(mplus mtimes) :test #'eq) 121 (factpluscomb (factcombplus e))) 122 ((eq (caar e) 'mexpt) 123 (simpexpt (list '(mexpt) (factcomb (cadr e)) 124 (factcomb (caddr e))) 125 1 nil)) 126 ((eq (caar e) 'mrat) 127 (factrat e)) 128 (t (cons (car e) (mapcar #'factcomb (cdr e)))))) 129 130(defun factrat (e) 131 (let (nn* dn*) 132 (setq e (factqsnt ($ratdisrep (cons (car e) (cons (cadr e) 1))) 133 ($ratdisrep (cons (car e) (cons (cddr e) 1))))) 134 (numden e) 135 (div* (factpluscomb nn*) (factpluscomb dn*)))) 136 137(defun factqsnt (num den) 138 (if (equal num 0) 0 139 (let (nn* dn* (e (factpluscomb (div* den num)))) 140 (numden e) 141 (factpluscomb (div* dn* nn*))))) 142 143(defun factcombplus (e) 144 (let (nn* dn*) 145 (do ((l1 (nplus e) (cdr l1)) 146 (l2)) 147 ((null l1) 148 (simplus (cons '(mplus) 149 (mapcar #'(lambda (q) (factqsnt (car q) (cdr q))) l2)) 150 1 nil)) 151 (numden (car l1)) 152 (do ((l3 l2 (cdr l3)) 153 (l4)) 154 ((null l3) (setq l2 (nconc l2 (list (cons nn* dn*))))) 155 (setq l4 (car l3)) 156 (cond ((not (free ($gcd dn* (cdr l4)) 'mfactorial)) 157 (numden (list '(mplus) (div* nn* dn*) 158 (div* (car l4) (cdr l4)))) 159 (setq l2 (delete l4 l2 :count 1 :test #'eq)))))))) 160 161(defun factpluscomb (e) 162 (prog (donel fact indl tt) 163 tag (setq e (factexpand e) 164 fact (getfactorial e)) 165 (or fact (return e)) 166 (setq indl (mapcar #'(lambda (q) (factplusdep q fact)) 167 (nplus e)) 168 tt (factpowerselect indl (nplus e) fact) 169 e (cond ((cdr tt) 170 (cons '(mplus) (mapcar #'(lambda (q) (factplus2 q fact)) 171 tt))) 172 (t (factplus2 (car tt) fact)))) 173 (go tag))) 174 175(defun nplus (e) 176 (if (eq (caar e) 'mplus) 177 (cdr e) 178 (list e))) 179 180(defun factexpand (e) 181 (cond ((atom e) e) 182 ((eq (caar e) 'mplus) 183 (simplus (cons '(mplus) (mapcar #'factexpand (cdr e))) 184 1 nil)) 185 ((free e 'mfactorial) e) 186 (t ($expand e)))) 187 188(defun getfactorial (e) 189 (cond ((atom e) nil) 190 ((member (caar e) '(mplus mtimes) :test #'eq) 191 (do ((e (cdr e) (cdr e)) 192 (a)) 193 ((null e) nil) 194 (setq a (getfactorial (car e))) 195 (and a (return a)))) 196 ((eq (caar e) 'mexpt) 197 (getfactorial (cadr e))) 198 ((eq (caar e) 'mfactorial) 199 (and (null (memalike (cadr e) donel)) 200 (list '(mfactorial) 201 (car (setq donel (cons (cadr e) donel)))))))) 202 203(defun factplusdep (e fact) 204 (cond ((alike1 e fact) 1) 205 ((atom e) nil) 206 ((eq (caar e) 'mtimes) 207 (do ((l (cdr e) (cdr l)) 208 (e) (out)) 209 ((null l) nil) 210 (setq e (car l)) 211 (and (setq out (factplusdep e fact)) 212 (return out)))) 213 ((eq (caar e) 'mexpt) 214 (let ((fto (factplusdep (cadr e) fact))) 215 (and fto (simptimes (list '(mtimes) fto 216 (caddr e)) 1 t)))) 217 ((eq (caar e) 'mplus) 218 (same (mapcar #'(lambda (q) (factplusdep q fact)) 219 (cdr e)))))) 220 221(defun same (l) 222 (do ((ca (car l)) 223 (cd (cdr l) (cdr cd)) 224 (cad)) 225 ((null cd) ca) 226 (setq cad (car cd)) 227 (or (alike1 ca cad) 228 (return nil)))) 229 230(defun factpowerselect (indl e fact) 231 (let (l fl) 232 (do ((i indl (cdr i)) 233 (j e (cdr j)) 234 (expt) (exp)) 235 ((null i) l) 236 (setq expt (car i) 237 exp (cond (expt 238 (setq exp ($divide (car j) `((mexpt) ,fact ,expt))) 239 ;; (car j) need not involve fact^expt since 240 ;; fact^expt may be the gcd of the num and denom 241 ;; of (car j) and $divide will cancel this out. 242 (if (not (equal (cadr exp) 0)) 243 (cadr exp) 244 (progn 245 (setq expt '()) 246 (caddr exp)))) 247 (t (car j)))) 248 (cond ((null l) (setq l (list (list expt exp)))) 249 ((setq fl (assolike expt l)) 250 (nconc fl (list exp))) 251 (t (nconc l (list (list expt exp)))))))) 252 253(defun factplus2 (l fact) 254 (let ((expt (car l))) 255 (cond (expt (factplus0 (cond ((cddr l) (rplaca l '(mplus))) 256 (t (cadr l))) 257 expt (cadr fact))) 258 (t (rplaca l '(mplus)))))) 259 260(defun factplus0 (r e fact) 261 (do ((i -1 (1- i)) 262 (fpn fact (list '(mplus) fact i)) 263 (j -1) (exp) (rfpn) (div)) 264 (nil) 265 (setq rfpn (simpexpt (list '(mexpt) fpn -1) 1 nil)) 266 (setq div (dypheyed r (simpexpt (list '(mexpt) rfpn e) 1 nil))) 267 (cond ((or (null (or $sumsplitfact (equal (cadr div) 0))) 268 (equal (car div) 0)) 269 (return (simplus (cons '(mplus) (mapcar 270 #'(lambda (q) 271 (incf j) 272 (list '(mtimes) q (list '(mexpt) 273 (list '(mfactorial) (list '(mplus) fpn j)) e))) 274 (factplus1 (cons r exp) e fpn))) 275 1 nil))) 276 (t (setq r (car div)) 277 (setq exp (cons (cadr div) exp)))))) 278 279(defun factplus1 (exp e fact) 280 (do ((l exp (cdr l)) 281 (i 2 (1+ i)) 282 (fpn (list '(mplus) fact 1) (list '(mplus) fact i)) 283 (div)) 284 ((null l) exp) 285 (setq div (dypheyed (car l) (list '(mexpt) fpn e))) 286 (and (or $sumsplitfact (equal (cadr div) 0)) 287 (null (equal (car div) 0)) 288 (rplaca l (cadr div)) 289 (rplacd l (cons (cond ((cadr l) 290 (simplus (list '(mplus) (car div) (cadr l)) 291 1 nil)) 292 (t 293 (setq donel 294 (cons (simplus fpn 1 nil) donel)) 295 (car div))) 296 (cddr l)))))) 297 298(defun dypheyed (r f) 299 (let (r1 p1 p2) 300 (newvar r) 301 (setq r1 (ratf f) 302 p1 (pdegreevector (cadr r1)) 303 p2 (pdegreevector (cddr r1))) 304 (do ((i p1 (cdr i)) 305 (j p2 (cdr j)) 306 (k (caddar r1) (cdr k))) 307 ((null k) (kansel r (cadr r1) (cddr r1))) 308 (cond ((> (car i) (car j)) 309 (return (cdr ($divide r f (car k))))))))) 310 311(defun kansel (r n d) 312 (let (r1 p1 p2) 313 (setq r1 (ratf r) 314 p1 (testdivide (cadr r1) n) 315 p2 (testdivide (cddr r1) d)) 316 (if (and p1 p2) 317 (cons (rdis (cons p1 p2)) '(0)) 318 (cons '0 (list r))))) 319 320;; euler and bernoulli stuff 321 322(defvar *bn* (make-array 17 :adjustable t :element-type 'integer 323 :initial-contents '(0 -1 1 -1 5. -691. 7. -3617. 43867. -174611. 854513. 324 -236364091. 8553103. -23749461029. 8615841276005. 325 -7709321041217. 2577687858367.))) 326 327(defvar *bd* (make-array 17 :adjustable t :element-type 'integer 328 :initial-contents '(0 30. 42. 30. 66. 2730. 6. 510. 798. 330. 138. 2730. 329 6. 870. 14322. 510. 6.))) 330 331(defvar *eu* (make-array 11 :adjustable t :element-type 'integer 332 :initial-contents '(-1 5. -61. 1385. -50521. 2702765. -199360981. 19391512145. 333 -2404879675441. 370371188237525. -69348874393137901.))) 334 335(putprop '*eu* 11 'lim) 336(putprop 'bern 16 'lim) 337 338(defmfun $euler (s) 339 (setq s 340 (let ((%n 0) $float) 341 (cond ((or (not (fixnump s)) (< s 0)) (list '($euler) s)) 342 ((zerop (setq %n s)) 1) 343 ($zerobern 344 (cond ((oddp %n) 0) 345 ((null (> (ash %n -1) (get '*eu* 'lim))) 346 (aref *eu* (1- (ash %n -1)))) 347 ((eq $zerobern '%$/#&) 348 (euler %n)) 349 ((setq *eu* (adjust-array *eu* (1+ (ash %n -1)))) 350 (euler %n)))) 351 ((<= %n (get '*eu* 'lim)) 352 (aref *eu* (1- %n))) 353 ((setq *eu* (adjust-array *eu* (1+ %n))) 354 (euler (* 2 %n)))))) 355 (simplify s)) 356 357(defun nxtbincoef (m nom) 358 (truncate (* nom (- *a* m)) m)) 359 360(defun euler (%a*) 361 (prog (nom %k e fl $zerobern *a*) 362 (setq nom 1 %k %a* fl nil e 0 $zerobern '%$/#& *a* (1+ %a*)) 363 a (cond ((zerop %k) 364 (setq e (- e)) 365 (setf (aref *eu* (1- (ash %a* -1))) e) 366 (putprop '*eu* (ash %a* -1) 'lim) 367 (return e))) 368 (setq nom (nxtbincoef (1+ (- %a* %k)) nom) %k (1- %k)) 369 (cond ((setq fl (null fl)) 370 (go a))) 371 (incf e (* nom ($euler %k))) 372 (go a))) 373 374(defun simpeuler (x vestigial z) 375 (declare (ignore vestigial)) 376 (oneargcheck x) 377 (let ((u (simpcheck (cadr x) z))) 378 (if (and (fixnump u) (>= u 0)) 379 ($euler u) 380 (eqtest (list '($euler) u) x)))) 381 382(defmfun $bern (s) 383 (setq s 384 (let ((%n 0) $float) 385 (cond ((or (not (fixnump s)) (< s 0)) (list '($bern) s)) 386 ((= (setq %n s) 0) 1) 387 ((= %n 1) '((rat) -1 2)) 388 ((= %n 2) '((rat) 1 6)) 389 ($zerobern 390 (cond ((oddp %n) 0) 391 ((null (> (setq %n (1- (ash %n -1))) (get 'bern 'lim))) 392 (list '(rat) (aref *bn* %n) (aref *bd* %n))) 393 ((eq $zerobern '$/#&) (bern (* 2 (1+ %n)))) 394 (t 395 (setq *bn* (adjust-array *bn* (setq %n (1+ %n)))) 396 (setq *bd* (adjust-array *bd* %n)) 397 (bern (* 2 %n))))) 398 ((null (> %n (get 'bern 'lim))) 399 (list '(rat) (aref *bn* (- %n 2)) (aref *bd* (- %n 2)))) 400 (t 401 (setq *bn* (adjust-array *bn* (1+ %n))) 402 (setq *bd* (adjust-array *bd* (1+ %n))) 403 (bern (* 2 (1- %n))))))) 404 (simplify s)) 405 406(defun bern (%a*) 407 (prog (nom %k bb a b $zerobern l *a*) 408 (setq %k 0 409 l (1- %a*) 410 %a* (1+ %a*) 411 nom 1 412 $zerobern '$/#& 413 a 1 414 b 1 415 *a* (1+ %a*)) 416 a (cond ((= %k l) 417 (setq bb (*red a (* -1 b %a*))) 418 (putprop 'bern (setq %a* (1- (ash %a* -1))) 'lim) 419 (setf (aref *bn* %a*) (cadr bb)) 420 (setf (aref *bd* %a*) (caddr bb)) 421 (return bb))) 422 (incf %k) 423 (setq a (+ (* b (setq nom (nxtbincoef %k nom)) 424 (num1 (setq bb ($bern %k)))) 425 (* a (denom1 bb)))) 426 (setq b (* b (denom1 bb))) 427 (setq a (*red a b) b (denom1 a) a (num1 a)) 428 (go a))) 429 430(defun simpbern (x vestigial z) 431 (declare (ignore vestigial)) 432 (oneargcheck x) 433 (let ((u (simpcheck (cadr x) z))) 434 (if (and (fixnump u) (not (< u 0))) 435 ($bern u) 436 (eqtest (list '($bern) u) x)))) 437 438;;; ---------------------------------------------------------------------------- 439;;; Bernoulli polynomials 440;;; 441;;; The following explicit formula is directly implemented: 442;;; 443;;; n 444;;; ==== 445;;; \ n - k 446;;; B (x) = > b binomial(n, k) x 447;;; n / k 448;;; ==== 449;;; k = 0 450;;; 451;;; The coeffizients b[k] are the Bernoulli numbers. The algorithm does not 452;;; skip over Beroulli numbers, which are zero. We have to ensure that 453;;; $zerobern is bound to true. 454;;; ---------------------------------------------------------------------------- 455 456(defmfun $bernpoly (x s) 457 (let ((%n 0) ($zerobern t)) 458 (cond ((not (fixnump s)) (list '($bernpoly) x s)) 459 ((> (setq %n s) -1) 460 (do ((sum (cons (if (and (= %n 0) (zerop1 x)) 461 (add 1 x) 462 (power x %n)) 463 nil) 464 (cons (mul (binocomp %n %k) 465 ($bern %k) 466 (if (and (= %n %k) (zerop1 x)) 467 (add 1 x) 468 (power x (- %n %k)))) 469 sum)) 470 (%k 1 (1+ %k))) 471 ((> %k %n) (addn sum t)))) 472 (t (list '($bernpoly) x %n))))) 473 474;;; ---------------------------------------------------------------------------- 475;;; Euler polynomials 476;;; 477;;; The following explicit formula is directly implemented: 478;;; 479;;; n 1 n - k 480;;; ==== E binomial(n, k) (x - -) 481;;; \ k 2 482;;; E (x) = > ------------------------------ 483;;; n / k 484;;; ==== 2 485;;; k = 0 486;;; 487;;; The coeffizients E[k] are the Euler numbers. 488;;; ---------------------------------------------------------------------------- 489 490(defmfun $eulerpoly (x s) 491 (let ((n 0) ($zerobern t) (y 0)) 492 (cond ((not (fixnump s)) (list '($eulerpoly) x s)) 493 ((> (setq n s) -1) 494 (do ((sum (cons (if (and (zerop1 (setq y (sub x (div 1 2)))) 495 (= n 0)) 496 (add 1 y) 497 (power y n)) 498 nil) 499 (cons (mul (binocomp n k) 500 ($euler k) 501 (power 2 (mul -1 k)) 502 (if (and (zerop1 (setq y (sub x (div 1 2)))) 503 (= n k)) 504 (add 1 y) 505 (power y (- n k)))) 506 sum)) 507 (k 1 (1+ k))) 508 ((> k n) ($expand (addn sum t))))) 509 (t (list '($eulerpoly) x n))))) 510 511;; zeta and fibonacci stuff 512 513;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 514;;; 515;;; Implementation of the Riemann Zeta function as a simplifying function 516;;; 517;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 518 519(defmfun $zeta (z) 520 (simplify (list '(%zeta) z))) 521 522;;; Set properties to give full support to the parser and display 523 524(defprop $zeta %zeta alias) 525(defprop $zeta %zeta verb) 526 527(defprop %zeta $zeta reversealias) 528(defprop %zeta $zeta noun) 529 530;;; The Riemann Zeta function is a simplifying function 531 532(defprop %zeta simp-zeta operators) 533 534;;; The Riemann Zeta function has mirror symmetry 535 536(defprop %zeta t commutes-with-conjugate) 537 538;;; The Riemann Zeta function distributes over lists, matrices, and equations 539 540(defprop %zeta (mlist $matrix mequal) distribute_over) 541 542;;; We support a simplim%function. The function is looked up in simplimit and 543;;; handles specific values of the function. 544 545(defprop %zeta simplim%zeta simplim%function) 546 547(defun simplim%zeta (expr var val) 548 ;; Look for the limit of the argument 549 (let* ((arg (limit (cadr expr) var val 'think)) 550 (dir (limit (add (cadr expr) (neg arg)) var val 'think))) 551 (cond 552 ;; Handle an argument 1 at this place 553 ((onep1 arg) 554 (cond ((eq dir '$zeroa) 555 '$inf) 556 ((eq dir '$zerob) 557 '$minf) 558 (t '$infinity))) 559 (t 560 ;; All other cases are handled by the simplifier of the function. 561 (simplify (list '(%zeta) arg)))))) 562 563;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 564 565(defun simp-zeta (expr z simpflag) 566 (oneargcheck expr) 567 (setq z (simpcheck (cadr expr) simpflag)) 568 (cond 569 570 ;; Check for special values 571 ((eq z '$inf) 1) 572 ((alike1 z '((mtimes) -1 $minf)) 1) 573 ((zerop1 z) 574 (cond (($bfloatp z) ($bfloat '((rat) -1 2))) 575 ((floatp z) -0.5) 576 (t '((rat simp) -1 2)))) 577 ((onep1 z) 578 (simp-domain-error (intl:gettext "zeta: zeta(~:M) is undefined.") z)) 579 580 ;; Check for numerical evaluation 581 ((or (bigfloat-numerical-eval-p z) 582 (complex-bigfloat-numerical-eval-p z) 583 (float-numerical-eval-p z) 584 (complex-float-numerical-eval-p z)) 585 (to (float-zeta z))) 586 ;; Check for transformations and argument simplifications 587 ((integerp z) 588 (cond 589 ((oddp z) 590 (cond ((> z 1) 591 (eqtest (list '(%zeta) z) expr)) 592 ((setq z (sub 1 z)) 593 (mul -1 (div ($bern z) z))))) 594 ((minusp z) 0) 595 ((not $zeta%pi) (eqtest (list '(%zeta) z) expr)) 596 (t (let ($numer $float) 597 (mul (power '$%pi z) 598 (mul (div (power 2 (1- z)) 599 (take '(mfactorial) z)) 600 (take '(mabs) ($bern z)))))))) 601 (t 602 (eqtest (list '(%zeta) z) expr)))) 603 604;; See http://numbers.computation.free.fr/Constants/constants.html 605;; and, in particular, 606;; http://numbers.computation.free.fr/Constants/Miscellaneous/zetaevaluations.pdf. 607;; We use the algorithm from Proposition 2: 608;; 609;; zeta(s) = 1/(1-2^(1-s)) * 610;; (sum((-1)^(k-1)/k^s,k,1,n) + 611;; 1/2^n*sum((-1)^(k-1)*e(k-n)/k^s,k,n+1,2*n)) 612;; + g(n,s) 613;; 614;; where e(k) = sum(binomial(n,j), j, k, n). Writing s = sigma + %i*t, when 615;; sigma is positive you get an error estimate of 616;; 617;; |g(n,s)| <= 1/8^n * h(s) 618;; 619;; where 620;; 621;; h(s) = ((1 + abs (t / sigma)) exp (abs (t) * %pi / 2)) / abs (1 - 2^(1 - s)) 622;; 623;; We need to figure out how many terms are required to make |g(n,s)| 624;; sufficiently small. The answer is 625;; 626;; n = (log h(s) - log (eps)) / log (8) 627;; 628;; and 629;; 630;; log (h (s)) = (%pi/2 * abs (t)) + log (1 + t/sigma) - log (abs (1 - 2^(1 - s))) 631;; 632;; Notice that this bound is a bit rubbish when sigma is near zero. In that 633;; case, use the expansion zeta(s) = -1/2-1/2*log(2*pi)*s. 634(defun float-zeta (s) 635 ;; If s is a rational (real or complex), convert to a float. This 636 ;; is needed so we can compute a sensible epsilon value. (What is 637 ;; the epsilon value for an exact rational?) 638 (setf s (bigfloat:to s)) 639 (typecase s 640 (rational 641 (setf s (float s))) 642 ((complex rational) 643 (setf s (coerce s '(complex flonum))))) 644 645 (let ((sigma (bigfloat:realpart s))) 646 (cond 647 ;; abs(s)^2 < epsilon, use the expansion zeta(s) = -1/2-1/2*log(2*%pi)*s 648 ((bigfloat:< (bigfloat:abs (bigfloat:* s s)) (bigfloat:epsilon s)) 649 (bigfloat:+ -1/2 650 (bigfloat:* -1/2 651 (bigfloat:log (bigfloat:* 2 (bigfloat:%pi s))) 652 s))) 653 654 ;; Reflection formula 655 ((bigfloat:minusp sigma) 656 (let ((n (bigfloat:floor sigma))) 657 ;; If sigma is a negative even integer, zeta(sigma) is zero, 658 ;; from the reflection formula because sin(%pi*n/2) is 0. 659 (when (and (bigfloat:= n sigma) (evenp n)) 660 (return-from float-zeta (bigfloat:float 0.0 sigma)))) 661 (bigfloat:* (bigfloat:expt 2 s) 662 (bigfloat:expt (bigfloat:%pi s) 663 (bigfloat:- s 1)) 664 (bigfloat:sin (bigfloat:* (bigfloat:/ (bigfloat:%pi s) 665 2) 666 s)) 667 (bigfloat:to ($gamma (to (bigfloat:- 1 s)))) 668 (float-zeta (bigfloat:- 1 s)))) 669 670 ;; The general formula from above. Call the imaginary part "tau" rather 671 ;; than the "t" above, because that isn't a CL keyword... 672 (t 673 (let* ((tau (bigfloat:imagpart s)) 674 (logh 675 (bigfloat:- 676 (if (bigfloat:zerop tau) 0 677 (bigfloat:+ 678 (bigfloat:* 1.6 (bigfloat:abs tau)) 679 (bigfloat:log (bigfloat:1+ 680 (bigfloat:abs 681 (bigfloat:/ tau sigma)))))) 682 (bigfloat:log 683 (bigfloat:abs 684 (bigfloat:- 1 (bigfloat:expt 2 (bigfloat:- 1 s))))))) 685 686 (logeps (bigfloat:log (bigfloat:epsilon s))) 687 688 (n (max (bigfloat:ceiling 689 (bigfloat:/ (bigfloat:- logh logeps) (bigfloat:log 8))) 690 1)) 691 692 (sum1 0) 693 (sum2 0)) 694 (flet ((binsum (k n) 695 ;; sum(binomial(n,j), j, k, n) = sum(binomial(n,j), j, n, k) 696 (let ((sum 0) 697 (term 1)) 698 (loop for j from n downto k 699 do 700 (progn 701 (incf sum term) 702 (setf term (/ (* term j) (+ n 1 (- j)))))) 703 sum))) 704 ;; (format t "n = ~D terms~%" n) 705 ;; sum1 = sum((-1)^(k-1)/k^s,k,1,n) 706 ;; sum2 = sum((-1)^(k-1)/e(n,k-n)/k^s, k, n+1, 2*n) 707 ;; = (-1)^n*sum((-1)^(m-1)*e(n,m)/(n+k)^s, k, 1, n) 708 (loop for k from 1 to n 709 for d = (bigfloat:expt k s) 710 do (progn 711 (bigfloat:incf sum1 (bigfloat:/ (cl:expt -1 (- k 1)) d)) 712 (bigfloat:incf sum2 (bigfloat:/ (* (cl:expt -1 (- k 1)) 713 (binsum k n)) 714 (bigfloat:expt (+ k n) s)))) 715 finally (return (values sum1 sum2))) 716 (when (oddp n) 717 (setq sum2 (bigfloat:- sum2))) 718 (bigfloat:/ (bigfloat:+ sum1 719 (bigfloat:/ sum2 (bigfloat:expt 2 n))) 720 (bigfloat:- 1 (bigfloat:expt 2 (bigfloat:- 1 s)))))))))) 721 722;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 723 724(defmfun $fib (n) 725 (cond ((fixnump n) (ffib n)) 726 (t (setq $prevfib `(($fib) ,(add2* n -1))) 727 `(($fib) ,n)))) 728 729(defun ffib (%n) 730 (declare (fixnum %n)) 731 (cond ((= %n -1) 732 (setq $prevfib -1) 733 1) 734 ((zerop %n) 735 (setq $prevfib 1) 736 0) 737 (t 738 (let* ((f2 (ffib (ash (logandc2 %n 1) -1))) ; f2 = fib(n/2) or fib((n-1)/2) 739 (x (+ f2 $prevfib)) 740 (y (* $prevfib $prevfib)) 741 (z (* f2 f2))) 742 (setq f2 (- (* x x) y) 743 $prevfib (+ y z)) 744 (when (oddp %n) 745 (psetq $prevfib f2 746 f2 (+ f2 $prevfib))) 747 f2)))) 748 749(defmfun $lucas (n) 750 (cond 751 ((fixnump n) (lucas n)) 752 (t (setq $next_lucas `(($lucas) ,(add2* n 1))) 753 `(($lucas) ,n) ))) 754 755(defun lucas (n) 756 (declare (fixnum n)) 757 (let ((w 2) (x 2) (y 1) u v (sign (signum n))) (declare (fixnum sign)) 758 (setq n (abs n)) 759 (do ((i (1- (integer-length n)) (1- i))) 760 ((< i 0)) 761 (declare (fixnum i)) 762 (setq u (* x x) v (* y y)) 763 (if (logbitp i n) 764 (setq y (+ v w) x (+ y (- u) w) w -2) 765 (setq x (- u w) y (+ v w (- x)) w 2) )) 766 (cond 767 ((or (= 1 sign) (not (logbitp 0 n))) 768 (setq $next_lucas y) 769 x ) 770 (t 771 (setq $next_lucas (neg y)) 772 (neg x) )))) 773 774;; continued fraction stuff 775 776(defmfun $cfdisrep (a) 777 (cond ((not ($listp a)) 778 (merror (intl:gettext "cfdisrep: argument must be a list; found ~M") a)) 779 ((null (cddr a)) (cadr a)) 780 ((equal (cadr a) 0) 781 (list '(mexpt) (cfdisrep1 (cddr a)) -1)) 782 ((cfdisrep1 (cdr a))))) 783 784(defun cfdisrep1 (a) 785 (cond ((cdr a) 786 (list '(mplus simp cf) (car a) 787 (prog2 (setq a (cfdisrep1 (cdr a))) 788 (cond ((integerp a) (list '(rat simp) 1 a)) 789 (t (list '(mexpt simp) a -1)))))) 790 ((car a)))) 791 792(defun cfmak (a) 793 (setq a (meval a)) 794 (cond ((integerp a) (list a)) 795 ((eq (caar a) 'mlist) (cdr a)) 796 ((eq (caar a) 'rat) (ratcf (cadr a) (caddr a))) 797 ((merror (intl:gettext "cf: continued fractions must be lists or integers; found ~M") a)))) 798 799(defun makcf (a) 800 (cond ((null (cdr a)) (car a)) 801 ((cons '(mlist simp cf) a)))) 802 803;;; Translation properties for $CF defined in MAXSRC;TRANS5 > 804 805(defmspec $cf (a) 806 (cfratsimp (let ($listarith) 807 (cfeval (meval (fexprcheck a)))))) 808 809;; Definition of cfratsimp as given in SF bug report # 620928. 810(defun cfratsimp (a) 811 (cond ((atom a) a) 812 ((member 'cf (car a) :test #'eq) a) 813 (t (cons '(mlist cf simp) 814 (apply 'find-cf (cf-back-recurrence (cdr a))))))) 815 816; Code to expand nth degree roots of integers into continued fraction 817; approximations. E.g. cf(2^(1/3)) 818; Courtesy of Andrei Zorine (feniy@mail.nnov.ru) 2005/05/07 819 820(defun cfnroot(b) 821 (let ((ans (list '(mlist xf))) ent ($algebraic $true)) 822 (dotimes (i $cflength (nreverse ans)) 823 (setq ent (meval `(($floor) ,b)) 824 ans (cons ent ans) 825 b ($ratsimp (m// (m- b ent))))))) 826 827(defun cfeval (a) 828 (let (temp $ratprint) 829 (cond ((integerp a) (list '(mlist cf) a)) 830 ((floatp a) 831 (let ((a (maxima-rationalize a))) 832 (cons '(mlist cf) (ratcf (car a) (cdr a))))) 833 (($bfloatp a) 834 (let (($bftorat t)) 835 (setq a (bigfloat2rat a)) 836 (cons '(mlist cf) (ratcf (car a) (cdr a))))) 837 ((atom a) 838 (merror (intl:gettext "cf: ~:M is not a continued fraction.") a)) 839 ((eq (caar a) 'rat) 840 (cons '(mlist cf) (ratcf (cadr a) (caddr a)))) 841 ((eq (caar a) 'mlist) 842 (cfratsimp a)) 843 ;;the following doesn't work for non standard form 844 ;; (cfplus a '((mlist) 0))) 845 ((and (mtimesp a) (cddr a) (null (cdddr a)) 846 (fixnump (cadr a)) 847 (mexptp (caddr a)) 848 (fixnump (cadr (caddr a))) 849 (alike1 (caddr (caddr a)) '((rat) 1 2))) 850 (cfsqrt (cfeval (* (expt (cadr a) 2) (cadr (caddr a)))))) 851 ((eq (caar a) 'mexpt) 852 (cond ((alike1 (caddr a) '((rat) 1 2)) 853 (cfsqrt (cfeval (cadr a)))) 854 ((integerp (m* 2 (caddr a))) ; a^(n/2) was sqrt(a^n) 855 (cfsqrt (cfeval (cfexpt (cadr a) (m* 2 (caddr a)))))) 856 ((integerp (cadr a)) (cfnroot a)) ; <=== new case x 857 ((cfexpt (cfeval (cadr a)) (caddr a))))) 858 ((setq temp (assoc (caar a) '((mplus . cfplus) (mtimes . cftimes) (mquotient . cfquot) 859 (mdifference . cfdiff) (mminus . cfminus)) :test #'eq)) 860 (cf (cfeval (cadr a)) (cddr a) (cdr temp))) 861 ((eq (caar a) 'mrat) 862 (cfeval ($ratdisrep a))) 863 (t (merror (intl:gettext "cf: ~:M is not a continued fraction.") a))))) 864 865(defun cf (a l fun) 866 (cond ((null l) a) 867 ((cf (funcall fun a (meval (list '($cf) (car l)))) (cdr l) fun)))) 868 869(defun cfplus (a b) 870 (setq a (cfmak a) b (cfmak b)) 871 (makcf (cffun '(0 1 1 0) '(0 0 0 1) a b))) 872 873(defun cftimes (a b) 874 (setq a (cfmak a) b (cfmak b)) 875 (makcf (cffun '(1 0 0 0) '(0 0 0 1) a b))) 876 877(defun cfdiff (a b) 878 (setq a (cfmak a) b (cfmak b)) 879 (makcf (cffun '(0 1 -1 0) '(0 0 0 1) a b))) 880 881(defun cfmin (a) 882 (setq a (cfmak a)) 883 (makcf (cffun '(0 0 -1 0) '(0 0 0 1) a '(0)))) 884 885(defun cfquot (a b) 886 (setq a (cfmak a) b (cfmak b)) 887 (makcf (cffun '(0 1 0 0) '(0 0 1 0) a b))) 888 889(defun cfexpt (b e) 890 (setq b (cfmak b)) 891 (cond ((null (integerp e)) 892 (merror (intl:gettext "cf: can't raise continued fraction to non-integral power ~M") e)) 893 ((let ((n (abs e))) 894 (do ((n (ash n -1) (ash n -1)) 895 (s (cond ((oddp n) b) 896 (t '(1))))) 897 ((zerop n) 898 (makcf 899 (cond ((signp g e) 900 s) 901 ((cffun '(0 0 0 1) '(0 1 0 0) b '(1)))))) 902 (setq b (cffun '(1 0 0 0) '(0 0 0 1) b b)) 903 (and (oddp n) 904 (setq s (cffun '(1 0 0 0) '(0 0 0 1) s b)))))))) 905 906 907(defun conf1 (f g a b &aux (den (conf2 g a b))) 908 (cond ((zerop den) 909 (* (signum (conf2 f a b )) ; (/ most-positive-fixnum (^ 2 4)) 910 #.(expt 2 31))) 911 (t (truncate (conf2 f a b) den)))) 912 913(defun conf2 (n a b) ;2*(abn_0+an_1+bn_2+n_3) 914 (* 2 (+ (* (car n) a b) 915 (* (cadr n) a) 916 (* (caddr n) b) 917 (cadddr n)))) 918 919;;(cffun '(0 1 1 0) '(0 0 0 1) '(1 2) '(1 1 1 2)) gets error 920;;should give (3 10) 921 922(defun cf-convergents-p-q (cf &optional (n (length cf)) &aux pp qq) 923 "returns two lists such that pp_i/qq_i is the quotient of the first i terms 924 of cf" 925 (case (length cf) 926 (0 1) 927 (1 cf(list 1)) 928 (t 929 (setq pp (list (1+ (* (first cf) (second cf))) (car cf))) 930 (setq qq (list (second cf) 1)) 931 (show pp qq) 932 (setq cf (cddr cf)) 933 (loop for i from 2 to n 934 while cf 935 do 936 (push (+ (* (car cf) (car pp)) 937 (second pp)) pp) 938 (push (+ (* (car cf) (car qq)) 939 (second qq)) qq) 940 (setq cf (cdr cf)) 941 finally (return (list (reverse pp) (reverse qq))))))) 942 943 944(defun find-cf1 (p q so-far) 945 (multiple-value-bind (quot rem) (truncate p q) 946 (cond ((< rem 0) (incf rem q) (incf quot -1)) 947 ((zerop rem) (return-from find-cf1 (cons quot so-far)))) 948 (setq so-far (cons quot so-far)) 949 (find-cf1 q rem so-far))) 950 951(defun find-cf (p q) 952 "returns the continued fraction for p and q integers, q not zero" 953 (cond ((zerop q) (maxima-error "find-cf: quotient by zero")) 954 ((< q 0) (setq p (- p)) (setq q (- q)))) 955 (nreverse (find-cf1 p q ()))) 956 957(defun cf-back-recurrence (cf &aux tem (num-gg 0)(den-gg 1)) 958 "converts CF (a continued fraction list) to a list of numerator 959 denominator using recurrence from end 960 and not calculating intermediate quotients. 961 The numerator and denom are relatively 962 prime" 963 (loop for v in (reverse cf) 964 do (setq tem (* den-gg v)) 965 (setq tem (+ tem num-gg)) 966 (setq num-gg den-gg) 967 (setq den-gg tem) 968 finally 969 (return 970 (cond ((and (<= den-gg 0) (< num-gg 0)) 971 (list (- den-gg) (- num-gg))) 972 (t(list den-gg num-gg)))))) 973 974(declare-top (unspecial w)) 975 976;;(cffun '(0 1 1 0) '(0 0 0 1) '(1 2) '(1 1 1 2)) gets error 977;;should give (3 10) 978 979(defun cffun (f g a b) 980 (prog (c v w) 981 (declare (special v)) 982 a (and (zerop (cadddr g)) 983 (zerop (caddr g)) 984 (zerop (cadr g)) 985 (zerop (car g)) 986 (return (reverse c))) 987 (and (equal (setq w (conf1 f g (car a) (1+ (car b)))) 988 (setq v (conf1 f g (car a) (car b)))) 989 (equal (conf1 f g (1+ (car a)) (car b)) v) 990 (equal (conf1 f g (1+ (car a)) (1+ (car b))) v) 991 (setq g (mapcar #'(lambda (a b) 992 (declare (special v)) 993 (- a (* v b))) 994 f (setq f g))) 995 (setq c (cons v c)) 996 (go a)) 997 (cond ((< (abs (- (conf1 f g (1+ (car a)) (car b)) v)) 998 (abs (- w v))) 999 (cond ((setq v (cdr b)) 1000 (setq f (conf6 f b)) 1001 (setq g (conf6 g b)) 1002 (setq b v)) 1003 (t (setq f (conf7 f b)) (setq g (conf7 g b))))) 1004 (t 1005 (cond ((setq v (cdr a)) 1006 (setq f (conf4 f a)) 1007 (setq g (conf4 g a)) 1008 (setq a v)) 1009 (t (setq f (conf5 f a)) (setq g (conf5 g a)))))) 1010 (go a))) 1011 1012(defun conf4 (n a) ;n_0*a_0+n_2,n_1*a_0+n_3,n_0,n_1 1013 (list (+ (* (car n) (car a)) (caddr n)) 1014 (+ (* (cadr n) (car a)) (cadddr n)) 1015 (car n) 1016 (cadr n))) 1017 1018(defun conf5 (n a) ;0,0, n_0*a_0,n_2 1019 (list 0 0 1020 (+ (* (car n) (car a)) (caddr n)) 1021 (+ (* (cadr n) (car a)) (cadddr n)))) 1022 1023(defun conf6 (n b) 1024 (list (+ (* (car n) (car b)) (cadr n)) 1025 (car n) 1026 (+ (* (caddr n) (car b)) (cadddr n)) 1027 (caddr n))) 1028 1029(defun conf7 (n b) 1030 (list 0 (+ (* (car n) (car b)) (cadr n)) 1031 0 (+ (* (caddr n) (car b)) (cadddr n)))) 1032 1033(defun cfsqrt (n) 1034 (cond ((cddr n) ;A non integer 1035 (merror (intl:gettext "cf: argument of sqrt must be an integer; found ~M") n)) 1036 ((setq n (cadr n)))) 1037 (setq n (sqcont n)) 1038 (cond ((= $cflength 1) 1039 (cons '(mlist simp) n)) 1040 ((do ((i 2 (1+ i)) 1041 (a (copy-tree (cdr n)))) 1042 ((> i $cflength) (cons '(mlist simp) n)) 1043 (setq n (nconc n (copy-tree a))))))) 1044 1045(defmfun $qunit (n) 1046 (let ((isqrtn ($isqrt n))) 1047 (when (or (not (integerp n)) 1048 (minusp n) 1049 (= (* isqrtn isqrtn) n)) 1050 (merror 1051 (intl:gettext "qunit: Argument must be a positive non quadratic integer."))) 1052 (let ((l (sqcont n))) 1053 (list '(mplus) (pelso1 l 0 1) 1054 (list '(mtimes) 1055 (list '(mexpt) n '((rat) 1 2)) 1056 (pelso1 l 1 0)))))) 1057 1058(defun pelso1 (l a b) 1059 (do ((i l (cdr i))) (nil) 1060 (and (null (cdr i)) (return b)) 1061 (setq b (+ a (* (car i) (setq a b)))))) 1062 1063(defun sqcont (n) 1064 (prog (q q1 q2 m m1 a0 a l) 1065 (setq a0 ($isqrt n) a (list a0) q2 1 m1 a0 1066 q1 (- n (* m1 m1)) l (* 2 a0)) 1067 a (setq a (cons (truncate (+ m1 a0) q1) a)) 1068 (cond ((equal (car a) l) 1069 (return (nreverse a)))) 1070 (setq m (- (* (car a) q1) m1) 1071 q (+ q2 (* (car a) (- m1 m))) 1072 q2 q1 q1 q m1 m) 1073 (go a))) 1074 1075(defun ratcf (x y) 1076 (prog (a b) 1077 a (cond ((equal y 1) (return (nreverse (cons x a)))) 1078 ((minusp x) 1079 (setq b (+ y (rem x y)) 1080 a (cons (1- (truncate x y)) a) 1081 x y y b)) 1082 ((> y x) 1083 (setq a (cons 0 a)) 1084 (setq b x x y y b)) 1085 ((equal x y) (return (nreverse (cons 1 a)))) 1086 ((setq b (rem x y)) 1087 (setq a (cons (truncate x y) a) x y y b))) 1088 (go a))) 1089 1090(defmfun $cfexpand (x) 1091 (cond ((null ($listp x)) x) 1092 ((cons '($matrix) (cfexpand (cdr x)))))) 1093 1094(defun cfexpand (ll) 1095 (do ((p1 0 p2) 1096 (p2 1 (simplify (list '(mplus) (list '(mtimes) (car l) p2) p1))) 1097 (q1 1 q2) 1098 (q2 0 (simplify (list '(mplus) (list '(mtimes) (car l) q2) q1))) 1099 (l ll (cdr l))) 1100 ((null l) (list (list '(mlist) p2 p1) (list '(mlist) q2 q1))))) 1101 1102;; Summation stuff 1103 1104(defun adsum (e) 1105 (push (simplify e) sum)) 1106 1107(defun adusum (e) 1108 (push (simplify e) usum)) 1109 1110(defun simpsum2 (exp i lo hi) 1111 (prog (*plus *times $simpsum u) 1112 (setq *plus (list 0) *times 1) 1113 (when (or (and (eq hi '$inf) (eq lo '$minf)) 1114 (equal 0 (m+ hi lo))) 1115 (setq $simpsum t lo 0) 1116 (setq *plus (cons (m* -1 *times (maxima-substitute 0 i exp)) *plus)) 1117 (setq exp (m+ exp (maxima-substitute (m- i) i exp)))) 1118 (cond ((eq ($sign (setq u (m- hi lo))) '$neg) 1119 (if (equal u -1) 1120 (return 0) 1121 (merror (intl:gettext "sum: lower bound ~M greater than upper bound ~M") lo hi))) 1122 ((free exp i) 1123 (return (m+l (cons (freesum exp lo hi *times) *plus)))) 1124 1125 ((setq exp (sumsum exp i lo hi)) 1126 (setq exp (m* *times (dosum (cadr exp) (caddr exp) 1127 (cadddr exp) (cadr (cdddr exp)) t :evaluate-summand nil)))) 1128 (t (return (m+l *plus)))) 1129 (return (m+l (cons exp *plus))))) 1130 1131(defun sumsum (e *var* lo hi) 1132 (let (sum usum) 1133 (cond ((eq hi '$inf) 1134 (cond (*infsumsimp (isum e lo)) 1135 ((setq usum (list e))))) 1136 ((finite-sum e 1 lo hi))) 1137 (cond ((eq sum nil) 1138 (return-from sumsum (list '(%sum) e *var* lo hi)))) 1139 (setq *plus 1140 (nconc (mapcar 1141 #'(lambda (q) (simptimes (list '(mtimes) *times q) 1 nil)) 1142 sum) 1143 *plus)) 1144 (and usum (setq usum (list '(%sum) (simplus (cons '(plus) usum) 1 t) *var* lo hi))))) 1145 1146(defun finite-sum (e y lo hi) 1147 (cond ((null e)) 1148 ((free e *var*) 1149 (adsum (m* y e (m+ hi 1 (m- lo))))) 1150 ((poly? e *var*) 1151 (adsum (m* y (fpolysum e lo hi)))) 1152 ((eq (caar e) '%binomial) (fbino e y lo hi)) 1153 ((eq (caar e) 'mplus) 1154 (mapc #'(lambda (q) (finite-sum q y lo hi)) (cdr e))) 1155 ((and (or (mtimesp e) (mexptp e) (mplusp e)) 1156 (fsgeo e y lo hi))) 1157 (t 1158 (adusum e) 1159 nil))) 1160 1161(defun isum-giveup (e) 1162 (cond ((atom e) nil) 1163 ((eq (caar e) 'mexpt) 1164 (not (or (free (cadr e) *var*) 1165 (ratp (caddr e) *var*)))) 1166 ((member (caar e) '(mplus mtimes) :test #'eq) 1167 (some #'identity (mapcar #'isum-giveup (cdr e)))) 1168 (t))) 1169 1170(defun isum (e lo) 1171 (cond ((isum-giveup e) 1172 (setq sum nil usum (list e))) 1173 ((eq (catch 'isumout (isum1 e lo)) 'divergent) 1174 (merror (intl:gettext "sum: sum is divergent."))))) 1175 1176(defun isum1 (e lo) 1177 (cond ((free e *var*) 1178 (unless (eq (asksign e) '$zero) 1179 (throw 'isumout 'divergent))) 1180 ((ratp e *var*) 1181 (adsum (ipolysum e lo))) 1182 ((eq (caar e) 'mplus) 1183 (mapc #'(lambda (x) (isum1 x lo)) (cdr e))) 1184 ( (isgeo e lo)) 1185 ((adusum e)))) 1186 1187(defun ipolysum (e lo) 1188 (ipoly1 ($expand e) lo)) 1189 1190(defun ipoly1 (e lo) 1191 (cond ((smono e *var*) 1192 (ipoly2 *a *n lo (asksign (simplify (list '(mplus) *n 1))))) 1193 ((mplusp e) 1194 (cons '(mplus) (mapcar #'(lambda (x) (ipoly1 x lo)) (cdr e)))) 1195 (t (adusum e) 1196 0))) 1197 1198(defun ipoly2 (a n lo sign) 1199 (cond ((member (asksign lo) '($zero $negative) :test #'eq) 1200 (throw 'isumout 'divergent))) 1201 (unless (equal lo 1) 1202 (let (($simpsum t)) 1203 (adsum `((%sum) 1204 ((mtimes) ,a -1 ((mexpt) ,*var* ,n)) 1205 ,*var* 1 ((mplus) -1 ,lo))))) 1206 (cond ((eq sign '$negative) 1207 (list '(mtimes) a ($zeta (meval (list '(mtimes) -1 n))))) 1208 ((throw 'isumout 'divergent)))) 1209 1210(defun fsgeo (e y lo hi) 1211 (let ((r ($ratsimp (div* (maxima-substitute (list '(mplus) *var* 1) *var* e) e)))) 1212 (cond ((equal r 1) 1213 (adsum 1214 (list '(mtimes) 1215 (list '(mplus) 1 hi (list '(mtimes) -1 lo)) 1216 (maxima-substitute lo *var* e)))) 1217 ((free r *var*) 1218 (adsum 1219 (list '(mtimes) y 1220 (maxima-substitute 0 *var* e) 1221 (list '(mplus) 1222 (list '(mexpt) r (list '(mplus) hi 1)) 1223 (list '(mtimes) -1 (list '(mexpt) r lo))) 1224 (list '(mexpt) (list '(mplus) r -1) -1))))))) 1225 1226(defun isgeo (e lo) 1227 (let ((r ($ratsimp (div* (maxima-substitute (list '(mplus) *var* 1) *var* e) e)))) 1228 (and (free r *var*) 1229 (isgeo1 (maxima-substitute lo *var* e) 1230 r (asksign (simplify (list '(mplus) (list '(mabs) r) -1))))))) 1231 1232(defun isgeo1 (a r sign) 1233 (cond ((eq sign '$positive) 1234 (throw 'isumout 'divergent)) 1235 ((eq sign '$zero) 1236 (throw 'isumout 'divergent)) 1237 ((eq sign '$negative) 1238 (adsum (list '(mtimes) a 1239 (list '(mexpt) (list '(mplus) 1 (list '(mtimes) -1 r)) -1)))))) 1240 1241 1242;; Sums of polynomials using 1243;; bernpoly(x+1, n) - bernpoly(x, n) = n*x^(n-1) 1244;; which implies 1245;; sum(k^n, k, A, B) = 1/(n+1)*(bernpoly(B+1, n+1) - bernpoly(A, n+1)) 1246;; 1247;; fpoly1 returns 1/(n+1)*(bernpoly(foo+1, n+1) - bernpoly(0, n+1)) for each power 1248;; in the polynomial e 1249 1250(defun fpolysum (e lo hi) ;returns *ans* 1251 (let ((a (fpoly1 (setq e ($expand ($ratdisrep ($rat e *var*)))) lo)) 1252 ($prederror)) 1253 (cond ((null a) 0) 1254 ((member lo '(0 1)) 1255 (maxima-substitute hi 'foo a)) 1256 (t 1257 (list '(mplus) (maxima-substitute hi 'foo a) 1258 (list '(mtimes) -1 (maxima-substitute (list '(mplus) lo -1) 'foo a))))))) 1259 1260(defun fpoly1 (e lo) 1261 (cond ((smono e *var*) 1262 (fpoly2 *a *n e lo)) 1263 ((eq (caar e) 'mplus) 1264 (cons '(mplus) (mapcar #'(lambda (x) (fpoly1 x lo)) (cdr e)))) 1265 (t (adusum e) 0))) 1266 1267(defun fpoly2 (a n e lo) 1268 (cond ((null (and (integerp n) (> n -1))) (adusum e) 0) 1269 ((equal n 0) 1270 (m* (cond ((signp e lo) 1271 (m1+ 'foo)) 1272 (t 'foo)) 1273 a)) 1274 (($ratsimp 1275 (m* a (list '(rat) 1 (1+ n)) 1276 (m- ($bernpoly (m+ 'foo 1) (1+ n)) 1277 ($bern (1+ n)))))))) 1278 1279;; fbino can do these sums: 1280;; a) sum(binomial(n,k),k,0,n) -> 2^n 1281;; b) sum(binomial(n-k,k,k,0,n) -> fib(n+1) 1282;; c) sum(binomial(n,2k),k,0,n) -> 2^(n-1) 1283;; d) sum(binomial(a+k,b),k,l,h) -> binomial(h+a+1,b+1) - binomial(l+a,b+1) 1284(defun fbino (e y lo hi) 1285 ;; e=binomial(n,d) 1286 (prog (n d l h) 1287 ;; check that n and d are linear in *var* 1288 (when (null (setq n (m2 (cadr e) (list 'n 'linear* *var*)))) 1289 (return (adusum e))) 1290 (setq n (cdr (assoc 'n n :test #'eq))) 1291 (when (null (setq d (m2 (caddr e) (list 'd 'linear* *var*)))) 1292 (return (adusum e))) 1293 (setq d (cdr (assoc 'd d :test #'eq))) 1294 1295 ;; binomial(a+b*k,c+b*k) -> binomial(a+b*k, a-c) 1296 (when (equal (cdr n) (cdr d)) 1297 (setq d (cons (m- (car n) (car d)) 0))) 1298 1299 (cond 1300 ;; substitute k with -k in sum(binomial(a+b*k, c-d*k)) 1301 ;; and sum(binomial(a-b*k,c)) 1302 ((and (numberp (cdr d)) 1303 (or (minusp (cdr d)) 1304 (and (zerop (cdr d)) 1305 (numberp (cdr n)) 1306 (minusp (cdr n))))) 1307 (rplacd d (- (cdr d))) 1308 (rplacd n (- (cdr n))) 1309 (setq l (m- hi) 1310 h (m- lo))) 1311 (t (setq l lo h hi))) 1312 1313 (cond 1314 1315 ;; sum(binomial(a+k,c),k,l,h) 1316 ((and (equal 0 (cdr d)) (equal 1 (cdr n))) 1317 (adsum (m* y (m- (list '(%binomial) (m+ h (car n) 1) (m+ (car d) 1)) 1318 (list '(%binomial) (m+ l (car n)) (m+ (car d) 1)))))) 1319 1320 ;; sum(binomial(n,k),k,0,n)=2^n 1321 ((and (equal 1 (cdr d)) (equal 0 (cdr n))) 1322 ;; sum(binomial(n,k+c),k,l,h)=sum(binomial(n,k+c+l),k,0,h-l) 1323 (let ((h1 (m- h l)) 1324 (c (m+ (car d) l))) 1325 (if (and (integerp (m- (car n) h1)) 1326 (integerp c)) 1327 (progn 1328 (adsum (m* y (m^ 2 (car n)))) 1329 (when (member (asksign (m- (m+ h1 c) (car n))) '($zero $negative) :test #'eq) 1330 (adsum (m* -1 y (dosum (list '(%binomial) (car n) *var*) 1331 *var* (m+ h1 c 1) (car n) t :evaluate-summand nil)))) 1332 (when (> c 0) 1333 (adsum (m* -1 y (dosum (list '(%binomial) (car n) *var*) 1334 *var* 0 (m- c 1) t :evaluate-summand nil))))) 1335 (adusum e)))) 1336 1337 ;; sum(binomial(b-k,k),k,0,floor(b/2))=fib(b+1) 1338 ((and (equal -1 (cdr n)) (equal 1 (cdr d))) 1339 ;; sum(binomial(a-k,b+k),k,l,h)=sum(binomial(a+b-k,k),k,l+b,h+b) 1340 (let ((h1 (m+ h (car d))) 1341 (l1 (m+ l (car d))) 1342 (a1 (m+ (car n) (car d)))) 1343 ;; sum(binomial(a1-k,k),k,0,floor(a1/2))=fib(a1+1) 1344 ;; we only do sums with h>floor(a1/2) 1345 (if (and (integerp l1) 1346 (member (asksign (m- h1 (m// a1 2))) '($zero $positive) :test #'eq)) 1347 (progn 1348 (adsum (m* y ($fib (m+ a1 1)))) 1349 (when (> l1 0) 1350 (adsum (m* -1 y (dosum (list '(%binomial) (m- a1 *var*) *var*) 1351 *var* 0 (m- l1 1) t :evaluate-summand nil))))) 1352 (adusum e)))) 1353 1354 ;; sum(binomial(n,2*k),k,0,floor(n/2))=2^(n-1) 1355 ;; sum(binomial(n,2*k+1),k,0,floor((n-1)/2))=2^(n-1) 1356 ((and (equal 0 (cdr n)) (equal 2 (cdr d))) 1357 ;; sum(binomial(a,2*k+b),k,l,h)=sum(binomial(a,2*k),k,l+b/2,h+b/2), b even 1358 ;; sum(binomial(a,2*k+b),k,l,h)=sum(binomial(a,2*k+1),k,l+(b-1)/2,h+(b-1)/2), b odd 1359 (let ((a (car n)) 1360 (r1 (if (oddp (car d)) 1 0)) 1361 (l1 (if (oddp (car d)) 1362 (m+ l (truncate (1- (car d)) 2)) 1363 (m+ l (truncate (car d) 2))))) 1364 (when (and (integerp l1) 1365 (member (asksign (m- a hi)) '($zero $positive) :test #'eq)) 1366 (adsum (m* y (m^ 2 (m- a 1)))) 1367 (when (> l1 0) 1368 (adsum (m* -1 y (dosum (list '(%binomial) a (m+ *var* *var* r1)) 1369 *var* 0 (m- l1 1) t :evaluate-summand nil))))))) 1370 1371 ;; other sums we can't do 1372 (t 1373 (adusum e))))) 1374 1375;; product routines 1376 1377(defmspec $product (l) 1378 (setq l (cdr l)) 1379 (cond ((not (= (length l) 4)) (merror (intl:gettext "product: expected exactly four arguments."))) 1380 ((dosum (car l) (cadr l) (meval (caddr l)) (meval (cadddr l)) nil :evaluate-summand t)))) 1381 1382(declare-top (special $ratsimpexpons)) 1383 1384;; Is this guy actually looking at the value of its middle arg? 1385 1386(defun simpprod (x y z) 1387 (let (($ratsimpexpons t)) 1388 (cond ((equal y 1) 1389 (setq y (simplifya (cadr x) z))) 1390 ((setq y (simptimes (list '(mexpt) (cadr x) y) 1 z))))) 1391 (simpprod1 y (caddr x) 1392 (simplifya (cadddr x) z) 1393 (simplifya (cadr (cdddr x)) z))) 1394 1395(defmfun $taytorat (e) 1396 (cond ((mbagp e) (cons (car e) (mapcar #'$taytorat (cdr e)))) 1397 ((or (atom e) (not (member 'trunc (cdar e) :test #'eq))) (ratf e)) 1398 ((catch 'srrat (srrat e))) 1399 (t (ratf ($ratdisrep e))))) 1400 1401(defun srrat (e) 1402 (cons (list 'mrat 'simp (caddar e) (cadddr (car e))) 1403 (srrat2 (cdr e)))) 1404 1405(defun srrat2 (e) 1406 (if (pscoefp e) e (srrat3 (terms e) (gvar e)))) 1407 1408(defun srrat3 (l *var*) 1409 (cond ((null l) '(0 . 1)) 1410 ((null (=1 (cdr (le l)))) 1411 (throw 'srrat nil)) 1412 ((null (n-term l)) 1413 (rattimes (cons (list *var* (car (le l)) 1) 1) 1414 (srrat2 (lc l)) 1415 t)) 1416 ((ratplus 1417 (rattimes (cons (list *var* (car (le l)) 1) 1) 1418 (srrat2 (lc l)) 1419 t) 1420 (srrat3 (n-term l) *var*))))) 1421 1422 1423(declare-top (special $props *i)) 1424 1425(defmspec $deftaylor (l) 1426 (prog (fun series param op ops) 1427 a (when (null (setq l (cdr l))) (return (cons '(mlist) ops))) 1428 (setq fun (meval (car l)) series (meval (cadr l)) l (cdr l) param () ) 1429 (when (or (atom fun) 1430 (if (eq (caar fun) 'mqapply) 1431 (or (cdddr fun) ; must be one parameter 1432 (null (cddr fun)) ; must have exactly one 1433 (do ((subs (cdadr fun) (cdr subs))) 1434 ((null subs) 1435 (setq op (caaadr fun)) 1436 (when (cddr fun) 1437 (setq param (caddr fun))) 1438 '()) 1439 (unless (atom (car subs)) (return 't)))) 1440 (progn 1441 (setq op (caar fun)) 1442 (when (cdr fun) (setq param (cadr fun))) 1443 (or (and (zl-get op 'op) (not (eq op 'mfactorial))) 1444 (not (atom (cadr fun))) 1445 (not (= (length fun) 2)))))) 1446 (merror (intl:gettext "deftaylor: don't know how to handle this function: ~M") fun)) 1447 (when (zl-get op 'sp2) 1448 (mtell (intl:gettext "deftaylor: redefining ~:M.~%") op)) 1449 (when param (setq series (subst 'sp2var param series))) 1450 (setq series (subsum '*index series)) 1451 (putprop op series 'sp2) 1452 (when (eq (caar fun) 'mqapply) 1453 (putprop op (cdadr fun) 'sp2subs)) 1454 (add2lnc op $props) 1455 (push op ops) 1456 (go a))) 1457 1458(defun subsum (*i e) (susum1 e)) 1459 1460(defun susum1 (e) 1461 (cond ((atom e) e) 1462 ((eq (caar e) '%sum) 1463 (if (null (smonop (cadr e) 'sp2var)) 1464 (merror (intl:gettext "deftaylor: argument must be a power series at 0.")) 1465 (subst *i (caddr e) e))) 1466 (t (recur-apply #'susum1 e)))) 1467 1468(declare-top (special varlist genvar $factorflag $ratfac *p* *var* *x*)) 1469 1470(defmfun $polydecomp (e v) 1471 (let ((varlist (list v)) 1472 (genvar nil) 1473 *var* p den $factorflag $ratfac) 1474 (setq p (cdr (ratf (ratdisrep e))) 1475 *var* (cdr (ratf v))) 1476 (cond ((or (null (cdr *var*)) 1477 (null (equal (cdar *var*) '(1 1)))) 1478 (merror (intl:gettext "polydecomp: second argument must be an atom; found ~M") v)) 1479 (t (setq *var* (caar *var*)))) 1480 (cond ((or (pcoefp (cdr p)) 1481 (null (eq (cadr p) *var*))) 1482 (setq den (cdr p) 1483 p (car p))) 1484 (t (merror (intl:gettext "polydecomp: cannot apply 'polydecomp' to a rational function.")))) 1485 (cons '(mlist) 1486 (cond ((or (pcoefp p) 1487 (null (eq (car p) *var*))) 1488 (list (rdis (cons p den)))) 1489 (t (setq p (pdecomp p *var*)) 1490 (do ((l 1491 (setq p (mapcar #'(lambda (q) (cons q 1)) p)) 1492 (cdr l)) 1493 (a)) 1494 ((null l) 1495 (cons (rdis (cons (caar p) 1496 (ptimes (cdar p) den))) 1497 (mapcar #'rdis (cdr p)))) 1498 (cond ((setq a (pdecpow (car l) *var*)) 1499 (rplaca l (car a)) 1500 (cond ((cdr l) 1501 (rplacd l 1502 (cons (ratplus 1503 (rattimes 1504 (cadr l) 1505 (cons (ptterm (cdaadr a) 1) 1506 (cdadr a)) 1507 t) 1508 (cons 1509 (ptterm (cdaadr a) 0) 1510 (cdadr a))) 1511 (cddr l)))) 1512 ((equal (cadr a) 1513 (cons (list *var* 1 1) 1))) 1514 (t (rplacd l (list (cadr a))))))))))))) 1515 1516 1517;;; POLYDECOMP is like $POLYDECOMP except it takes a poly in *POLY* format (as 1518;;; defined in SOLVE) (numerator of a RAT form) and returns a list of 1519;;; headerless rat forms. In otherwords, it is $POLYDECOMP minus type checking 1520;;; and conversions to/from general representation which SOLVE doesn't 1521;;; want/need on a general basis. 1522;;; It is used in the SOLVE package and as such it should have an autoload 1523;;; property 1524 1525(defun polydecomp (p *var*) 1526 (let ($factorflag $ratfac) 1527 (cond ((or (pcoefp p) 1528 (null (eq (car p) *var*))) 1529 (cons p nil)) 1530 (t (setq p (pdecomp p *var*)) 1531 (do ((l (setq p (mapcar #'(lambda (q) (cons q 1)) p)) 1532 (cdr l)) 1533 (a)) 1534 ((null l) 1535 (cons (cons (caar p) 1536 (cdar p)) 1537 (cdr p))) 1538 (cond ((setq a (pdecpow (car l) *var*)) 1539 (rplaca l (car a)) 1540 (cond ((cdr l) 1541 (rplacd l 1542 (cons (ratplus 1543 (rattimes 1544 (cadr l) 1545 (cons (ptterm (cdaadr a) 1) 1546 (cdadr a)) 1547 t) 1548 (cons 1549 (ptterm (cdaadr a) 0) 1550 (cdadr a))) 1551 (cddr l)))) 1552 ((equal (cadr a) 1553 (cons (list *var* 1 1) 1))) 1554 (t (rplacd l (list (cadr a)))))))))))) 1555 1556 1557 1558(defun pdecred (f h *var*) ;f = g(h(*var*)) 1559 (cond ((or (pcoefp h) (null (eq (car h) *var*)) 1560 (equal (cadr h) 1) 1561 (null (zerop (rem (cadr f) (cadr h)))) 1562 (and (null (pzerop (caadr (setq f (pdivide f h))))) 1563 (equal (cdadr f) 1))) 1564 nil) 1565 (t (do ((q (pdivide (caar f) h) (pdivide (caar q) h)) 1566 (i 1 (1+ i)) 1567 (*ans*)) 1568 ((pzerop (caar q)) 1569 (cond ((and (equal (cdadr q) 1) 1570 (or (pcoefp (caadr q)) 1571 (null (eq (caar (cadr q)) *var*)))) 1572 (psimp *var* (cons i (cons (caadr q) *ans*)))))) 1573 (cond ((and (equal (cdadr q) 1) 1574 (or (pcoefp (caadr q)) 1575 (null (eq (caar (cadr q)) *var*)))) 1576 (and (null (pzerop (caadr q))) 1577 (setq *ans* (cons i (cons (caadr q) *ans*))))) 1578 (t (return nil))))))) 1579 1580(defun pdecomp (p *var*) 1581 (let ((c (ptterm (cdr p) 0)) 1582 (a) (*x* (list *var* 1 1))) 1583 (cons (pcplus c (car (setq a (pdecomp* (pdifference p c))))) 1584 (cdr a)))) 1585 1586(defun pdecomp* (*p*) 1587 (let ((a) 1588 (l (pdecgdfrm (pfactor (pquotient *p* *x*))))) 1589 (cond ((or (pdecprimep (cadr *p*)) 1590 (null (setq a (pdecomp1 *x* l)))) 1591 (list *p*)) 1592 (t (append (pdecomp* (car a)) (cdr a)))))) 1593 1594(defun pdecomp1 (prod l) 1595 (cond ((null l) 1596 (and (null (equal (cadr prod) (cadr *p*))) 1597 (setq l (pdecred *p* prod *var*)) 1598 (list l prod))) 1599 ((pdecomp1 prod (cdr l))) 1600 (t (pdecomp1 (ptimes (car l) prod) (cdr l))))) 1601 1602(defun pdecgdfrm (l) ;Get list of divisors 1603 (do ((l (copy-list l )) 1604 (ll (list (car l)) 1605 (cons (car l) ll))) 1606 (nil) 1607 (rplaca (cdr l) (1- (cadr l))) 1608 (cond ((signp e (cadr l)) 1609 (setq l (cddr l)))) 1610 (cond ((null l) (return ll))))) 1611 1612(defun pdecprimep (x) 1613 (setq x (cfactorw x)) 1614 (and (null (cddr x)) (equal (cadr x) 1))) 1615 1616(defun pdecpow (p *var*) 1617 (setq p (car p)) 1618 (let ((p1 (pderivative p *var*)) 1619 p2 p1p p1c a lin p2p) 1620 (setq p1p (oldcontent p1) 1621 p1c (car p1p) p1p (cadr p1p)) 1622 (setq p2 (pderivative p1 *var*)) 1623 (setq p2p (cadr (oldcontent p2))) 1624 (and (setq lin (testdivide p1p p2p)) 1625 (null (pcoefp lin)) 1626 (eq (car lin) *var*) 1627 (list (ratplus 1628 (rattimes (cons (list *var* (cadr p) 1) 1) 1629 (setq a (ratreduce p1c 1630 (ptimes (cadr p) 1631 (caddr lin)))) 1632 t) 1633 (ratdif (cons p 1) 1634 (rattimes a (cons (pexpt lin (cadr p)) 1) 1635 t))) 1636 (cons lin 1))))) 1637 1638(declare-top (unspecial *mfactl *factlist donel nn* dn* 1639 *var* *ans* *n *a* 1640 *infsumsimp *times *plus sum usum makef)) 1641