src/combin.lisp

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