home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / combin.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  42.5 KB  |  1,403 lines
  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. (macsyma-module combin)
  13.  
  14.  
  15. (declare-top (special *mfactl *factlist donel nn* dn* splist *ans* *var* dict ans var
  16.           a* $zerobern *a *n $cflength *a* $prevfib hi lo
  17.           *infsumsimp *times *plus sum usum makef
  18.           varlist genvar $sumsplitfact gensim $ratfac $simpsum
  19.           $prederror $listarith $sumhack $prodhack
  20.           $ratprint $zeta%pi $bftorat)
  21.      (*lexpr $ratcoef $divide)
  22.      (fixnum %n %a* %i %m)
  23.      (genprefix comb))
  24.  
  25. (load-macsyma-macros mhayat rzmac ratmac)
  26.  
  27.  
  28. (declare-top (splitfile minfct))
  29.  
  30. (comment minfactorial and factcomb stuff)
  31.  
  32. (DEFMFUN $makefact (e)
  33.   (let ((makef t)) (if (atom e) e (simplify (makefact1 e)))))
  34.  
  35. (defun makefact1 (e)
  36.     (cond ((atom e) e)
  37.       ((eq (caar e) '%binomial)
  38.        (subst (makefact1 (cadr e)) 'x
  39.           (subst (makefact1 (caddr e)) 'y
  40.              '((mtimes) ((mfactorial) x)
  41.                     ((mexpt) ((mfactorial) y) -1)
  42.                     ((mexpt) ((mfactorial) ((mplus) x ((mtimes) -1 y)))
  43.                          -1)))))
  44.       ((eq (caar e) '%gamma)
  45.        (list '(mfactorial) (list '(mplus) -1 (makefact1 (cadr e)))))
  46.       ((eq (caar e) '$beta)
  47.        (makefact1 (subst (cadr e) 'x
  48.                  (subst (caddr e) 'y
  49.                     '((mtimes) ((%gamma) x)
  50.                            ((%gamma) y)
  51.                            ((mexpt) ((%gamma) ((mplus) x y)) -1))))))
  52.       (t (recur-apply #'makefact1 e))))
  53.  
  54. (DEFMFUN $makegamma (e)
  55.     (if (atom e) e (simplify (makegamma1 ($makefact e)))))
  56.  
  57. (DEFMFUN $minfactorial (e)
  58.   (let (*mfactl *factlist)
  59.     (if (specrepp e) (setq e (specdisrep e)))
  60.     (getfact e)
  61.     ;(maplist #'evfac1 *factlist)        ;this is to save consing
  62.     ;If it was to save consing, why was MAPLIST used?
  63.     (mapl #'evfac1 *factlist)
  64.     (setq e (evfact e))))
  65.  
  66. (defun evfact (e)
  67.     (cond ((atom e) e)
  68.       ((eq (caar e) 'mfactorial)
  69.        (cdr (zl-ASSOC (cadr e) *factlist)))
  70.       ((memq  (caar e) '(%sum %derivative %integrate %product))
  71.        (cons (list (caar e)) (cons (evfact (cadr e)) (cddr e))))
  72.       (t (recur-apply #'evfact e))))
  73.  
  74. (defun adfactl (e l)
  75.   (let (n)
  76.     (cond ((null l) (setq *mfactl (cons (list e) *mfactl)))
  77.       ((numberp (setq n ($ratsimp `((mplus) ,e ((mtimes) -1 ,(caar l))))))
  78.          (cond ((plusp n)
  79.             ;(signp g n)
  80.             (rplacd (car l) (cons e (cdar l))))
  81.            ((rplaca l (cons e (car l))))))
  82.       ((adfactl e (cdr l))))))
  83.  
  84. (defun getfact (e)
  85.     (cond ((atom e) nil)
  86.       ((eq (caar e) 'mfactorial)
  87.        (and (null (zl-MEMBER (cadr e) *factlist))
  88.         (prog2 (setq *factlist (cons (cadr e) *factlist))
  89.                (adfactl (cadr e) *mfactl))))
  90.       ((memq (caar e) '(%sum %derivative %integrate %product))
  91.        (getfact (cadr e)))
  92.       ((mapc #'getfact (cdr e)))))
  93.  
  94. (defun evfac1 (e)
  95.   (do ((al *mfactl (cdr al)))
  96.       ((zl-MEMBER (car e) (car al))
  97.        (rplaca e
  98.            (cons (car e)
  99.              (list '(mtimes)
  100.                (gfact (car e)
  101.                   ($ratsimp (list '(mplus) (car e)
  102.                           (list '(mtimes) -1 (caar al)))) 1)
  103.                (list '(mfactorial) (caar al))))))))
  104.  
  105. (DEFMFUN $factcomb (e)
  106.        (let ((varlist varlist ) genvar $ratfac (ratrep (and (not (atom e)) (eq (caar e) 'mrat))))
  107.         (and ratrep (setq e (ratdisrep e)))
  108.         (setq e (factcomb e)
  109.             e (cond ((atom e) e)
  110.               (t (simplify (cons (list (caar e))
  111.                          (mapcar #'factcomb1 (cdr e)))))))
  112.         (or $sumsplitfact (setq e ($minfactorial e)))
  113.         (cond (ratrep (ratf e))
  114.           (t e))))
  115.  
  116. (defun factcomb1 (e)
  117.     (cond ((free e 'mfactorial) e)
  118.       ((memq (caar e) '(mplus mtimes mexpt))
  119.        (cons (list (caar e)) (mapcar #'factcomb1 (cdr e))))
  120.       (t (setq e (factcomb e))
  121.          (cond ((atom e) e)
  122.            (t (cons (list (caar e))
  123.                 (mapcar #'factcomb1 (cdr e))))))))
  124.  
  125. (defun factcomb (e)
  126.     (cond ((atom e) e)
  127.       ((free e 'mfactorial) e)
  128.       ((memq (caar e) '(mplus mtimes))
  129.        (factpluscomb (factcombplus e)))
  130.       ((eq (caar e) 'mexpt)
  131.        (simpexpt (list '(mexpt) (factcomb (cadr e))
  132.                (factcomb (caddr e)))
  133.              1 nil))
  134.       ((eq (caar e) 'mrat)
  135.        (factrat e))
  136.       (t (cons (car e) (mapcar #'factcomb (cdr e))))))
  137.  
  138. (defun factrat (e)
  139.   (let (nn* dn*)
  140.        (setq e (factqsnt ($ratdisrep (cons (car e) (cons (cadr e) 1)))
  141.              ($ratdisrep (cons (car e) (cons (cddr e) 1)))))
  142.        (numden e)
  143.        (div* (factpluscomb nn*) (factpluscomb dn*))))
  144.  
  145. (defun factqsnt (num den)
  146.    (if (equal num 0) 0
  147.       (let (nn* dn* (e (factpluscomb (div* den num))))
  148.      (numden e)
  149.      (factpluscomb (div* dn* nn*)))))
  150.  
  151. (defun factcombplus (e)
  152.   (let (nn* dn*)
  153.        (do ((l1 (nplus e) (cdr l1))
  154.         (l2))
  155.        ((null l1)
  156.         (simplus (cons '(mplus)
  157.                (mapcar #'(lambda (q) (factqsnt (car q) (cdr q))) l2))
  158.              1 nil))
  159.        (numden (car l1))
  160.        (do ((l3 l2 (cdr l3))
  161.         (l4))
  162.            ((null l3) (setq l2 (nconc l2 (list (cons nn* dn*)))))
  163.            (setq l4 (car l3))
  164.            (cond ((not (free ($gcd dn* (cdr l4)) 'mfactorial))
  165.               (numden (list '(mplus) (div* nn* dn*)
  166.                     (div* (car l4) (cdr l4))))
  167.               (setq l2 (delq l4 l2 1))))))))
  168.  
  169. (defun factpluscomb (e)
  170.   (prog (donel fact indl tt)
  171.     tag (setq e (factexpand e)
  172.           fact (getfactorial e))
  173.     (or fact (return e))
  174.     (setq indl (mapcar #'(lambda (q) (factplusdep q fact))
  175.                (nplus e))
  176.           tt (factpowerselect indl (nplus e) fact)
  177.           e (cond ((cdr tt)
  178.                (cons '(mplus) (mapcar #'(lambda (q) (factplus2 q fact))
  179.                           tt)))
  180.               (t (factplus2 (car tt) fact))))
  181.     (go tag)))
  182.  
  183. (defun nplus (e)
  184.   (cond ((eq (caar e) 'mplus) (cdr e))
  185.     (t (list e))))
  186.  
  187. (defun factexpand (e)
  188.        (cond ((atom e) e)
  189.          ((eq (caar e) 'mplus)
  190.           (simplus (cons '(mplus) (mapcar #'factexpand (cdr e)))
  191.                1 nil))
  192.          ((free e 'mfactorial) e)
  193.          (t ($expand e))))
  194.  
  195. (defun getfactorial (e)
  196.     (cond ((atom e) nil)
  197.       ((memq (caar e) '(mplus mtimes))
  198.        (do ((e (cdr e) (cdr e))
  199.         (a))
  200.            ((null e) nil)
  201.            (setq a (getfactorial (car e)))
  202.            (and a (return a))))
  203.       ((eq (caar e) 'mexpt)
  204.        (getfactorial (cadr e)))
  205.       ((eq (caar e) 'mfactorial)
  206.        (and (null (memalike (cadr e) donel))
  207.         (list '(mfactorial)
  208.               (car (setq donel (cons (cadr e) donel))))))))
  209.  
  210. (defun factplusdep (e fact)
  211.     (cond ((alike1 e fact) 1)
  212.       ((atom e) nil)
  213.       ((eq (caar e) 'mtimes)
  214.        (do ((l (cdr e) (cdr l))
  215.         (e) (out))
  216.            ((null l) nil)
  217.            (setq e (car l))
  218.            (and (setq out (factplusdep e fact))
  219.             (return out))))
  220.       ((eq (caar e) 'mexpt)
  221.        (let ((fto (factplusdep (cadr e) fact)))
  222.         (and fto (simptimes (list '(mtimes) fto
  223.                       (caddr e)) 1 t))))
  224.       ((eq (caar e) 'mplus)
  225.        (same (mapcar #'(lambda (q) (factplusdep q fact))
  226.              (cdr e))))))
  227.  
  228. (defun same (l)
  229.     (do ((ca (car l))
  230.      (cd (cdr l) (cdr cd))
  231.      (cad))
  232.     ((null cd) ca)
  233.     (setq cad (car cd))
  234.     (or (alike1 ca cad) (return nil))))
  235.  
  236. (defun factpowerselect (indl e fact)
  237.   (let (l fl)
  238.        (do ((i indl (cdr i))
  239.         (j e (cdr j))
  240.         (expt) (exp))
  241.        ((null i) l)
  242.        (setq expt (car i)
  243.          exp (cond (expt
  244.                 (setq exp ($divide (car j) `((MEXPT) ,fact ,expt)))
  245.                 ;; (car j) need not involve fact^expt since
  246.                 ;; fact^expt may be the gcd of the num and denom
  247.                 ;; of (car j) and $divide will cancel this out.
  248.                 (if (not (equal (cadr exp) 0)) (cadr exp)
  249.                    (setq expt ()) (caddr exp)))
  250.                (t (car j))))
  251.        (cond ((null l) (setq l (list (list expt exp))))
  252.          ((setq fl (assolike expt l))
  253.           (nconc fl (list exp)))
  254.          (t (nconc l (list (list expt exp))))))))
  255.  
  256. (defun factplus2 (l fact)
  257.   (let ((expt (car l)))
  258.        (cond (expt (factplus0 (cond ((cddr l) (rplaca l '(mplus)))
  259.                     (t (cadr l)))
  260.                   expt (cadr fact)))
  261.          (t (rplaca l '(mplus))))))
  262.  
  263. (defun factplus0 (r e fact)
  264.     (do ((i -1 (f1- i))
  265.      (fpn fact (list '(mplus) fact i))
  266.      (j -1) (exp) (rfpn) (div))
  267.     (nil)
  268.     (setq rfpn (simpexpt (list '(mexpt) fpn -1) 1 nil))
  269.     (setq div (dypheyed r (simpexpt (list '(mexpt) rfpn e) 1 nil)))
  270.     (cond ((or (null (or $sumsplitfact (equal (cadr div) 0)))
  271.            (equal (car div) 0))
  272.            (return (simplus (cons '(mplus) (mapcar
  273.          #'(lambda (q)
  274.            (setq j (f1+ j))
  275.            (list '(mtimes) q (list '(mexpt)
  276.             (list '(mfactorial) (list '(mplus) fpn j)) e)))
  277.           (factplus1 (cons r exp) e fpn)))
  278.                 1 nil)))
  279.           (t (setq r (car div))
  280.          (setq exp (cons (cadr div) exp))))))
  281.  
  282. (defun factplus1 (exp e fact)
  283.     (do ((l exp (cdr l))
  284.      (i 2 (f1+ i))
  285.      (fpn (list '(mplus) fact 1) (list '(mplus) fact i))
  286.      (div))
  287.     ((null l) exp)
  288.     (setq div (dypheyed (car l) (list '(mexpt) fpn e)))
  289.     (and (or $sumsplitfact (equal (cadr div) 0))
  290.          (null (equal (car div) 0))
  291.          (rplaca l (cadr div))
  292.          (rplacd l (cons (cond ((cadr l)
  293.              (simplus (list '(mplus) (car div) (cadr l))
  294.                   1 nil))
  295.                    (t
  296.                     (setq donel
  297.                       (cons (simplus fpn 1 nil) donel))
  298.                     (car div)))
  299.                  (cddr l))))))
  300.  
  301. (defun dypheyed (r f)
  302.   (let (r1 p1 p2)
  303.        (newvar r)
  304.        (setq r1 (ratf f)
  305.          p1 (pdegreevector (cadr r1))
  306.          p2 (pdegreevector (cddr r1)))
  307.        (do ((i p1 (cdr i))
  308.         (j p2 (cdr j))
  309.         (k (caddar r1) (cdr k)))
  310.        ((null k) (kansel r (cadr r1) (cddr r1)))
  311.        (cond ((> (car i) (car j))
  312.           (return (cdr ($divide r f (car k)))))))))
  313.  
  314. (defun kansel (r n d)
  315.   (let (r1 p1 p2)
  316.        (setq r1 (ratf r)
  317.          p1 (testdivide (cadr r1) n)
  318.          p2 (testdivide (cddr r1) d))
  319.        (cond ((and p1 p2) (cons (rdis (cons p1 p2)) '(0)))
  320.          (t (cons '0 (list r))))))
  321.  
  322. (declare-top (splitfile eulbrn)
  323. ;     (array* (notype *eu* 1 *bn* 1 *bd* 1)
  324.          )
  325.  
  326. (comment euler and bernoulli stuff)
  327.  
  328. #+cl
  329. (eval-when (compile eval load)
  330. (defvar *bn* (vector nil -1. 1. -1. 5. -691. 7. -3617. 43867. -174611. 854513.
  331.            -236364091.
  332.         8553103. -23749461029. 8615841276005. -7709321041217.
  333.         2577687858367.))
  334. (defvar *bd* (vector nil 30. 42. 30. 66. 2730. 6. 510. 798. 330. 138. 2730. 6. 870. 14322.
  335.         510. 6.))
  336. (defvar *eu* (vector -1. 5. -61. 1385. -50521. 2702765. -199360981. 19391512145.
  337.         -2404879675441. 370371188237525. -69348874393137901.))
  338. )
  339.  
  340. #-cl
  341. (progn 'compile
  342. (defvar *bn* (*array nil t 20.))
  343. (defvar *bd* (*array nil t 20.))
  344. (defvar *eu* (*array nil t 20.))
  345.  
  346. (setq i 0)
  347. (mapc #'(lambda (q)
  348.       (store (aref *bn* (setq i (f1+ i))) q))
  349.       '(-1. 1. -1. 5. -691. 7. -3617. 43867. -174611. 854513. -236364091.
  350.         8553103. -23749461029. 8615841276005. -7709321041217. 2577687858367.))
  351.  
  352.  
  353. (setq i 0)
  354. (mapc #'(lambda (a)
  355.        (store (aref *bd* (setq i (f1+ i))) a))
  356.       '(30. 42. 30. 66. 2730. 6. 510. 798. 330. 138. 2730. 6. 870. 14322.
  357.         510. 6.))
  358.  
  359. (setq i -1)
  360. (mapc #'(lambda (a)
  361.       (store (aref *eu* (setq i (f1+ i))) a))
  362.       '(-1. 5. -61. 1385. -50521. 2702765. -199360981. 19391512145.
  363.         -2404879675441. 370371188237525. -69348874393137901.))
  364. )
  365.  
  366.  
  367. (putprop '*eu* 11 'lim)
  368. (putprop 'bern 16 'lim)
  369.  
  370. (defmfun $euler (s)
  371.   (setq s
  372.     (let ((%n 0) $float)
  373.       (cond ((or (not (fixnump s)) (< s 0)) (list '($euler) s))
  374.         ((= (setq %n s) 0) 1)
  375.         ($zerobern
  376.          (cond ((oddp %n) 0)
  377.            ((null (> (// %n 2) (get '*eu* 'lim)))
  378.             (aref *eu* (f1- (// %n 2))))
  379.            ((eq $zerobern '%$/#&)
  380.             (euler %n))
  381.            ((*rearray '*eu* t (f1+ (// %n 2))) 
  382.             (euler %n))))
  383.         ((null (> %n (get '*eu* 'lim)))
  384.          (aref *eu* (f1- %n)))
  385.         ((*rearray '*eu* t (f1+ %n))
  386.          (euler (f* 2 %n))))))
  387.   (simplify s))
  388.  
  389. (defun euler (%a*)
  390.     (prog (nom %k e fl $zerobern *a*)
  391.         (setq nom 1 %k %a* fl nil e 0 $zerobern '%$/#& *a* (f1+ %a*))
  392.     a    (cond ((= %k 0)
  393.                (setq e (minus e))
  394.                ;(eval (list 'store (list '*eu* (f1- (// %a* 2))) e))
  395.                (store (aref *eu* (f1- (// %a* 2))) e)
  396.                (putprop '*eu* (// %a* 2) 'lim)
  397.                (return e)))
  398.         (setq nom (nxtbincoef (f1+ (- %a* %k)) nom) %k (f1- %k))
  399.         (cond ((setq fl (null fl)) (go a)))
  400.         (setq e (plus e (times nom ($euler %k))))
  401.         (go a)))
  402.  
  403. (defmfun simpeuler
  404.  (x vestigial z) 
  405.  vestigial ;ignored.
  406.  (oneargcheck x)
  407.  (let ((u (simpcheck (cadr x) z)))
  408.       (if (and (fixnump u) (not (< u 0)))
  409.       ($euler u)
  410.       (eqtest (list '($euler) u) x))))
  411.  
  412. (defmfun $bern (s)
  413.  (setq s 
  414.    (let ((%n 0) $float)
  415.     (cond ((or (not (fixnump s)) (< s 0)) (list '($bern) s))
  416.           ((= (setq %n s) 0) 1)
  417.           ((= %n 1) '((rat) -1 2))
  418.           ((= %n 2) '((rat) 1 6))
  419.           ($zerobern
  420.            (cond ((oddp %n) 0)
  421.              ((null (> (setq %n (f1- (// %n 2))) (get 'bern 'lim)))
  422.               (list '(rat) (aref *bn* %n) (aref *bd* %n)))
  423.              ((eq $zerobern '$/#&) (bern  (f* 2 (f1+ %n))))
  424.              (t (*rearray '*bn* t (setq %n (f1+ %n)))
  425.                 (*rearray '*bd* t %n) (bern  (f* 2 %n)))))
  426.           ((null (> %n (get 'bern 'lim)))
  427.            (list '(rat) (aref *bn* %n) (aref *bd* %n)))
  428.           (t (*rearray '*bn* t (f1+ %n)) (*rearray '*bd* t (f1+ %n))
  429.              (bern %n)))))
  430.  (simplify s))
  431.  
  432. (defun bern (%a*)
  433.     (prog (nom %k bb a b $zerobern l *a*)
  434.         (setq %k 0 l (f1- %a*) %a* (f1+ %a*) nom 1 $zerobern '$/#& a 1 b 1 *a* (f1+ %a*))
  435.     a    (cond ((= %k l)
  436.                (setq bb (*red a (times -1 b %a*)))
  437.                (putprop 'bern (setq %a* (f1- (// %a* 2))) 'lim)
  438.                ;if bb is a list, then it will always be
  439.                ;((RAT SIMP) number number).  ergo, evaluation
  440.                ; by store is no-op.
  441.                ;(eval (list 'store (list '*bn* %a*) (cadr bb)))
  442.                ;(eval (list 'store (list '*bd* %a*) (caddr bb)))
  443.                (store (aref *bn* %a*) (cadr bb))
  444.                (store (aref *bd* %a*) (caddr bb))
  445.                (return bb)))
  446.         (setq %k (f1+ %k))
  447.         (setq a (plus (times b (setq nom (nxtbincoef %k nom))
  448.                      (num1 (setq bb ($bern %k))))
  449.                   (times a (denom1 bb))))
  450.         (setq b (times b (denom1 bb)))
  451.         (setq a (*red a b) b (denom1 a) a (num1 a))
  452.         (go a)))
  453.  
  454. (defmfun simpbern (x vestigial z) 
  455.  vestigial ;ignored.
  456.  (oneargcheck x)
  457.  (let ((u (simpcheck (cadr x) z)))
  458.       (if (and (fixnump u) (not (< u 0)))
  459.       ($bern u)
  460.       (eqtest (list '($bern) u) x))))
  461.  
  462. (DEFMFUN $bernpoly (x s)
  463.   (let ((%n 0))
  464.        (cond ((not (fixnump s)) (list '($bernpoly) x s))
  465.          ((> (setq %n s) -1)
  466.           (do ((sum (cons (if (and (= %n 0) (zerop1 x))
  467.                    (addk 1 x)
  468.                    (power x %n)) nil)
  469.             (cons (m* (timesk (binocomp %n %k) ($bern %k))
  470.                   (if (and (= %n %k) (zerop1 x))
  471.                       (addk 1 x)
  472.                       (m^ x (f- %n %k))))
  473.                   sum))
  474.            (%k 1 (f1+ %k)))
  475.           ((> %k %n) (addn sum t))))
  476.          (t (list '($bernpoly) x %n)))))
  477.  
  478.  
  479. (declare-top (splitfile zeta))
  480.  
  481. (comment zeta and fibonacci stuff)
  482.  
  483. (DEFMFUN $zeta (s)
  484.  (cond ((null (fixnump s)) (list '($zeta) s))
  485.        ((oddp s)
  486.         (cond ((greaterp s 0)
  487.            (list '($zeta) s))
  488.           ((setq s (*dif 1 s))
  489.            (timesk -1
  490.             (timesk ($bern s) (list '(rat) (expt -1 (*quo s 2)) s))))))
  491.        ((equal s 0) '((rat simp) -1 2))
  492.        ((minusp s) 0)
  493.        ((not $zeta%pi) (list '($zeta) s))
  494.        (t (let ($numer $float)
  495.            (setq s (mul2 (power '$%pi s)
  496.                  (timesk (*red (expt 2 (sub1 s)) (factorial s))
  497.                      (simpabs (list 'mabs ($bern s)) 1 nil)))))
  498.       (resimplify s))))
  499.  
  500. (DEFMFUN $fib (n) 
  501.        (cond ((fixnump n) (ffib n))
  502.          (t (setq $prevfib (list '($fib) (add2* n -1)))
  503.         (list '($fib) n))))
  504.  
  505. (defun ffib (%n) 
  506.        (cond ((or (equal %n -1.) (zerop %n))
  507.           (setq $prevfib (boole boole-ior %n 1.) *a (f- %n)))
  508.          (t (let ((x (plus (ffib (// (boole  boole-andc2 %n 1.) 2.)) $prevfib)) (y (times $prevfib $prevfib)) (z (times *a *a))) (setq *a (*dif (times x x) y) 
  509.                                                                  $prevfib (plus y z)))
  510.          (cond ((oddp %n)
  511.             (setq *a (prog2 nil
  512.                     (plus *a $prevfib)
  513.                     (setq $prevfib *a))))
  514.                (*a))))) 
  515.  
  516. (declare-top (splitfile cffun))
  517.  
  518. (comment continued fraction stuff)
  519.  
  520. (DEFMFUN $cfdisrep (a)
  521.     (cond ((not ($listp a))
  522.            (merror "Arg to CFDISREP not a list: ~M" A))
  523.           ((null (cddr a)) (cadr a))
  524.           ((zerop (cadr a))
  525.            (list '(mexpt) (cfdisrep1 (cddr a)) -1))
  526.           ((cfdisrep1 (cdr a)))))
  527.  
  528. (defun cfdisrep1 (a)
  529.     (cond ((cdr a)
  530.            (list '(mplus simp cf) (car a)
  531.              (prog2 (setq a (cfdisrep1 (cdr a)))
  532.                 (cond ((integerp a) (list '(rat simp) 1 a))
  533.                   (t (list '(mexpt simp) a -1))))))
  534.           ((car a))))
  535.  
  536. (defun cfmak (a)
  537.     (setq a (meval a))
  538.     (cond ((integerp a) (list a))
  539.           ((eq (caar a) 'mlist) (cdr a))
  540.           ((eq (caar a) 'rat) (ratcf (cadr a) (caddr a)))
  541.           ((merror "Continued fractions must be lists or integers"))))
  542.  
  543. (defun makcf (a)
  544.     (cond ((null (cdr a)) (car a))
  545.           ((cons '(mlist simp cf) a))))
  546.  
  547. ;;; Translation properties for $CF defined in MAXSRC;TRANS5 >
  548.  
  549. (defmacro bind-status-divov-t (&rest forms)
  550.       (cond
  551.         #-cl
  552.         ((status feature maclisp)
  553.          `(let ((divov (status divov)))
  554.                (unwind-protect
  555.             (progn (sstatus divov t)
  556.                    ,@forms)
  557.             (sstatus divov divov))))
  558.         (t
  559.          `(progn ,@forms))))
  560.  
  561. (DEFMSPEC $cf (a)
  562.  (cfratsimp  (let ($listarith)
  563.        (bind-status-divov-t (cfeval   (meval (fexprcheck a)))))))
  564.  
  565. (defun cfratsimp (a)
  566.   (cond ((memq (car a) '(cf)) a)
  567.     (t
  568.      (cons '(mlist cf simp)(apply 'find-cf (cf-back-recurrence (cdr a)))))))
  569.  
  570. (defun cfeval (a)
  571.     (let (temp $ratprint)
  572.      (cond ((integerp a) (list '(mlist cf) a))
  573.            ((floatp a)
  574.         (let ((a (MAXIMA-RATIONALIZE a)))
  575.              (cons '(mlist cf) (ratcf (car a) (cdr a)))))
  576.            (($bfloatp a)
  577.         (let (($bftorat t))
  578.              (setq a (bigfloat2rat a))
  579.              (cons '(mlist cf) (ratcf (car a) (cdr a)))))
  580.            ((atom a)
  581.         (merror "~:M - not a continued fraction" a))
  582.            ((eq (caar a) 'rat)
  583.         (cons '(mlist cf) (ratcf (cadr a) (caddr a))))
  584.            ((eq (caar a) 'mlist)
  585.         (cfratsimp a))
  586.         ;;the following doesn't work for non standard form
  587. ;        (cfplus a '((mlist) 0)))
  588.            ((and (mtimesp a) (cddr a) (null (cdddr a))
  589.              (fixnump (cadr a))
  590.              (mexptp (caddr a))
  591.              (fixnump (cadr (caddr a)))
  592.              (alike1 (caddr (caddr a)) '((rat) 1 2)))
  593.         (cfsqrt (cfeval (times (expt (cadr a) 2) (cadr (caddr a))))))
  594.            ((eq (caar a) 'mexpt)
  595.         (cond ((alike1 (caddr a) '((rat) 1 2))
  596.                (cfsqrt (cfeval (cadr a))))
  597.               ((integerp (m- (caddr a) '((rat) 1 2)))
  598.                (cftimes (cfsqrt (cfeval (cadr a)))
  599.                 (cfexpt (cfeval (cadr a))
  600.                     (m- (caddr a) '((rat) 1 2)))))
  601.               ((cfexpt (cfeval (cadr a)) (caddr a)))))
  602.            ((setq temp (assq (caar a) '((mplus . cfplus) (mtimes . cftimes) (mquotient . cfquot)
  603.                                  (mdifference . cfdiff) (mminus . cfminus))))
  604.         (cf (cfeval (cadr a)) (cddr a) (cdr temp)))
  605.            ((eq (caar a) 'mrat)
  606.         (cfeval ($ratdisrep a)))
  607.            (t (merror "Not a continued fraction:~%~M" a)))))
  608.  
  609. (defun cf (a l fun)
  610.     (cond ((null l) a)
  611.           ((cf (FUNCALL fun a (Meval (list '($cf) (car l)))) (cdr l) fun))))
  612.  
  613. (defun cfplus (a b)
  614.     (setq a (cfmak a) b (cfmak b))
  615.     (makcf (cffun '(0 1 1 0) '(0 0 0 1) a b)))
  616.  
  617. (defun cftimes (a b)
  618.     (setq a (cfmak a) b (cfmak b))
  619.     (makcf (cffun '(1 0 0 0) '(0 0 0 1) a b)))
  620.  
  621. (defun cfdiff (a b)
  622.     (setq a (cfmak a) b (cfmak b))
  623.     (makcf (cffun '(0 1 -1 0) '(0 0 0 1) a b)))
  624.  
  625. (defun cfmin (a)
  626.     (setq a (cfmak a))
  627.     (makcf (cffun '(0 0 -1 0) '(0 0 0 1) a '(0))))
  628.  
  629. (defun cfquot (a b)
  630.     (setq a (cfmak a) b (cfmak b))
  631.     (makcf (cffun '(0 1 0 0) '(0 0 1 0) a b)))
  632.  
  633. (defun cfexpt (b e)
  634.        (setq b (cfmak b))
  635.        (cond ((null (integerp e))
  636.           (merror "Can't raise continued fraction to non-integral powers"))
  637.          ((let ((n (abs e)))
  638.             (do ((n (// n 2) (// n 2))
  639.              (s (cond ((oddp n) b)
  640.                   (t '(1)))))
  641.             ((zerop n)
  642.              (makcf
  643.               (cond ((signp g e)
  644.                  s)
  645.                    ((cffun '(0 0 0 1) '(0 1 0 0) b '(1))))))
  646.             (setq b (cffun '(1 0 0 0) '(0 0 0 1) b b))
  647.             (and (oddp n)
  648.                  (setq s (cffun '(1 0 0 0) '(0 0 0 1) s b))))))))
  649.  
  650.  
  651.  
  652. ;(defun conf1 (f g a b)
  653. ;  (*quo (conf2 f a b ) (conf2 g a b)))
  654.  
  655. (defun conf1 (f g a b &aux (den (conf2 g a b)))
  656.   (cond ((zerop den)
  657.      (* (signum (conf2 f a b ))  ;#+cl (/ most-positive-fixnum (^ 2 4))
  658.           #.(expt 2 31)))
  659.     (t (*quo (conf2 f a b) den))))
  660.  
  661. (defun conf2 (n a b) ;2*(abn_0+an_1+bn_2+n_3)
  662.     (times 2 (plus (times (car n) a b)
  663.                (times (cadr n) a)
  664.                (times (caddr n) b)
  665.                (cadddr n))))
  666.  
  667. ;;(cffun '(0 1 1 0) '(0 0 0 1) '(1 2) '(1 1 1 2)) gets error
  668. ;;should give (3 10)
  669.  
  670. (defun cf-convergents-p-q (cf &optional (n (length cf)) &aux pp qq)
  671.   "returns two lists such that pp_i/qq_i is the quotient of the first i terms
  672.    of cf"
  673.   (case (length cf)
  674.        (0 1)
  675.        (1  cf(list 1))
  676.        (t
  677.          (setq pp (list (add1
  678.                   (* (first cf) (second cf)))
  679.                 (car cf)))
  680.          (setq qq (list (second cf) 1))
  681.          (show pp qq)
  682.          (setq cf (cddr cf))
  683.          (sloop for i from 2 to n
  684.               while cf
  685.            do 
  686.               (push (+  (* (car cf) (car pp))
  687.                 (second pp)) pp)
  688.               (push (+  (* (car cf) (car qq))
  689.                 (second qq)) qq)
  690.                (setq cf (cdr cf))
  691.            finally (return (list (reverse pp) (reverse qq)))))))
  692.  
  693.  
  694. (defun find-cf1 (p q &optional so-far quot zl-REM )
  695.   (setq quot (quotient p q))
  696.   (setq zl-REM (- p (* quot q) ))
  697.   (prog ()
  698.     (cond ((< zl-REM 0) (incf zl-REM q) (incf quot -1))
  699.           ((zerop zl-REM) (return (cons quot so-far))))
  700.     (setq so-far (cons quot so-far))
  701.     (return (find-cf1 q zl-REM so-far))))
  702.  
  703. ;(compare-functions 
  704. ;;the following are equivalent but find-cf faster.
  705.  
  706. (defun find-cf (p q)
  707.   "returns the continued fraction for p and q integers, q not zero"
  708.   (cond  ((zerop q)(error "quotient by zero"))
  709.      ((< q 0) (setq p (f- p)) (setq q (f- q))))
  710.   (nreverse (find-cf1 p q))) 
  711.  
  712. ;this works but since it is slower than find-cf
  713. ;(defun standard-cf (p q &aux start  tem new-q)
  714. ; "returns the continued fraction for p and q integers, q not zero"
  715. ;  (setq start (floor p q))
  716. ;  (setq p (f- p (f* start q)))
  717. ;  (sloop until (zerop q)
  718. ;    when start
  719. ;      collecting start
  720. ;    and
  721. ;    do(setq tem (floor p q))(setq start nil)
  722. ;    else
  723. ;      collecting (setq tem (floor p q))
  724. ;    do
  725. ;           (setq new-q (f- p (f* tem q)))
  726. ;           (setq p q)
  727. ;           (setq q new-q)))
  728. ;)
  729.  
  730.  
  731. (defun cf-back-recurrence (cf &aux tem (num-gg 0)(den-gg 1))
  732.   "converts CF (a continued fraction list) to a list of numerator
  733.   denominator using  recurrence from end
  734.   and not calculating intermediate quotients.
  735.   The numerator and denom are relatively
  736.    prime"
  737.   (sloop for v in (reverse cf)
  738.     do (setq tem (* den-gg v))
  739.        (setq tem (+ tem num-gg))
  740.        (setq num-gg den-gg)
  741.        (setq den-gg tem)
  742.     finally
  743.       (return
  744.         (cond ((and (<= den-gg 0) (< num-gg 0))
  745.            (list  (- den-gg) (- num-gg)))
  746.           (t(list den-gg num-gg))))))
  747.  
  748. (declare-top(unspecial w))
  749. ;;(cffun '(0 1 1 0) '(0 0 0 1) '(1 2) '(1 1 1 2)) gets error
  750. ;;should give (3 10)
  751.  
  752. (defun cffun (f g a b)
  753.     (prog (c v w)
  754.       (declare (special v))
  755.     a   (and (zerop (cadddr g))
  756.          (zerop (caddr g))
  757.          (zerop (cadr g))
  758.          (zerop (car g))
  759.          (return (reverse c)))
  760.     (and (equal (setq w (conf1 f g (car a) (add1 (car b))))
  761.             (setq v (conf1 f g (car a) (car b))))
  762.          (equal (conf1 f g (add1 (car a)) (car b)) v)
  763.          (equal (conf1 f g (add1 (car a)) (add1 (car b))) v)
  764.          (setq g (mapcar #'(lambda (a b)
  765.                  (declare (special v))
  766.                  (*dif a (times v b)))
  767.                   f (setq f g)))
  768.          (setq c (cons v c))
  769. ;         (progn (show v w) t)
  770.          (go a))
  771. ;    (show v w)
  772.     (cond ((lessp (abs (*dif (conf1 f g (add1 (car a)) (car b)) v))
  773.               (abs (*dif w v)))
  774.            (cond ((setq v (cdr b))
  775.               (setq f (conf6 f b))
  776.               (setq g (conf6 g b))
  777.               (setq b v))
  778.              (t (setq f (conf7 f b)) (setq g  (conf7 g b)))))
  779.           (t
  780.            (cond ((setq v (cdr a))
  781.               (setq f (conf4 f a))
  782.               (setq g (conf4 g a))
  783.               (setq a v))
  784.              (t (setq f (conf5 f a)) (setq g  (conf5 g a))))))
  785.                (go a)))
  786.  
  787. (defun conf4 (n a) ;n_0*a_0+n_2,n_1*a_0+n_3,n_0,n_1
  788.     (list (plus (times (car n) (car a)) (caddr n))
  789.       (plus (times (cadr n) (car a)) (cadddr n))
  790.       (car n)
  791.       (cadr n)))
  792.  
  793. (defun conf5 (n a) ;0,0, n_0*a_0,n_2
  794.     (list 0 0
  795.       (plus (times (car n) (car a)) (caddr n))
  796.       (plus (times (cadr n) (car a)) (cadddr n))))
  797.  
  798. (defun conf6 (n b)
  799.     (list (plus (times (car n) (car b)) (cadr n))
  800.       (car n)
  801.       (plus (times (caddr n) (car b)) (cadddr n))
  802.       (caddr n)))
  803.  
  804. (defun conf7 (n b)
  805.     (list 0 (plus (times (car n) (car b)) (cadr n))
  806.       0 (plus (times (caddr n) (car b)) (cadddr n))))
  807.  
  808. (defun cfsqrt (n)
  809.     (cond ((cddr n)                    ;A non integer
  810.        (merror "Can't take square roots of non-integers yet"))
  811.       ((setq n (cadr n))))
  812.     (setq n (sqcont n))
  813.     (cond ((= $cflength 1)
  814.        (cons '(mlist simp) n))
  815.       ((do ((i 2 (f1+ i))
  816.         (a (copy (cdr n))))
  817.            ((> i $cflength) (cons '(mlist simp) n))
  818.            (setq n (nconc n (copy a)))))))    
  819.  
  820. (DEFMFUN $qunit (n)
  821.     (let ((l (sqcont n))) (list '(mplus) (pelso1 l 0 1) 
  822.                     (list '(mtimes) 
  823.                       (list '(mexpt) n '((rat) 1 2))
  824.                       (pelso1 l 1 0)))))
  825.  
  826. (defun pelso1 (l a b)
  827.     (do ((i l (cdr i))) (nil)
  828.         (and (null (cdr i)) (return b))
  829.         (setq b (plus a (times (car i) (setq a b))))))
  830.  
  831. (defun sqcont (n)
  832.     (prog (q q1 q2 m m1 a0 a l)
  833.         (setq a0 ($isqrt n) a (list a0) q2 1 m1 a0 
  834.                q1 (*dif n (times m1 m1)) l (times 2 a0))
  835.     a    (setq a (cons (*quo (plus m1 a0) q1) a))
  836.         (cond ((equal (car a) l)
  837.                (return (nreverse a))))
  838.         (setq m (*dif (times (car a) q1) m1)
  839.               q (plus q2 (times (car a) (*dif m1 m)))
  840.               q2 q1 q1 q m1 m)
  841.         (go a)))
  842.  
  843. (defun ratcf (x y)
  844.     (prog (a b)
  845.     a    (cond ((equal  y 1) (return (nreverse (cons x a))))
  846.               ((minusp x)
  847.                (setq b (plus y (remainder x y))
  848.                  a (cons (sub1 (*quo x y)) a)
  849.                  x y y b))
  850.               ((greaterp y x)
  851.                (setq a (cons 0 a))
  852.                (setq b x x y y b))
  853.               ((equal x y) (return (nreverse (cons 1 a))))
  854.               ((setq b (remainder x y))
  855.             (setq a (cons (*quo x y) a) x y y b)))
  856.         (go a)))
  857.  
  858. (DEFMFUN $cfexpand (x)
  859.     (cond ((null ($listp x)) x)
  860.           ((cons '($matrix) (cfexpand (cdr x))))))
  861.  
  862. (defun cfexpand (ll)
  863.     (do ((p1 0 p2)
  864.          (p2 1 (simplify (list '(mplus) (list '(mtimes) (car l) p2) p1)))
  865.          (q1 1 q2)
  866.          (q2 0 (simplify (list '(mplus) (list '(mtimes) (car l) q2) q1)))
  867.          (l ll (cdr l)))
  868.         ((null l) (list (list '(mlist) p2 p1) (list '(mlist) q2 q1)))))
  869.  
  870. (declare-top (splitfile sum)
  871.      (*expr sum))
  872.  
  873. (comment Summation stuff)
  874.  
  875. (defun adsum (e) (setq sum (cons (simplify e) sum)))
  876.  
  877. (defun adusum (e) (setq usum (cons (simplify e) usum)))
  878.  
  879. (defmfun simpsum2 (exp i lo hi)
  880.   (prog (*plus *times $simpsum u)
  881.     (setq *plus (list 0) *times 1)
  882.         (when (or (and (eq hi '$inf) (eq lo '$minf))
  883.           (equal 0 (m+ hi lo)))
  884.       (setq $simpsum t lo 0)
  885.       (setq *plus (cons (m* -1 *times (MAXIMA-SUBSTITUTE 0 i exp)) *plus))
  886.       (setq exp (m+ exp (MAXIMA-SUBSTITUTE (m- i) i exp))))
  887.     (cond ((and (null $sumhack)
  888.             (eq ($sign (setq u (m- hi lo))) '$neg))
  889.            (cond ((equal u -1) (return 0))
  890.              (t (merror "Lower bound to sum is > upper bound"))))
  891.           ((free exp i)
  892.            (return (m+l (cons (freesum exp lo hi *times) *plus))))
  893.  
  894.           ((setq exp (sumsum exp i lo hi))
  895.            (setq exp (m* *times (dosum (cadr exp) (caddr exp)
  896.                       (cadddr exp) (cadr (cdddr exp)) t))))
  897.           (t (return (m+l *plus))))
  898.        (return (m+l (cons exp *plus)))))
  899.  
  900. (defun sumsum (e *var* lo hi)
  901.    (let (sum usum)
  902.     (cond ((eq hi '$inf)
  903.            (cond (*infsumsimp (isum e))
  904.              ((setq usum (list e))))) 
  905.           ((sum e 1)))
  906.     (cond ((eq sum nil)
  907.            (return-from sumsum (list '(%sum) e *var* lo hi))))
  908.     (setq *plus
  909.           (nconc (mapcar 
  910.               #'(lambda (q) (simptimes (list '(mtimes) *times q) 1 nil))
  911.               sum)
  912.              *plus))
  913.     (and usum (setq usum (list '(%sum) (simplus (cons '(plus) usum) 1 t) *var* lo hi)))))
  914.  
  915. (defun sum (e y)
  916.   (cond ((null e))
  917.     ((free e *var*)
  918.      (adsum (m* y e (m+ hi 1 (m- lo)))))
  919.     ((poly? e *var*)
  920.      (adsum (m* y (fpolysum e))))
  921.     ((eq (caar e) '%binomial) (fbino e y))
  922.     ((eq (caar e) 'mplus)
  923.      (mapc #'(lambda (q) (sum q y)) (cdr e)))
  924.     ((and (or (mtimesp e) (mexptp e) (mplusp e))
  925.           (fsgeo e y)))
  926.     (t
  927.      nil
  928.      #-cl
  929.      (let (*a *n)
  930.          (cond ((prog2 (m2 e '((mtimes) ((coefftt) (var* (set) *a freevar))
  931.                         ((coefftt) (var* (set) *n true)))
  932.                    nil)
  933.                (not (equal *a 1)))
  934.             (sum *n (list '(mtimes) y *a)))
  935.            ((and (not (atom
  936.                    (setq *n
  937.                      ($ratdisrep
  938.                       (let (genvar (varlist (cons *var* nil)))
  939.                     (ratrep* *n))))))
  940.              (not (equal *n e))
  941.              (not (eq (caar *n) 'mtimes)))
  942.             (sum *n (list '(mtimes) y *a)))
  943.            (t (adusum (list '(mtimes) e y))))))))
  944.  
  945. (defun isum (e)
  946.     (cond ((memq (setq e (catch 'isumout (isum1 e))) '($inf $undefined $minf))
  947.            (setq sum (list e) usum nil))))
  948.  
  949. (defun isum1 (e)
  950.     (cond ((or (free e *var*) (atom e))
  951.            (throw 'isumout '$inf))
  952.           ((ratp e *var*)
  953.            (adsum (ipolysum e)))
  954.           ((eq (caar e) 'mplus)
  955.            (mapc #'isum1 (cdr e)))
  956.           ( (isgeo e))
  957.           ((adusum e))))
  958.  
  959. (defun ipolysum (e)
  960.     (ipoly1 ($expand e)))
  961.  
  962. (defun ipoly1 (e)
  963.     (cond ((smono e *var*)
  964.            (ipoly2 *a *n (asksign (simplify (list '(mplus) *n 1)))))
  965.           ((mplusp e)
  966.            (cons '(mplus) (mapcar #'ipoly1 (cdr e))))
  967.           (t (adusum e)
  968.          0)))
  969.  
  970. (defun ipoly2 (a n sign)
  971.     (cond ((memq (asksign lo) '($zero $negative))
  972.            (throw 'isumout '$inf)))
  973.        (and (null (equal lo 1))
  974.     (adsum `((%sum) 
  975.           ((mtimes) ,a -1 ((mexpt) ,*var* ,n))
  976.           ,*var* 1 ((mplus) -1 ,lo))))
  977.     (cond ((eq sign '$negative)
  978.            (list '(mtimes) a ($zeta (meval (list '(mtimes) -1 n)))))
  979.           ((throw 'isumout '$inf))))
  980.  
  981. (defun fsgeo (e y)
  982.    (let ((r ($ratsimp (div* (MAXIMA-SUBSTITUTE (list '(mplus) *var* 1) *var* e) e))))
  983.     (cond ((free r *var*)
  984.            (adsum 
  985.         (list '(mtimes) y
  986.               (MAXIMA-SUBSTITUTE 0 *var* e)
  987.               (list '(mplus)
  988.                 (list '(mexpt) r (list '(mplus) hi 1))
  989.                 (list '(mtimes) -1 (list '(mexpt) r lo)))
  990.               (list '(mexpt) (list '(mplus) r -1) -1)))))))
  991.  
  992. (defun isgeo (e)
  993.    (let ((r ($ratsimp (div* (MAXIMA-SUBSTITUTE (list '(mplus) *var* 1) *var* e) e))))
  994.     (and (free r *var*)
  995.          (isgeo1 (MAXIMA-SUBSTITUTE lo *var* e)
  996.              r (asksign (simplify (list '(mplus) (list '(mabs) r) -1)))))))
  997.  
  998. (defun isgeo1 (a r sign)
  999.     (cond ((eq sign '$positive)
  1000.            (cond ((mgrp a 0) (throw 'isumout '$inf))
  1001.              ((throw 'isumout '$minf))))
  1002.           ((eq sign '$zero)
  1003.            (throw 'isumout '$undefined))
  1004.           ((eq sign '$negative)
  1005.            (adsum (list '(mtimes) a
  1006.                  (list '(mexpt) (list '(mplus) 1 (list '(mtimes) -1 r)) -1))))))
  1007.  
  1008. (defun fpolysum (e)    ;returns *ans*
  1009.    (let ((a (fpoly1 (setq e ($expand ($ratdisrep ($rat e *var*))))))
  1010.       (b) ($prederror))
  1011.      (cond ((null a) 0)
  1012.            ((zl-MEMBER lo '(0 1))
  1013.         (MAXIMA-SUBSTITUTE hi 'foo a))
  1014.            ((or (equal t (mevalp (list '(mgeqp) lo 0)))
  1015.             (memq (asksign lo) '($zero $positive)))
  1016.         (list '(mplus) (MAXIMA-SUBSTITUTE hi 'foo a)
  1017.               (list '(mtimes) -1 (MAXIMA-SUBSTITUTE (list '(mplus) lo -1) 'foo a))))
  1018.           (t
  1019.            (setq b (fpoly1 (MAXIMA-SUBSTITUTE (list '(mtimes) -1 *var*) *var* e)))
  1020.            (list '(mplus) (MAXIMA-SUBSTITUTE hi 'foo a) (MAXIMA-SUBSTITUTE lo 'foo b))))))
  1021.  
  1022. (defun fpoly1 (e)
  1023.     (cond ((smono e *var*)
  1024.            (fpoly2 *a *n e))
  1025.           ((eq (caar e) 'mplus)
  1026.            (cons '(mplus) (mapcar #'fpoly1 (cdr e))))
  1027.           (t (adusum e) 0)))
  1028.  
  1029. (defun fpoly2 (a n e)
  1030.     (cond ((null (and (integerp n) (greaterp n -1))) (adusum e) 0)
  1031.           ((equal n 0)
  1032.            (m* (cond ((signp e lo)
  1033.               (m1+ 'foo))
  1034.              (t 'foo))
  1035.            a))
  1036.           (($ratsimp
  1037.         (m* a (list '(rat) 1 (f1+ n))
  1038.         (m- ($bernpoly (m+ 'foo 1) (f1+ n))
  1039.             ($bern (f1+ n))))))))
  1040.  
  1041. (defun fbino (e y)
  1042.     (prog (n d l h fl)
  1043.         (cond ((null (setq n (m2 (cadr e) (list 'n 'linear* *var*) nil)))
  1044.                (return (adusum e))))
  1045.         (setq n (cdr (assq 'n n)))
  1046.         (cond ((null (setq d (m2 (caddr e) (list 'd 'linear* *var*) nil)))
  1047.                (return (adusum e))))
  1048.         (setq d (cdr (assq 'd d)))
  1049.         (cond ((equal (cdr n) (cdr d))
  1050.                (setq d (cons (simplus (list '(mplus) (car n) 
  1051.                        (list '(mtimes) -1 (car d))) 1 nil) 0))))
  1052.                 (cond ((and (numberp (cdr d)) (or (minusp (cdr d)) (and (zerop (cdr d))
  1053.                         (numberp (cdr n)) (minusp (cdr n)))))
  1054.                        (rplacd d (minus (cdr d)))
  1055.                        (rplacd n (minus (cdr n)))
  1056.                        (setq l (simptimes (list '(mtimes) hi -1) 1 nil)
  1057.                              h (simptimes (list '(mtimes) lo -1) 1 nil)))
  1058.                       (t (setq l lo  h hi)))
  1059.         (cond ((null (or (zl-MEMBER (cdr n) '(0 -1)) (equal 0 (cdr d))))
  1060.                (return (adusum e)))
  1061.                       ((and (equal 0 (cdr d)) (equal 1 (cdr n)))
  1062.                        (return (adsum (list '(mplus) (list '(%binomial)
  1063.                                        (list '(mplus) h (car n) 1) (list '(mplus) (car d) 1))
  1064.                                        (list '(mtimes) -1 (list '(%binomial)
  1065.                                         (list '(mplus) l (car n)) (list '(mplus) (car d) 1)))))))
  1066.               ((equal 1 (cdr d))
  1067.                (cond ((equal -1 (cdr n))
  1068.                   (setq fl 0))))
  1069.               ((and (equal 2 (cdr d)) (equal 0 (cdr n)))
  1070.                (setq fl 1))
  1071.               ((return (adusum e))))
  1072.         (setq n (car n))
  1073.         (cond ((equal (cdr d) -1)
  1074.                (setq d (simplus (list '(mplus) n (list '(mtimes) -1 (car d))) 1 nil)))
  1075.               ((setq d (car d))))
  1076.         (cond ((equal fl 1) (setq d (*quo d 2))))
  1077.         (setq l (simplus (list '(mplus) l d) 1 nil)
  1078.               h (simplus (list '(mplus) h d) 1 nil))
  1079.         (cond ((or (null (numberp l)) 
  1080.                (null
  1081.                 (or (numberp (setq d (simplus (list '(mplus) h (list '(mtimes) -1 n)) 1 nil)))
  1082.                 (and (numberp (simplus (list '(mplus) n (list '(mtimes) -2 h)) 1 nil))
  1083.                      (progn (setq fl 2) t)))))
  1084.                (return (adusum e))))
  1085.         (setq e (list '(%binomial) n *var*))
  1086.         (and (greaterp l 0) (adsum (dosum (list '(mtimes) y -1 e) *var* 0 (sub1 l) t)))
  1087.         (and (lessp d 0)
  1088.              (adsum (dosum (list '(mtimes) y -1 e) *var* (simplus (list '(mplus) h 1) 1 nil) n t)))
  1089.         (and (greaterp d 0)
  1090.              (adsum (dosum (list '(mtimes) y e) *var* (simplus (list '(mplus) n 1) 1 nil) h t)))
  1091.         (setq fl
  1092.          (cond ((null fl) (list '(mexpt) 2 n))
  1093.                ((zerop fl) ($fib (simplus (list '(mplus) n 1) 1 nil)))
  1094.                ((list '(mexpt) 2 (list '(mplus) n -1)))))
  1095.         (adsum (list '(mtimes) y fl))))
  1096.  
  1097. (declare-top (splitfile prodct))
  1098.  
  1099. (comment product routines)
  1100.  
  1101. (DEFMSPEC $product (l) (SETQ l (CDR l))
  1102.     (cond ((not (= (length l) 4)) (merror "Wrong no. of args to product"))
  1103.       ((dosum (car l) (cadr l)   (meval (caddr l)) (meval (cadddr l)) nil))))
  1104.  
  1105. (declare-top(special $ratsimpexpons))
  1106.  
  1107. ;; Is this guy actually looking at the value of its middle arg?
  1108.  
  1109. (defun simpprod (x y z)
  1110.    (let (($ratsimpexpons t))
  1111.     (cond ((equal y 1)
  1112.            (setq y (simplifya (cadr x) z)))
  1113.           ((setq y (simptimes (list '(mexpt) (cadr x) y) 1 z)))))
  1114.    (simpprod1 y (caddr x)
  1115.            (simplifya (cadddr x) z)
  1116.            (simplifya (cadr (cdddr x)) z)))
  1117.  
  1118. (defun simpprod1 (exp i lo hi)
  1119.   (let (u)
  1120.        (cond ((not (symbolp i)) (merror "Bad index to product:~%~M" i))
  1121.          ((alike1 lo hi)
  1122.           (let ((valist (list i)))
  1123.            (mbinding (valist (list hi))
  1124.                  (meval exp))))
  1125.          ((and (null $prodhack)
  1126.            (eq ($sign (setq u (m- hi lo))) '$neg))
  1127.           (cond ((eq ($sign (add2 u 1)) '$zero) 1)
  1128.             (t (merror "Lower bound to product is > upper bound."))))
  1129.          ((atom exp)
  1130.           (cond ((null (eq exp i))
  1131.              (power* exp (list '(mplus) hi 1 (list '(mtimes) -1 lo))))
  1132.             ((let ((lot (asksign lo)))
  1133.               (cond ((equal lot '$zero) 0)
  1134.                 ((eq lot '$positive)
  1135.                  (m// (list '(mfactorial) hi)
  1136.                       (list '(mfactorial) (list '(mplus) lo -1))))
  1137.                 ((m* (list '(mfactorial)
  1138.                        (list '(mabs) lo))
  1139.                      (cond ((memq (asksign hi) '($zero $positive))
  1140.                         0)
  1141.                        (t (prog2 0
  1142.                              (m^ -1 (m+ hi lo 1))
  1143.                              (setq hi (list '(mabs) hi)))))
  1144.                      (list '(mfactorial) hi))))))))
  1145.          ((list '(%product simp) exp i lo hi)))))
  1146.  
  1147. (declare-top (splitfile tayrat))
  1148.  
  1149. (defmfun $taytorat (e)
  1150.   (cond ((mbagp e) (cons (car e) (mapcar #'$taytorat (cdr e))))
  1151.     ((or (atom e) (not (memq 'trunc (cdar e)))) (ratf e))
  1152.     ((catch 'srrat (srrat e))) 
  1153.     (t (ratf ($ratdisrep e)))))
  1154.  
  1155. (defun srrat (e)
  1156.        (cons (list 'mrat 'simp (caddar e) (cadddr (car e)))
  1157.          (srrat2 (cdr e))))
  1158.  
  1159. (defun srrat2 (e)
  1160.        (if (pscoefp e) e (srrat3 (terms e) (gvar e))))
  1161.  
  1162. (defun srrat3 (l *var*)
  1163.     (cond ((null l) '(0 . 1))
  1164.           ((null (=1 (cdr (le l))))
  1165.            (throw 'srrat nil))
  1166.           ((null (n-term l))
  1167.            (rattimes (cons (list *var* (car (le l)) 1) 1)
  1168.                  (srrat2 (lc l))
  1169.              t))
  1170.           ((ratplus
  1171.         (rattimes (cons (list *var* (car (le l)) 1) 1)
  1172.               (srrat2 (lc l))
  1173.               t)
  1174.         (srrat3 (n-term l) *var*)))))
  1175.  
  1176.  
  1177. (declare-top (splitfile deftay) (special $props *i))
  1178.  
  1179. (defmspec $deftaylor (l)
  1180.   (prog (fun series param op ops)
  1181.    a    (when (null (setq l (cdr l))) (return (cons '(MLIST) ops)))
  1182.     (setq fun (meval (car l)) series (meval (cadr l)) l (cdr l) param () )
  1183.     (when (or (atom fun) 
  1184.           (if (eq (caar fun) 'mqapply)
  1185.               (or (cdddr fun)    ; must be one parameter
  1186.               (null (cddr fun))    ; must have exactly one
  1187.               (do ((subs (cdadr fun) (cdr subs)))
  1188.                   ((null subs)
  1189.                    (setq op (caaadr fun))
  1190.                    (when (cddr fun) (setq param (caddr fun)))
  1191.                    () )
  1192.                  (unless (atom (car subs)) (return 'T))))
  1193.              (setq op (caar fun))
  1194.              (when (cdr fun) (setq param (cadr fun)))
  1195.              (or (and (oldget op 'op) (not (eq op 'mfactorial)))
  1196.              (not (atom (cadr fun)))
  1197.              (not (= (length fun) 2)))))
  1198.        (merror "Bad argument to DEFTAYLOR:~%~M" fun))
  1199.     (when (oldget op 'sp2)
  1200.        (mtell "~:M being redefined in DEFTAYLOR.~%" op))
  1201.     (when param (setq series (subst 'sp2var param series)))
  1202.     (setq series (subsum '*index series))
  1203.     (putprop op series 'sp2)
  1204.     (when (eq (caar fun) 'mqapply)
  1205.        (putprop op (cdadr fun) 'sp2subs))
  1206.     (add2lnc op $props)
  1207.     (push op ops)
  1208.     (go a)))
  1209.  
  1210. (defun subsum (*i e) (susum1 e))
  1211.  
  1212. (defun susum1 (e)
  1213.     (cond ((atom e) e)
  1214.           ((eq (caar e) '%sum)
  1215.            (if (null (smonop (cadr e) 'sp2var))
  1216.            (merror "Argument to DEFTAYLOR must be power series at 0.")
  1217.            (subst *i (caddr e) e)))
  1218.           (t (recur-apply #'susum1 e))))
  1219.  
  1220. (declare-top (splitfile decomp)
  1221.      (special varlist genvar $factorflag $ratfac *p* *var* *l* *x*))
  1222.  
  1223. (DEFMFUN $polydecomp (e v)
  1224.        (let ((varlist (list v))
  1225.           (genvar nil)
  1226.           *var* p den $FACTORFLAG $RATFAC)
  1227.          (setq p (cdr (ratf (ratdisrep e)))
  1228.            *var* (cdr (ratf v)))
  1229.          (cond ((or (null (cdr *var*))
  1230.             (null (equal (cdar *var*) '(1 1))))
  1231.             (merror "Second arg to POLYDECOMP must be an atom"))
  1232.            (t (setq *var* (caar *var*))))
  1233.          (cond ((or (pcoefp (cdr p))
  1234.             (null (eq (cadr p) *var*)))
  1235.             (setq den (cdr p)
  1236.               p (car p)))
  1237.            (t (merror "Cannot POLYDECOMP a rational function")))
  1238.          (cons '(mlist)
  1239.            (cond ((or (pcoefp p)
  1240.                   (null (eq (car p) *var*)))
  1241.               (list (rdis (cons p den))))
  1242.              (t (setq p (pdecomp p *var*))
  1243.                 (do ((l
  1244.                   (setq p (mapcar #'(lambda (q) (cons q 1)) p))
  1245.                   (cdr l))
  1246.                  (a))
  1247.                 ((null l)
  1248.                  (cons (rdis (cons (caar p)
  1249.                            (ptimes (cdar p) den)))
  1250.                        (mapcar #'rdis (cdr p))))
  1251.                 (cond ((setq a (pdecpow (car l) *var*))
  1252.                        (rplaca l (car a))
  1253.                        (cond ((cdr l)
  1254.                           (rplacd l
  1255.                               (cons (ratplus
  1256.                                  (rattimes
  1257.                                   (cadr l)
  1258.                                   (cons (pterm (cdaadr a) 1)
  1259.                                     (cdadr a))
  1260.                                   t)
  1261.                                  (cons
  1262.                                   (pterm (cdaadr a) 0)
  1263.                                   (cdadr a)))
  1264.                                 (cddr l))))
  1265.                          ((equal (cadr a)
  1266.                              (cons (list *var* 1 1) 1)))
  1267.                          (t (rplacd l (list (cadr a)))))))))))))
  1268.  
  1269.  
  1270. ;;; POLYDECOMP is like $POLYDECOMP except it takes a poly in *POLY* format (as
  1271. ;;; defined in SOLVE) (numerator of a RAT form) and returns a list of
  1272. ;;; headerless rat forms.  In otherwords, it is $POLYDECOMP minus type checking
  1273. ;;; and conversions to/from general representation which SOLVE doesn't
  1274. ;;; want/need on a general basis.
  1275. ;;; It is used in the SOLVE package and as such it should have an autoload
  1276. ;;; property 
  1277.  
  1278. (defun polydecomp (p *var*)
  1279.        (let ($FACTORFLAG $RATFAC)
  1280.          (cond ((or (pcoefp p)
  1281.             (null (eq (car p) *var*)))
  1282.             (cons p nil))
  1283.            (t (setq p (pdecomp p *var*))
  1284.               (do ((l (setq p
  1285.                     (mapcar #'(lambda (q) (cons q 1))
  1286.                         p))
  1287.                   (cdr l))
  1288.                (a))
  1289.               ((null l)
  1290.                (cons (cons (caar p)
  1291.                        (cdar p))
  1292.                  (cdr p)))
  1293.               (cond ((setq a (pdecpow (car l) *var*))
  1294.                  (rplaca l (car a))
  1295.                  (cond ((cdr l)
  1296.                     (rplacd l
  1297.                         (cons (ratplus
  1298.                                (rattimes
  1299.                             (cadr l)
  1300.                             (cons (pterm (cdaadr a) 1)
  1301.                                   (cdadr a))
  1302.                             t)
  1303.                                (cons
  1304.                             (pterm (cdaadr a) 0)
  1305.                             (cdadr a)))
  1306.                               (cddr l))))
  1307.                        ((equal (cadr a)
  1308.                            (cons (list *var* 1 1) 1)))
  1309.                        (t (rplacd l (list (cadr a))))))))))))
  1310.  
  1311.  
  1312.  
  1313. (defun pdecred (f h *var*)                ;f = g(h(*var*))
  1314.        (cond ((or (pcoefp h) (null (eq (car h) *var*))
  1315.           (equal (cadr h) 1)
  1316.           (null (zerop (remainder (cadr f) (cadr h))))
  1317.           (and (null (pzerop (caadr (setq f (pdivide f h)))))
  1318.                (equal (cdadr f) 1)))
  1319.           nil)
  1320.          (t (do ((q (pdivide (caar f) h) (pdivide (caar q) h))
  1321.              (i 1 (f1+ i))
  1322.              (*ans*))
  1323.             ((pzerop (caar q))
  1324.              (cond ((and (equal (cdadr q) 1)
  1325.                  (or (pcoefp (caadr q))
  1326.                      (null (eq (caar (cadr q)) *var*))))
  1327.                 (psimp *var* (cons i (cons (caadr q) *ans*))))))
  1328.             (cond ((and (equal (cdadr q) 1)
  1329.                 (or (pcoefp (caadr q))
  1330.                     (null (eq (caar (cadr q)) *var*))))
  1331.                (and (null (pzerop (caadr q)))
  1332.                 (setq *ans* (cons i (cons (caadr q) *ans*)))))
  1333.               (t (return nil)))))))
  1334.  
  1335. (defun pdecomp (p *var*)
  1336.        (let ((c (pterm (cdr p) 0))
  1337.           (a) (*x* (list *var* 1 1)))
  1338.          (cons (pcplus c (car (setq a (pdecomp* (pdifference p c)))))
  1339.            (cdr a))))
  1340.  
  1341. (defun pdecomp* (*p*)
  1342.        (let ((a)
  1343.           (l (pdecgdfrm (pfactor (pquotient *p* *x*)))))
  1344.          (cond ((or (pdecprimep (cadr *p*))
  1345.             (null (setq a (pdecomp1 *x* l))))
  1346.             (list *p*))
  1347.            (t (append (pdecomp* (car a)) (cdr a))))))
  1348.  
  1349. (defun pdecomp1 (prod l)
  1350.        (cond ((null l)
  1351.           (and (null (equal (cadr prod) (cadr *p*)))
  1352.            (setq l (pdecred *p* prod *var*))
  1353.            (list l prod)))
  1354.          ((pdecomp1 prod (cdr l)))
  1355.          (t (pdecomp1 (ptimes (car l) prod) (cdr l)))))
  1356.  
  1357. (defun pdecgdfrm (l)                    ;Get list of divisors
  1358.        (do ((l (copy-top-level l ))
  1359.         (ll (list (car l))
  1360.         (cons (car l) ll)))
  1361.        (nil)
  1362.        (rplaca (cdr l) (f1- (cadr l)))
  1363.        (cond ((signp e (cadr l))
  1364.           (setq l (cddr l))))
  1365.        (cond ((null l) (return ll)))))
  1366.  
  1367. (defun pdecprimep (x)
  1368.        (setq x (cfactorw x))
  1369.        (and (null (cddr x)) (equal (cadr x) 1))) 
  1370.        
  1371. (defun pdecpow (p *var*)
  1372.        (setq p (car p))
  1373.        (let ((p1 (pderivative p *var*))
  1374.           p2 p1p p1c a lin p2p)
  1375.          (setq p1p (oldcontent p1)
  1376.            p1c (car p1p) p1p (cadr p1p))
  1377.          (setq p2 (pderivative p1 *var*))
  1378.          (setq p2p (cadr (oldcontent p2)))
  1379.          (and (setq lin (testdivide p1p p2p))
  1380.           (null (pcoefp lin))
  1381.           (eq (car lin) *var*)
  1382.           (list (ratplus
  1383.              (rattimes (cons (list *var* (cadr p) 1) 1)
  1384.                    (setq a (ratreduce p1c
  1385.                               (ptimes (cadr p)
  1386.                                   (caddr lin))))
  1387.                    t)
  1388.              (ratdif (cons p 1)
  1389.                  (rattimes a (cons (pexpt lin (cadr p)) 1)
  1390.                        t)))
  1391.             (cons lin 1)))))
  1392.  
  1393. #-NIL
  1394. (declare-top(unspecial *mfactl *factlist donel nn* dn* splist ans var dict *var* *ans*
  1395.             a* *a *n *a*  hi lo
  1396.             *infsumsimp *times *plus sum usum makef gensim))
  1397.  
  1398. ;; (declare (unspecial *mfactl *factlist donel nn* dn* splist *ans* *var*  dict
  1399. ;;             a* *a *n *a* $prevfib hi lo
  1400. ;;             *infsumsimp *times *plus sum usum makef gensim))
  1401.  
  1402.  
  1403.