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