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