?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fix the inefficient evaluation of translated predicates
[maxima.git] / src / combin.lisp
blob62aebb7bba95ac00a7670ef1d89e56bbfc85dbe9
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 ;;; The Riemann Zeta function has mirror symmetry
523
524 (defprop %zeta t commutes-with-conjugate)
525
526 ;;; The Riemann Zeta function distributes over lists, matrices, and equations
527
528 (defprop %zeta (mlist $matrix mequal) distribute_over)
529
530 ;;; We support a simplim%function. The function is looked up in simplimit and
531 ;;; handles specific values of the function.
532
533 (defprop %zeta simplim%zeta simplim%function)
534
535 (defun simplim%zeta (expr var val)
536 ;; Look for the limit of the argument
537 (let* ((arg (limit (cadr expr) var val 'think))
538 (dir (limit (add (cadr expr) (neg arg)) var val 'think)))
539 (cond
540 ;; Handle an argument 1 at this place
541 ((onep1 arg)
542 (cond ((eq dir '$zeroa)
543 '$inf)
544 ((eq dir '$zerob)
545 '$minf)
546 (t '$infinity)))
547 (t
548 ;; All other cases are handled by the simplifier of the function.
549 (simplify (list '(%zeta) arg))))))
550
551 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
552
553 (def-simplifier zeta (z)
554 (cond
555
556 ;; Check for special values
557 ((eq z '$inf) 1)
558 ((alike1 z '((mtimes) -1 $minf)) 1)
559 ((zerop1 z)
560 (cond (($bfloatp z) ($bfloat '((rat) -1 2)))
561 ((floatp z) -0.5)
562 (t '((rat simp) -1 2))))
563 ((onep1 z)
564 (simp-domain-error (intl:gettext "zeta: zeta(~:M) is undefined.") z))
565
566 ;; Check for numerical evaluation
567 ((or (bigfloat-numerical-eval-p z)
568 (complex-bigfloat-numerical-eval-p z)
569 (float-numerical-eval-p z)
570 (complex-float-numerical-eval-p z))
571 (to (float-zeta z)))
572 ;; Check for transformations and argument simplifications
573 ((integerp z)
574 (cond
575 ((oddp z)
576 (cond ((> z 1)
577 (give-up))
578 ((setq z (sub 1 z))
579 (mul -1 (div ($bern z) z)))))
580 ((minusp z) 0)
581 ((not $zeta%pi)
582 (give-up))
583 (t (let ($numer $float)
584 (mul (power '$%pi z)
585 (mul (div (power 2 (1- z))
586 (take '(mfactorial) z))
587 (take '(mabs) ($bern z))))))))
588 ((and (mnegp z))
589 ;; z is a negative number. We can use the relationship (A&S
590 ;; 23.2.6):
591 ;;
592 ;; zeta(s) = 2^s*%pi^(s-1)*sin(%pi/2*s)*gamma(1-s)*zeta(1-s)
593 ;;
594 ;; to transform zeta of a negative number in terms of a zeta to a
595 ;; positive number.
596 (mul (power 2 z)
597 (power '$%pi (sub z 1))
598 (ftake '%sin (div (mul '$%pi z) 2))
599 (ftake '%gamma (sub 1 z))
600 (ftake '%zeta (sub 1 z))))
601 (t
602 (give-up))))
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 (tau (bigfloat:imagpart s)))
647 (cond
648 ;; abs(s)^2 < epsilon, use the expansion zeta(s) = -1/2-1/2*log(2*%pi)*s
649 ((bigfloat:< (bigfloat:abs (bigfloat:* s s)) (bigfloat:epsilon s))
650 (bigfloat:+ -1/2
651 (bigfloat:* -1/2
652 (bigfloat:log (bigfloat:* 2 (bigfloat:%pi s)))
653 s)))
654
655 ;; Reflection formula:
656 ;; zeta(s) = 2^s*%pi^(s-1)*sin(%pi*s/2)*gamma(1-s)*zeta(1-s)
657 ((not (bigfloat:plusp sigma))
658 (let ((n (bigfloat:floor sigma)))
659 ;; If s is a negative even integer, zeta(s) is zero,
660 ;; from the reflection formula because sin(%pi*s/2) is 0.
661 (when (and (bigfloat:zerop tau) (bigfloat:= n sigma) (evenp n))
662 (return-from float-zeta (bigfloat:float 0.0 sigma))))
663 (bigfloat:* (bigfloat:expt 2 s)
664 (bigfloat:expt (bigfloat:%pi s)
665 (bigfloat:- s 1))
666 (bigfloat:sin (bigfloat:* (bigfloat:/ (bigfloat:%pi s)
667 2)
668 s))
669 (bigfloat:to ($gamma (to (bigfloat:- 1 s))))
670 (float-zeta (bigfloat:- 1 s))))
671
672 ;; The general formula from above. Call the imaginary part "tau" rather
673 ;; than the "t" above, because that isn't a CL keyword...
674 (t
675 (let* ((logh
676 (bigfloat:-
677 (if (bigfloat:zerop tau) 0
678 (bigfloat:+
679 (bigfloat:* 1.6 (bigfloat:abs tau))
680 (bigfloat:log (bigfloat:1+
681 (bigfloat:abs
682 (bigfloat:/ tau sigma))))))
683 (bigfloat:log
684 (bigfloat:abs
685 (bigfloat:- 1 (bigfloat:expt 2 (bigfloat:- 1 s)))))))
686
687 (logeps (bigfloat:log (bigfloat:epsilon s)))
688
689 (n (max (bigfloat:ceiling
690 (bigfloat:/ (bigfloat:- logh logeps) (bigfloat:log 8)))
691 1))
692
693 (sum1 0)
694 (sum2 0))
695 (flet ((binsum (k n)
696 ;; sum(binomial(n,j), j, k, n) = sum(binomial(n,j), j, n, k)
697 (let ((sum 0)
698 (term 1))
699 (loop for j from n downto k
700 do
701 (progn
702 (incf sum term)
703 (setf term (/ (* term j) (+ n 1 (- j))))))
704 sum)))
705 ;; (format t "n = ~D terms~%" n)
706 ;; sum1 = sum((-1)^(k-1)/k^s,k,1,n)
707 ;; sum2 = sum((-1)^(k-1)/e(n,k-n)/k^s, k, n+1, 2*n)
708 ;; = (-1)^n*sum((-1)^(m-1)*e(n,m)/(n+k)^s, k, 1, n)
709 (loop for k from 1 to n
710 for d = (bigfloat:expt k s)
711 do (progn
712 (bigfloat:incf sum1 (bigfloat:/ (cl:expt -1 (- k 1)) d))
713 (bigfloat:incf sum2 (bigfloat:/ (* (cl:expt -1 (- k 1))
714 (binsum k n))
715 (bigfloat:expt (+ k n) s))))
716 finally (return (values sum1 sum2)))
717 (when (oddp n)
718 (setq sum2 (bigfloat:- sum2)))
719 (bigfloat:/ (bigfloat:+ sum1
720 (bigfloat:/ sum2 (bigfloat:expt 2 n)))
721 (bigfloat:- 1 (bigfloat:expt 2 (bigfloat:- 1 s))))))))))
722
723 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
724
725 ;; Returns the N'th Fibonacci number.
726 (defmfun $fib (n)
727 (cond ((fixnump n)
728 (ffib n))
729 (t
730 (setq $prevfib `(($fib) ,(add2* n -1)))
731 `(($fib) ,n))))
732
733 ;; Main routine for computing the n'th Fibonacci number where n is an
734 ;; integer.
735 (defun ffib (%n)
736 (declare (fixnum %n))
737 (cond ((= %n -1)
738 (setq $prevfib -1)
739 1)
740 ((zerop %n)
741 (setq $prevfib 1)
742 0)
743 (t
744 (let* ((f2 (ffib (ash (logandc2 %n 1) -1))) ; f2 = fib(n/2) or fib((n-1)/2)
745 (x (+ f2 $prevfib))
746 (y (* $prevfib $prevfib))
747 (z (* f2 f2)))
748 (setq f2 (- (* x x) y)
749 $prevfib (+ y z))
750 (when (oddp %n)
751 (psetq $prevfib f2
752 f2 (+ f2 $prevfib)))
753 f2))))
754
755 ;; Returns the N'th Lucas number defined by the following recursion,
756 ;; where L(n) is the n'th Lucas number:
757 ;;
758 ;; L(0) = 2;
759 ;; L(1) = 1;
760 ;; L(n) = L(n-1) + L(n-2), n > 1;
761 (defmfun $lucas (n)
762 (cond
763 ((fixnump n)
764 (lucas n))
765 (t
766 (setq $next_lucas `(($lucas) ,(add2* n 1)))
767 `(($lucas) ,n))))
768
769 (defun lucas (n)
770 (declare (fixnum n))
771 (let ((w 2) (x 2) (y 1) u v (sign (signum n))) (declare (fixnum sign))
772 (setq n (abs n))
773 (do ((i (1- (integer-length n)) (1- i)))
774 ((< i 0))
775 (declare (fixnum i))
776 (setq u (* x x) v (* y y))
777 (if (logbitp i n)
778 (setq y (+ v w) x (+ y (- u) w) w -2)
779 (setq x (- u w) y (+ v w (- x)) w 2) ))
780 (cond
781 ((or (= 1 sign) (not (logbitp 0 n)))
782 (setq $next_lucas y)
783 x )
784 (t
785 (setq $next_lucas (neg y))
786 (neg x) ))))
787
788 ;; continued fraction stuff
789
790 (defmfun $cfdisrep (a)
791 (cond ((not ($listp a))
792 (merror (intl:gettext "cfdisrep: argument must be a list; found ~M") a))
793 ((null (cddr a)) (cadr a))
794 ((equal (cadr a) 0)
795 (list '(mexpt) (cfdisrep1 (cddr a)) -1))
796 ((cfdisrep1 (cdr a)))))
797
798 (defun cfdisrep1 (a)
799 (cond ((cdr a)
800 (list '(mplus simp cf) (car a)
801 (prog2 (setq a (cfdisrep1 (cdr a)))
802 (cond ((integerp a) (list '(rat simp) 1 a))
803 (t (list '(mexpt simp) a -1))))))
804 ((car a))))
805
806 (defun cfmak (a)
807 (setq a (meval a))
808 (cond ((integerp a) (list a))
809 ((eq (caar a) 'mlist) (cdr a))
810 ((eq (caar a) 'rat) (ratcf (cadr a) (caddr a)))
811 ((merror (intl:gettext "cf: continued fractions must be lists or integers; found ~M") a))))
812
813 (defun makcf (a)
814 (cond ((null (cdr a)) (car a))
815 ((cons '(mlist simp cf) a))))
816
817 ;;; Translation properties for $CF defined in MAXSRC;TRANS5 >
818
819 (defmspec $cf (a)
820 (cfratsimp (let ($listarith)
821 (cfeval (meval (fexprcheck a))))))
822
823 ;; Definition of cfratsimp as given in SF bug report # 620928.
824 (defun cfratsimp (a)
825 (cond ((atom a) a)
826 ((member 'cf (car a) :test #'eq) a)
827 (t (cons '(mlist cf simp)
828 (apply 'find-cf (cf-back-recurrence (cdr a)))))))
829
830 ; Code to expand nth degree roots of integers into continued fraction
831 ; approximations. E.g. cf(2^(1/3))
832 ; Courtesy of Andrei Zorine (feniy@mail.nnov.ru) 2005/05/07
833
834 (defun cfnroot(b)
835 (let ((ans (list '(mlist xf))) ent ($algebraic $true))
836 (dotimes (i $cflength (nreverse ans))
837 (setq ent (meval `(($floor) ,b))
838 ans (cons ent ans)
839 b ($ratsimp (m// (m- b ent)))))))
840
841 (defun cfeval (a)
842 (let (temp $ratprint)
843 (cond ((integerp a) (list '(mlist cf) a))
844 ((floatp a)
845 (let ((a (maxima-rationalize a)))
846 (cons '(mlist cf) (ratcf (car a) (cdr a)))))
847 (($bfloatp a)
848 (let (($bftorat t))
849 (setq a (bigfloat2rat a))
850 (cons '(mlist cf) (ratcf (car a) (cdr a)))))
851 ((atom a)
852 (merror (intl:gettext "cf: ~:M is not a continued fraction.") a))
853 ((eq (caar a) 'rat)
854 (cons '(mlist cf) (ratcf (cadr a) (caddr a))))
855 ((eq (caar a) 'mlist)
856 (cfratsimp a))
857 ;;the following doesn't work for non standard form
858 ;; (cfplus a '((mlist) 0)))
859 ((and (mtimesp a) (cddr a) (null (cdddr a))
860 (fixnump (cadr a))
861 (mexptp (caddr a))
862 (fixnump (cadr (caddr a)))
863 (alike1 (caddr (caddr a)) '((rat) 1 2)))
864 (cfsqrt (cfeval (* (expt (cadr a) 2) (cadr (caddr a))))))
865 ((eq (caar a) 'mexpt)
866 (cond ((alike1 (caddr a) '((rat) 1 2))
867 (cfsqrt (cfeval (cadr a))))
868 ((integerp (m* 2 (caddr a))) ; a^(n/2) was sqrt(a^n)
869 (cfsqrt (cfeval (cfexpt (cadr a) (m* 2 (caddr a))))))
870 ((integerp (cadr a)) (cfnroot a)) ; <=== new case x
871 ((cfexpt (cfeval (cadr a)) (caddr a)))))
872 ((setq temp (assoc (caar a) '((mplus . cfplus) (mtimes . cftimes) (mquotient . cfquot)
873 (mdifference . cfdiff) (mminus . cfminus)) :test #'eq))
874 (cf (cfeval (cadr a)) (cddr a) (cdr temp)))
875 ((eq (caar a) 'mrat)
876 (cfeval ($ratdisrep a)))
877 (t (merror (intl:gettext "cf: ~:M is not a continued fraction.") a)))))
878
879 (defun cf (a l fun)
880 (cond ((null l) a)
881 ((cf (funcall fun a (meval (list '($cf) (car l)))) (cdr l) fun))))
882
883 (defun cfplus (a b)
884 (setq a (cfmak a) b (cfmak b))
885 (makcf (cffun '(0 1 1 0) '(0 0 0 1) a b)))
886
887 (defun cftimes (a b)
888 (setq a (cfmak a) b (cfmak b))
889 (makcf (cffun '(1 0 0 0) '(0 0 0 1) a b)))
890
891 (defun cfdiff (a b)
892 (setq a (cfmak a) b (cfmak b))
893 (makcf (cffun '(0 1 -1 0) '(0 0 0 1) a b)))
894
895 (defun cfmin (a)
896 (setq a (cfmak a))
897 (makcf (cffun '(0 0 -1 0) '(0 0 0 1) a '(0))))
898
899 (defun cfquot (a b)
900 (setq a (cfmak a) b (cfmak b))
901 (makcf (cffun '(0 1 0 0) '(0 0 1 0) a b)))
902
903 (defun cfexpt (b e)
904 (setq b (cfmak b))
905 (cond ((null (integerp e))
906 (merror (intl:gettext "cf: can't raise continued fraction to non-integral power ~M") e))
907 ((let ((n (abs e)))
908 (do ((n (ash n -1) (ash n -1))
909 (s (cond ((oddp n) b)
910 (t '(1)))))
911 ((zerop n)
912 (makcf
913 (cond ((signp g e)
914 s)
915 ((cffun '(0 0 0 1) '(0 1 0 0) b '(1))))))
916 (setq b (cffun '(1 0 0 0) '(0 0 0 1) b b))
917 (and (oddp n)
918 (setq s (cffun '(1 0 0 0) '(0 0 0 1) s b))))))))
919
920
921 (defun conf1 (f g a b &aux (den (conf2 g a b)))
922 (cond ((zerop den)
923 (* (signum (conf2 f a b )) ; (/ most-positive-fixnum (^ 2 4))
924 #.(expt 2 31)))
925 (t (truncate (conf2 f a b) den))))
926
927 (defun conf2 (n a b) ;2*(abn_0+an_1+bn_2+n_3)
928 (* 2 (+ (* (car n) a b)
929 (* (cadr n) a)
930 (* (caddr n) b)
931 (cadddr n))))
932
933 ;;(cffun '(0 1 1 0) '(0 0 0 1) '(1 2) '(1 1 1 2)) gets error
934 ;;should give (3 10)
935
936 (defun cf-convergents-p-q (cf &optional (n (length cf)) &aux pp qq)
937 "returns two lists such that pp_i/qq_i is the quotient of the first i terms
938 of cf"
939 (case (length cf)
940 (0 1)
941 (1 cf(list 1))
942 (t
943 (setq pp (list (1+ (* (first cf) (second cf))) (car cf)))
944 (setq qq (list (second cf) 1))
945 (show pp qq)
946 (setq cf (cddr cf))
947 (loop for i from 2 to n
948 while cf
949 do
950 (push (+ (* (car cf) (car pp))
951 (second pp)) pp)
952 (push (+ (* (car cf) (car qq))
953 (second qq)) qq)
954 (setq cf (cdr cf))
955 finally (return (list (reverse pp) (reverse qq)))))))
956
957
958 (defun find-cf1 (p q so-far)
959 (multiple-value-bind (quot rem) (truncate p q)
960 (cond ((< rem 0) (incf rem q) (incf quot -1))
961 ((zerop rem) (return-from find-cf1 (cons quot so-far))))
962 (setq so-far (cons quot so-far))
963 (find-cf1 q rem so-far)))
964
965 (defun find-cf (p q)
966 "returns the continued fraction for p and q integers, q not zero"
967 (cond ((zerop q) (maxima-error "find-cf: quotient by zero"))
968 ((< q 0) (setq p (- p)) (setq q (- q))))
969 (nreverse (find-cf1 p q ())))
970
971 (defun cf-back-recurrence (cf &aux tem (num-gg 0)(den-gg 1))
972 "converts CF (a continued fraction list) to a list of numerator
973 denominator using recurrence from end
974 and not calculating intermediate quotients.
975 The numerator and denom are relatively
976 prime"
977 (loop for v in (reverse cf)
978 do (setq tem (* den-gg v))
979 (setq tem (+ tem num-gg))
980 (setq num-gg den-gg)
981 (setq den-gg tem)
982 finally
983 (return
984 (cond ((and (<= den-gg 0) (< num-gg 0))
985 (list (- den-gg) (- num-gg)))
986 (t(list den-gg num-gg))))))
987
988 (declare-top (unspecial w))
989
990 ;;(cffun '(0 1 1 0) '(0 0 0 1) '(1 2) '(1 1 1 2)) gets error
991 ;;should give (3 10)
992
993 (defun cffun (f g a b)
994 (prog (c v w)
995 (declare (special v))
996 a (and (zerop (cadddr g))
997 (zerop (caddr g))
998 (zerop (cadr g))
999 (zerop (car g))
1000 (return (reverse c)))
1001 (and (equal (setq w (conf1 f g (car a) (1+ (car b))))
1002 (setq v (conf1 f g (car a) (car b))))
1003 (equal (conf1 f g (1+ (car a)) (car b)) v)
1004 (equal (conf1 f g (1+ (car a)) (1+ (car b))) v)
1005 (setq g (mapcar #'(lambda (a b)
1006 (declare (special v))
1007 (- a (* v b)))
1008 f (setq f g)))
1009 (setq c (cons v c))
1010 (go a))
1011 (cond ((< (abs (- (conf1 f g (1+ (car a)) (car b)) v))
1012 (abs (- w v)))
1013 (cond ((setq v (cdr b))
1014 (setq f (conf6 f b))
1015 (setq g (conf6 g b))
1016 (setq b v))
1017 (t (setq f (conf7 f b)) (setq g (conf7 g b)))))
1018 (t
1019 (cond ((setq v (cdr a))
1020 (setq f (conf4 f a))
1021 (setq g (conf4 g a))
1022 (setq a v))
1023 (t (setq f (conf5 f a)) (setq g (conf5 g a))))))
1024 (go a)))
1025
1026 (defun conf4 (n a) ;n_0*a_0+n_2,n_1*a_0+n_3,n_0,n_1
1027 (list (+ (* (car n) (car a)) (caddr n))
1028 (+ (* (cadr n) (car a)) (cadddr n))
1029 (car n)
1030 (cadr n)))
1031
1032 (defun conf5 (n a) ;0,0, n_0*a_0,n_2
1033 (list 0 0
1034 (+ (* (car n) (car a)) (caddr n))
1035 (+ (* (cadr n) (car a)) (cadddr n))))
1036
1037 (defun conf6 (n b)
1038 (list (+ (* (car n) (car b)) (cadr n))
1039 (car n)
1040 (+ (* (caddr n) (car b)) (cadddr n))
1041 (caddr n)))
1042
1043 (defun conf7 (n b)
1044 (list 0 (+ (* (car n) (car b)) (cadr n))
1045 0 (+ (* (caddr n) (car b)) (cadddr n))))
1046
1047 (defun cfsqrt (n)
1048 (cond ((cddr n) ;A non integer
1049 (merror (intl:gettext "cf: argument of sqrt must be an integer; found ~M") n))
1050 ((setq n (cadr n))))
1051 (setq n (sqcont n))
1052 (cond ((= $cflength 1)
1053 (cons '(mlist simp) n))
1054 ((do ((i 2 (1+ i))
1055 (a (copy-tree (cdr n))))
1056 ((> i $cflength) (cons '(mlist simp) n))
1057 (setq n (nconc n (copy-tree a)))))))
1058
1059 (defmfun $qunit (n)
1060 (let ((isqrtn ($isqrt n)))
1061 (when (or (not (integerp n))
1062 (minusp n)
1063 (= (* isqrtn isqrtn) n))
1064 (merror
1065 (intl:gettext "qunit: Argument must be a positive non quadratic integer.")))
1066 (let ((l (sqcont n)))
1067 (list '(mplus) (pelso1 l 0 1)
1068 (list '(mtimes)
1069 (list '(mexpt) n '((rat) 1 2))
1070 (pelso1 l 1 0))))))
1071
1072 (defun pelso1 (l a b)
1073 (do ((i l (cdr i))) (nil)
1074 (and (null (cdr i)) (return b))
1075 (setq b (+ a (* (car i) (setq a b))))))
1076
1077 (defun sqcont (n)
1078 (prog (q q1 q2 m m1 a0 a l)
1079 (setq a0 ($isqrt n) a (list a0) q2 1 m1 a0
1080 q1 (- n (* m1 m1)) l (* 2 a0))
1081 a (setq a (cons (truncate (+ m1 a0) q1) a))
1082 (cond ((equal (car a) l)
1083 (return (nreverse a))))
1084 (setq m (- (* (car a) q1) m1)
1085 q (+ q2 (* (car a) (- m1 m)))
1086 q2 q1 q1 q m1 m)
1087 (go a)))
1088
1089 (defun ratcf (x y)
1090 (prog (a b)
1091 a (cond ((equal y 1) (return (nreverse (cons x a))))
1092 ((minusp x)
1093 (setq b (+ y (rem x y))
1094 a (cons (1- (truncate x y)) a)
1095 x y y b))
1096 ((> y x)
1097 (setq a (cons 0 a))
1098 (setq b x x y y b))
1099 ((equal x y) (return (nreverse (cons 1 a))))
1100 ((setq b (rem x y))
1101 (setq a (cons (truncate x y) a) x y y b)))
1102 (go a)))
1103
1104 (defmfun $cfexpand (x)
1105 (cond ((null ($listp x)) x)
1106 ((cons '($matrix) (cfexpand (cdr x))))))
1107
1108 (defun cfexpand (ll)
1109 (do ((p1 0 p2)
1110 (p2 1 (simplify (list '(mplus) (list '(mtimes) (car l) p2) p1)))
1111 (q1 1 q2)
1112 (q2 0 (simplify (list '(mplus) (list '(mtimes) (car l) q2) q1)))
1113 (l ll (cdr l)))
1114 ((null l) (list (list '(mlist) p2 p1) (list '(mlist) q2 q1)))))
1115
1116 ;; Summation stuff
1117
1118 (defun adsum (e)
1119 (push (simplify e) sum))
1120
1121 (defun adusum (e)
1122 (push (simplify e) usum))
1123
1124 (defun simpsum2 (exp i lo hi)
1125 (prog (*plus *times $simpsum u)
1126 (setq *plus (list 0) *times 1)
1127 (when (or (and (eq hi '$inf) (eq lo '$minf))
1128 (equal 0 (m+ hi lo)))
1129 (setq $simpsum t lo 0)
1130 (setq *plus (cons (m* -1 *times (maxima-substitute 0 i exp)) *plus))
1131 (setq exp (m+ exp (maxima-substitute (m- i) i exp))))
1132 (cond ((eq ($sign (setq u (m- hi lo))) '$neg)
1133 (if (equal u -1)
1134 (return 0)
1135 (merror (intl:gettext "sum: lower bound ~M greater than upper bound ~M") lo hi)))
1136 ((free exp i)
1137 (return (m+l (cons (freesum exp lo hi *times) *plus))))
1138
1139 ((setq exp (sumsum exp i lo hi))
1140 (setq exp (m* *times (dosum (cadr exp) (caddr exp)
1141 (cadddr exp) (cadr (cdddr exp)) t :evaluate-summand nil))))
1142 (t (return (m+l *plus))))
1143 (return (m+l (cons exp *plus)))))
1144
1145 (defun sumsum (e *var* lo hi)
1146 (let (sum usum)
1147 (cond ((eq hi '$inf)
1148 (cond (*infsumsimp (isum e lo))
1149 ((setq usum (list e)))))
1150 ((finite-sum e 1 lo hi)))
1151 (cond ((eq sum nil)
1152 (return-from sumsum (list '(%sum) e *var* lo hi))))
1153 (setq *plus
1154 (nconc (mapcar
1155 #'(lambda (q) (simptimes (list '(mtimes) *times q) 1 nil))
1156 sum)
1157 *plus))
1158 (and usum (setq usum (list '(%sum) (simplus (cons '(plus) usum) 1 t) *var* lo hi)))))
1159
1160 (defun finite-sum (e y lo hi)
1161 (cond ((null e))
1162 ((free e *var*)
1163 (adsum (m* y e (m+ hi 1 (m- lo)))))
1164 ((poly? e *var*)
1165 (adsum (m* y (fpolysum e lo hi))))
1166 ((eq (caar e) '%binomial) (fbino e y lo hi))
1167 ((eq (caar e) 'mplus)
1168 (mapc #'(lambda (q) (finite-sum q y lo hi)) (cdr e)))
1169 ((and (or (mtimesp e) (mexptp e) (mplusp e))
1170 (fsgeo e y lo hi)))
1171 (t
1172 (adusum e)
1173 nil)))
1174
1175 (defun isum-giveup (e)
1176 (cond ((atom e) nil)
1177 ((eq (caar e) 'mexpt)
1178 (not (or (free (cadr e) *var*)
1179 (ratp (caddr e) *var*))))
1180 ((member (caar e) '(mplus mtimes) :test #'eq)
1181 (some #'identity (mapcar #'isum-giveup (cdr e))))
1182 (t)))
1183
1184 (defun isum (e lo)
1185 (cond ((isum-giveup e)
1186 (setq sum nil usum (list e)))
1187 ((eq (catch 'isumout (isum1 e lo)) 'divergent)
1188 (merror (intl:gettext "sum: sum is divergent.")))))
1189
1190 (defun isum1 (e lo)
1191 (cond ((free e *var*)
1192 (unless (eq (asksign e) '$zero)
1193 (throw 'isumout 'divergent)))
1194 ((ratp e *var*)
1195 (adsum (ipolysum e lo)))
1196 ((eq (caar e) 'mplus)
1197 (mapc #'(lambda (x) (isum1 x lo)) (cdr e)))
1198 ( (isgeo e lo))
1199 ((adusum e))))
1200
1201 (defun ipolysum (e lo)
1202 (ipoly1 ($expand e) lo))
1203
1204 (defun ipoly1 (e lo)
1205 (cond ((smono e *var*)
1206 (ipoly2 *a *n lo (asksign (simplify (list '(mplus) *n 1)))))
1207 ((mplusp e)
1208 (cons '(mplus) (mapcar #'(lambda (x) (ipoly1 x lo)) (cdr e))))
1209 (t (adusum e)
1210 0)))
1211
1212 (defun ipoly2 (a n lo sign)
1213 (cond ((member (asksign lo) '($zero $negative) :test #'eq)
1214 (throw 'isumout 'divergent)))
1215 (unless (equal lo 1)
1216 (let (($simpsum t))
1217 (adsum `((%sum)
1218 ((mtimes) ,a -1 ((mexpt) ,*var* ,n))
1219 ,*var* 1 ((mplus) -1 ,lo)))))
1220 (cond ((eq sign '$negative)
1221 (list '(mtimes) a ($zeta (meval (list '(mtimes) -1 n)))))
1222 ((throw 'isumout 'divergent))))
1223
1224 (defun fsgeo (e y lo hi)
1225 (let ((r ($ratsimp (div* (maxima-substitute (list '(mplus) *var* 1) *var* e) e))))
1226 (cond ((equal r 1)
1227 (adsum
1228 (list '(mtimes)
1229 (list '(mplus) 1 hi (list '(mtimes) -1 lo))
1230 (maxima-substitute lo *var* e))))
1231 ((free r *var*)
1232 (adsum
1233 (list '(mtimes) y
1234 (maxima-substitute 0 *var* e)
1235 (list '(mplus)
1236 (list '(mexpt) r (list '(mplus) hi 1))
1237 (list '(mtimes) -1 (list '(mexpt) r lo)))
1238 (list '(mexpt) (list '(mplus) r -1) -1)))))))
1239
1240 (defun isgeo (e lo)
1241 (let ((r ($ratsimp (div* (maxima-substitute (list '(mplus) *var* 1) *var* e) e))))
1242 (and (free r *var*)
1243 (isgeo1 (maxima-substitute lo *var* e)
1244 r (asksign (simplify (list '(mplus) (list '(mabs) r) -1)))))))
1245
1246 (defun isgeo1 (a r sign)
1247 (cond ((eq sign '$positive)
1248 (throw 'isumout 'divergent))
1249 ((eq sign '$zero)
1250 (throw 'isumout 'divergent))
1251 ((eq sign '$negative)
1252 (adsum (list '(mtimes) a
1253 (list '(mexpt) (list '(mplus) 1 (list '(mtimes) -1 r)) -1))))))
1254
1255
1256 ;; Sums of polynomials using
1257 ;; bernpoly(x+1, n) - bernpoly(x, n) = n*x^(n-1)
1258 ;; which implies
1259 ;; sum(k^n, k, A, B) = 1/(n+1)*(bernpoly(B+1, n+1) - bernpoly(A, n+1))
1260 ;;
1261 ;; fpoly1 returns 1/(n+1)*(bernpoly(foo+1, n+1) - bernpoly(0, n+1)) for each power
1262 ;; in the polynomial e
1263
1264 (defun fpolysum (e lo hi) ;returns *ans*
1265 (let ((a (fpoly1 (setq e ($expand ($ratdisrep ($rat e *var*)))) lo))
1266 ($prederror))
1267 (cond ((null a) 0)
1268 ((member lo '(0 1))
1269 (maxima-substitute hi 'foo a))
1270 (t
1271 (list '(mplus) (maxima-substitute hi 'foo a)
1272 (list '(mtimes) -1 (maxima-substitute (list '(mplus) lo -1) 'foo a)))))))
1273
1274 (defun fpoly1 (e lo)
1275 (cond ((smono e *var*)
1276 (fpoly2 *a *n e lo))
1277 ((eq (caar e) 'mplus)
1278 (cons '(mplus) (mapcar #'(lambda (x) (fpoly1 x lo)) (cdr e))))
1279 (t (adusum e) 0)))
1280
1281 (defun fpoly2 (a n e lo)
1282 (cond ((null (and (integerp n) (> n -1))) (adusum e) 0)
1283 ((equal n 0)
1284 (m* (cond ((signp e lo)
1285 (m1+ 'foo))
1286 (t 'foo))
1287 a))
1288 (($ratsimp
1289 (m* a (list '(rat) 1 (1+ n))
1290 (m- ($bernpoly (m+ 'foo 1) (1+ n))
1291 ($bern (1+ n))))))))
1292
1293 ;; fbino can do these sums:
1294 ;; a) sum(binomial(n,k),k,0,n) -> 2^n
1295 ;; b) sum(binomial(n-k,k,k,0,n) -> fib(n+1)
1296 ;; c) sum(binomial(n,2k),k,0,n) -> 2^(n-1)
1297 ;; d) sum(binomial(a+k,b),k,l,h) -> binomial(h+a+1,b+1) - binomial(l+a,b+1)
1298 (defun fbino (e y lo hi)
1299 ;; e=binomial(n,d)
1300 (prog (n d l h)
1301 ;; check that n and d are linear in *var*
1302 (when (null (setq n (m2 (cadr e) (list 'n 'linear* *var*))))
1303 (return (adusum e)))
1304 (setq n (cdr (assoc 'n n :test #'eq)))
1305 (when (null (setq d (m2 (caddr e) (list 'd 'linear* *var*))))
1306 (return (adusum e)))
1307 (setq d (cdr (assoc 'd d :test #'eq)))
1308
1309 ;; binomial(a+b*k,c+b*k) -> binomial(a+b*k, a-c)
1310 (when (equal (cdr n) (cdr d))
1311 (setq d (cons (m- (car n) (car d)) 0)))
1312
1313 (cond
1314 ;; substitute k with -k in sum(binomial(a+b*k, c-d*k))
1315 ;; and sum(binomial(a-b*k,c))
1316 ((and (numberp (cdr d))
1317 (or (minusp (cdr d))
1318 (and (zerop (cdr d))
1319 (numberp (cdr n))
1320 (minusp (cdr n)))))
1321 (rplacd d (- (cdr d)))
1322 (rplacd n (- (cdr n)))
1323 (setq l (m- hi)
1324 h (m- lo)))
1325 (t (setq l lo h hi)))
1326
1327 (cond
1328
1329 ;; sum(binomial(a+k,c),k,l,h)
1330 ((and (equal 0 (cdr d)) (equal 1 (cdr n)))
1331 (adsum (m* y (m- (list '(%binomial) (m+ h (car n) 1) (m+ (car d) 1))
1332 (list '(%binomial) (m+ l (car n)) (m+ (car d) 1))))))
1333
1334 ;; sum(binomial(n,k),k,0,n)=2^n
1335 ((and (equal 1 (cdr d)) (equal 0 (cdr n)))
1336 ;; sum(binomial(n,k+c),k,l,h)=sum(binomial(n,k+c+l),k,0,h-l)
1337 (let ((h1 (m- h l))
1338 (c (m+ (car d) l)))
1339 (if (and (integerp (m- (car n) h1))
1340 (integerp c))
1341 (progn
1342 (adsum (m* y (m^ 2 (car n))))
1343 (when (member (asksign (m- (m+ h1 c) (car n))) '($zero $negative) :test #'eq)
1344 (adsum (m* -1 y (dosum (list '(%binomial) (car n) *var*)
1345 *var* (m+ h1 c 1) (car n) t :evaluate-summand nil))))
1346 (when (> c 0)
1347 (adsum (m* -1 y (dosum (list '(%binomial) (car n) *var*)
1348 *var* 0 (m- c 1) t :evaluate-summand nil)))))
1349 (adusum e))))
1350
1351 ;; sum(binomial(b-k,k),k,0,floor(b/2))=fib(b+1)
1352 ((and (equal -1 (cdr n)) (equal 1 (cdr d)))
1353 ;; sum(binomial(a-k,b+k),k,l,h)=sum(binomial(a+b-k,k),k,l+b,h+b)
1354 (let ((h1 (m+ h (car d)))
1355 (l1 (m+ l (car d)))
1356 (a1 (m+ (car n) (car d))))
1357 ;; sum(binomial(a1-k,k),k,0,floor(a1/2))=fib(a1+1)
1358 ;; we only do sums with h>floor(a1/2)
1359 (if (and (integerp l1)
1360 (member (asksign (m- h1 (m// a1 2))) '($zero $positive) :test #'eq))
1361 (progn
1362 (adsum (m* y ($fib (m+ a1 1))))
1363 (when (> l1 0)
1364 (adsum (m* -1 y (dosum (list '(%binomial) (m- a1 *var*) *var*)
1365 *var* 0 (m- l1 1) t :evaluate-summand nil)))))
1366 (adusum e))))
1367
1368 ;; sum(binomial(n,2*k),k,0,floor(n/2))=2^(n-1)
1369 ;; sum(binomial(n,2*k+1),k,0,floor((n-1)/2))=2^(n-1)
1370 ((and (equal 0 (cdr n)) (equal 2 (cdr d)))
1371 ;; sum(binomial(a,2*k+b),k,l,h)=sum(binomial(a,2*k),k,l+b/2,h+b/2), b even
1372 ;; 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
1373 (let ((a (car n))
1374 (r1 (if (oddp (car d)) 1 0))
1375 (l1 (if (oddp (car d))
1376 (m+ l (truncate (1- (car d)) 2))
1377 (m+ l (truncate (car d) 2)))))
1378 (when (and (integerp l1)
1379 (member (asksign (m- a hi)) '($zero $positive) :test #'eq))
1380 (adsum (m* y (m^ 2 (m- a 1))))
1381 (when (> l1 0)
1382 (adsum (m* -1 y (dosum (list '(%binomial) a (m+ *var* *var* r1))
1383 *var* 0 (m- l1 1) t :evaluate-summand nil)))))))
1384
1385 ;; other sums we can't do
1386 (t
1387 (adusum e)))))
1388
1389 ;; product routines
1390
1391 (defmspec $product (l)
1392 (arg-count-check 4 l)
1393 (setq l (cdr l))
1394 (dosum (car l) (cadr l) (meval (caddr l)) (meval (cadddr l)) nil :evaluate-summand t))
1395
1396 (declare-top (special $ratsimpexpons))
1397
1398 ;; Is this guy actually looking at the value of its middle arg?
1399
1400 (defun simpprod (x y z)
1401 (let (($ratsimpexpons t))
1402 (cond ((equal y 1)
1403 (setq y (simplifya (cadr x) z)))
1404 ((setq y (simptimes (list '(mexpt) (cadr x) y) 1 z)))))
1405 (simpprod1 y (caddr x)
1406 (simplifya (cadddr x) z)
1407 (simplifya (cadr (cdddr x)) z)))
1408
1409 (defmfun $taytorat (e)
1410 (cond ((mbagp e) (cons (car e) (mapcar #'$taytorat (cdr e))))
1411 ((or (atom e) (not (member 'trunc (cdar e) :test #'eq))) (ratf e))
1412 ((catch 'srrat (srrat e)))
1413 (t (ratf ($ratdisrep e)))))
1414
1415 (defun srrat (e)
1416 (unless (some (lambda (v) (switch 'multivar v)) (mrat-tlist e))
1417 (cons (list 'mrat 'simp (mrat-varlist e) (mrat-genvar e))
1418 (srrat2 (mrat-ps e)))))
1419
1420 (defun srrat2 (e)
1421 (if (pscoefp e) e (srrat3 (terms e) (gvar e))))
1422
1423 (defun srrat3 (l *var*)
1424 (cond ((null l) '(0 . 1))
1425 ((null (=1 (cdr (le l))))
1426 (throw 'srrat nil))
1427 ((null (n-term l))
1428 (rattimes (cons (list *var* (car (le l)) 1) 1)
1429 (srrat2 (lc l))
1430 t))
1431 ((ratplus
1432 (rattimes (cons (list *var* (car (le l)) 1) 1)
1433 (srrat2 (lc l))
1434 t)
1435 (srrat3 (n-term l) *var*)))))
1436
1437
1438 (declare-top (special $props *i))
1439
1440 (defmspec $deftaylor (l)
1441 (prog (fun series param op ops)
1442 a (when (null (setq l (cdr l))) (return (cons '(mlist) ops)))
1443 (setq fun (meval (car l)) series (meval (cadr l)) l (cdr l) param () )
1444 (when (or (atom fun)
1445 (if (eq (caar fun) 'mqapply)
1446 (or (cdddr fun) ; must be one parameter
1447 (null (cddr fun)) ; must have exactly one
1448 (do ((subs (cdadr fun) (cdr subs)))
1449 ((null subs)
1450 (setq op (caaadr fun))
1451 (when (cddr fun)
1452 (setq param (caddr fun)))
1453 '())
1454 (unless (atom (car subs)) (return 't))))
1455 (progn
1456 (setq op (caar fun))
1457 (when (cdr fun) (setq param (cadr fun)))
1458 (or (and (zl-get op 'op) (not (eq op 'mfactorial)))
1459 (not (atom (cadr fun)))
1460 (not (= (length fun) 2))))))
1461 (merror (intl:gettext "deftaylor: don't know how to handle this function: ~M") fun))
1462 (when (zl-get op 'sp2)
1463 (mtell (intl:gettext "deftaylor: redefining ~:M.~%") op))
1464 (when param (setq series (subst 'sp2var param series)))
1465 (setq series (subsum '*index series))
1466 (putprop op series 'sp2)
1467 (when (eq (caar fun) 'mqapply)
1468 (putprop op (cdadr fun) 'sp2subs))
1469 (add2lnc op $props)
1470 (push op ops)
1471 (go a)))
1472
1473 (defun subsum (*i e) (susum1 e))
1474
1475 (defun susum1 (e)
1476 (cond ((atom e) e)
1477 ((eq (caar e) '%sum)
1478 (if (null (smonop (cadr e) 'sp2var))
1479 (merror (intl:gettext "deftaylor: argument must be a power series at 0."))
1480 (subst *i (caddr e) e)))
1481 (t (recur-apply #'susum1 e))))
1482
1483 (declare-top (special varlist genvar $factorflag $ratfac *p* *var* *x*))
1484
1485 (defmfun $polydecomp (e v)
1486 (let ((varlist (list v))
1487 (genvar nil)
1488 *var* p den $factorflag $ratfac)
1489 (setq p (cdr (ratf (ratdisrep e)))
1490 *var* (cdr (ratf v)))
1491 (cond ((or (null (cdr *var*))
1492 (null (equal (cdar *var*) '(1 1))))
1493 (merror (intl:gettext "polydecomp: second argument must be an atom; found ~M") v))
1494 (t (setq *var* (caar *var*))))
1495 (cond ((or (pcoefp (cdr p))
1496 (null (eq (cadr p) *var*)))
1497 (setq den (cdr p)
1498 p (car p)))
1499 (t (merror (intl:gettext "polydecomp: cannot apply 'polydecomp' to a rational function."))))
1500 (cons '(mlist)
1501 (cond ((or (pcoefp p)
1502 (null (eq (car p) *var*)))
1503 (list (rdis (cons p den))))
1504 (t (setq p (pdecomp p *var*))
1505 (do ((l
1506 (setq p (mapcar #'(lambda (q) (cons q 1)) p))
1507 (cdr l))
1508 (a))
1509 ((null l)
1510 (cons (rdis (cons (caar p)
1511 (ptimes (cdar p) den)))
1512 (mapcar #'rdis (cdr p))))
1513 (cond ((setq a (pdecpow (car l) *var*))
1514 (rplaca l (car a))
1515 (cond ((cdr l)
1516 (rplacd l
1517 (cons (ratplus
1518 (rattimes
1519 (cadr l)
1520 (cons (ptterm (cdaadr a) 1)
1521 (cdadr a))
1522 t)
1523 (cons
1524 (ptterm (cdaadr a) 0)
1525 (cdadr a)))
1526 (cddr l))))
1527 ((equal (cadr a)
1528 (cons (list *var* 1 1) 1)))
1529 (t (rplacd l (list (cadr a)))))))))))))
1530
1531
1532 ;;; POLYDECOMP is like $POLYDECOMP except it takes a poly in *POLY* format (as
1533 ;;; defined in SOLVE) (numerator of a RAT form) and returns a list of
1534 ;;; headerless rat forms. In otherwords, it is $POLYDECOMP minus type checking
1535 ;;; and conversions to/from general representation which SOLVE doesn't
1536 ;;; want/need on a general basis.
1537 ;;; It is used in the SOLVE package and as such it should have an autoload
1538 ;;; property
1539
1540 (defun polydecomp (p *var*)
1541 (let ($factorflag $ratfac)
1542 (cond ((or (pcoefp p)
1543 (null (eq (car p) *var*)))
1544 (cons p nil))
1545 (t (setq p (pdecomp p *var*))
1546 (do ((l (setq p (mapcar #'(lambda (q) (cons q 1)) p))
1547 (cdr l))
1548 (a))
1549 ((null l)
1550 (cons (cons (caar p)
1551 (cdar p))
1552 (cdr p)))
1553 (cond ((setq a (pdecpow (car l) *var*))
1554 (rplaca l (car a))
1555 (cond ((cdr l)
1556 (rplacd l
1557 (cons (ratplus
1558 (rattimes
1559 (cadr l)
1560 (cons (ptterm (cdaadr a) 1)
1561 (cdadr a))
1562 t)
1563 (cons
1564 (ptterm (cdaadr a) 0)
1565 (cdadr a)))
1566 (cddr l))))
1567 ((equal (cadr a)
1568 (cons (list *var* 1 1) 1)))
1569 (t (rplacd l (list (cadr a))))))))))))
1570
1571
1572
1573 (defun pdecred (f h *var*) ;f = g(h(*var*))
1574 (cond ((or (pcoefp h) (null (eq (car h) *var*))
1575 (equal (cadr h) 1)
1576 (null (zerop (rem (cadr f) (cadr h))))
1577 (and (null (pzerop (caadr (setq f (pdivide f h)))))
1578 (equal (cdadr f) 1)))
1579 nil)
1580 (t (do ((q (pdivide (caar f) h) (pdivide (caar q) h))
1581 (i 1 (1+ i))
1582 (*ans*))
1583 ((pzerop (caar q))
1584 (cond ((and (equal (cdadr q) 1)
1585 (or (pcoefp (caadr q))
1586 (null (eq (caar (cadr q)) *var*))))
1587 (psimp *var* (cons i (cons (caadr q) *ans*))))))
1588 (cond ((and (equal (cdadr q) 1)
1589 (or (pcoefp (caadr q))
1590 (null (eq (caar (cadr q)) *var*))))
1591 (and (null (pzerop (caadr q)))
1592 (setq *ans* (cons i (cons (caadr q) *ans*)))))
1593 (t (return nil)))))))
1594
1595 (defun pdecomp (p *var*)
1596 (let ((c (ptterm (cdr p) 0))
1597 (a) (*x* (list *var* 1 1)))
1598 (cons (pcplus c (car (setq a (pdecomp* (pdifference p c)))))
1599 (cdr a))))
1600
1601 (defun pdecomp* (*p*)
1602 (let ((a)
1603 (l (pdecgdfrm (pfactor (pquotient *p* *x*)))))
1604 (cond ((or (pdecprimep (cadr *p*))
1605 (null (setq a (pdecomp1 *x* l))))
1606 (list *p*))
1607 (t (append (pdecomp* (car a)) (cdr a))))))
1608
1609 (defun pdecomp1 (prod l)
1610 (cond ((null l)
1611 (and (null (equal (cadr prod) (cadr *p*)))
1612 (setq l (pdecred *p* prod *var*))
1613 (list l prod)))
1614 ((pdecomp1 prod (cdr l)))
1615 (t (pdecomp1 (ptimes (car l) prod) (cdr l)))))
1616
1617 (defun pdecgdfrm (l) ;Get list of divisors
1618 (do ((l (copy-list l ))
1619 (ll (list (car l))
1620 (cons (car l) ll)))
1621 (nil)
1622 (rplaca (cdr l) (1- (cadr l)))
1623 (cond ((signp e (cadr l))
1624 (setq l (cddr l))))
1625 (cond ((null l) (return ll)))))
1626
1627 (defun pdecprimep (x)
1628 (setq x (cfactorw x))
1629 (and (null (cddr x)) (equal (cadr x) 1)))
1630
1631 (defun pdecpow (p *var*)
1632 (setq p (car p))
1633 (let ((p1 (pderivative p *var*))
1634 p2 p1p p1c a lin p2p)
1635 (setq p1p (oldcontent p1)
1636 p1c (car p1p) p1p (cadr p1p))
1637 (setq p2 (pderivative p1 *var*))
1638 (setq p2p (cadr (oldcontent p2)))
1639 (and (setq lin (testdivide p1p p2p))
1640 (null (pcoefp lin))
1641 (eq (car lin) *var*)
1642 (list (ratplus
1643 (rattimes (cons (list *var* (cadr p) 1) 1)
1644 (setq a (ratreduce p1c
1645 (ptimes (cadr p)
1646 (caddr lin))))
1647 t)
1648 (ratdif (cons p 1)
1649 (rattimes a (cons (pexpt lin (cadr p)) 1)
1650 t)))
1651 (cons lin 1)))))
1652
1653 (declare-top (unspecial *mfactl *factlist donel nn* dn*
1654 *var* *ans* *n *a*
1655 *infsumsimp *times *plus sum usum makef))
Maxima, a Computer Algebra System