src/combin.lisp

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