src/risch.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 risch)
  14
  15 (load-macsyma-macros rzmac ratmac)
  16
  17 (declare-top (special *mosesflag
  18                       context *in-risch-p*))
  19
  20 (defmvar $erfflag t "Controls whether `risch' generates `erfs'")
  21
  22 (defvar *changevp* t
  23   "When nil prevents changevar hack")
  24
  25 (defmacro pair (al bl) `(mapcar #'cons ,al ,bl))
  26
  27 ;; internal representation of risch expressions: list with canonical rational
  28 ;; expression (CRE) as first element, standard maxima expressions as remaining
  29 ;; elements.  risch expression is sum of CRE and remaining elements.
  30 (defmacro rischzero () ''((0 . 1) 0))
  31
  32 (defun rischnoun (exp1 risch-ratform risch-intvar &optional (exp2 exp1 exp2p))
  33   (unless exp2p (setq exp1 (rzero)))
  34   `(,exp1 ((%integrate) ,(disrep exp2 risch-ratform) ,risch-intvar)))
  35
  36 (defun getrischvar ()
  37   (do ((vl varlist (cdr vl))
  38        (gl genvar (cdr gl)))
  39       ((null (cdr vl)) (car gl))))
  40
  41 ;; test whether CRE p is constant with respect to variable of integration.
  42 ;; requires variables in varlist and genvar
  43 ;; to be ordered as by intsetup, with var of integration ordered before
  44 ;; any other expressions that contain it.
  45 (defun risch-pconstp (p risch-mainvar)
  46   (or (pcoefp p) (pointergp risch-mainvar (car p))))
  47
  48 (defun risch-constp (r risch-mainvar)
  49   (setq r (ratfix r))
  50   (and (risch-pconstp (car r) risch-mainvar)
  51        (risch-pconstp (cdr r) risch-mainvar)))
  52
  53 ;; adds two risch expressions (defined above).
  54 (defun rischadd (x y)
  55   (destructuring-let (((a . b) x) ((c . d) y))
  56     (cons (r+ a c) (append b d))))
  57
  58 (defmfun $risch (exp risch-var)
  59   (let ((*integrator-level* 0))
  60     (declare (special *integrator-level*))
  61     (with-new-context (context)
  62       (rischint exp risch-var))))
  63
  64 (defun spderivative (p risch-var)
  65   (cond ((pcoefp p) '(0 . 1))
  66         ((null (cdr p)) '(0 . 1))
  67         ((or (not (atom (car p))) (numberp (car p))) ;P IS A RATFORM
  68          (let ((denprime (spderivative (cdr p) risch-var)))
  69            (cond ((rzerop denprime)
  70                   (ratqu (spderivative (car p) risch-var) (cdr p)))
  71                  (t (ratqu (r- (r* (spderivative (car p) risch-var)
  72                                    (cdr p))
  73                                (r* (car p) denprime))
  74                            (r* (cdr p) (cdr p)))))))
  75         (t (r+ (spderivative1 (car p)
  76                               (cadr p)
  77                               (caddr p)
  78                               risch-var)
  79                (spderivative (cons (car p) (cdddr p))
  80                              risch-var)))))
  81
  82 (defun spderivative1 (var1 deg coeff risch-var)
  83   (cond ((eq var1 risch-var)
  84          (r* (ratexpt (cons (list risch-var 1 1) 1) (1- deg))
  85              (pctimes deg coeff)))
  86         ((pointergp risch-var var1) '(0 . 1))
  87         ((equal deg 0) (spderivative coeff risch-var))
  88         (t (r+ (r* (ratexpt (cons (list var1 1 1) 1) deg)
  89                    (spderivative coeff risch-var))
  90                (r* (cond ((equal deg 1) coeff)
  91                          (t (r* deg
  92                                 coeff
  93                                 (ratexpt (cons (list var1 1 1) 1)
  94                                          (1- deg)))))
  95                    (get var1 'rischdiff) )))))
  96
  97 (defun polylogp (exp &optional sub)
  98   (and (mqapplyp exp) (eq (subfunname exp) '$li)
  99        (or (null sub) (equal sub (car (subfunsubs exp))))))
 100
 101 (defun rischint (exp risch-intvar &aux ($logarc nil) ($exponentialize nil)
 102                      ($gcd '$algebraic) ($algebraic t) (implicit-real t)
 103                      ($float nil) ($numer nil)
 104                      ;; The risch integrator expects $logexpand T. Otherwise,
 105                      ;; the integrator hangs for special types of integrals
 106                      ;; (See bug report ID:3039452)
 107                      ($logexpand t))
 108   (prog ($%e_to_numlog $logsimp risch-y z risch-var risch-ratform risch-liflag
 109          risch-mainvar varlist genvar $ratfac $ratalgdenom risch-degree
 110          rischform-value risch-trigint risch-hypertrigint risch-operator)
 111      (if (specrepp exp)
 112          (setq exp (specdisrep exp)))
 113      (if (specrepp risch-intvar)
 114          (setq risch-intvar (specdisrep risch-intvar)))
 115      (if (mnump risch-intvar)
 116          (merror (intl:gettext "risch: attempt to integrate wrt a number: ~:M") risch-intvar))
 117      (if (and (atom risch-intvar)
 118               (isinop exp risch-intvar))
 119          (go noun))
 120      (multiple-value-setq (rischform-value risch-trigint risch-hypertrigint risch-operator)
 121        (rischform exp risch-intvar))
 122      (cond (risch-trigint
 123             (return (trigin1 exp risch-intvar)))
 124            (risch-hypertrigint
 125             (return (hypertrigint1 exp risch-intvar t)))
 126            (risch-operator
 127             (go noun)))
 128      (multiple-value-setq (risch-y risch-operator)
 129        (intsetup exp risch-intvar))
 130      (if risch-operator
 131          (go noun))
 132      (setq risch-ratform (car risch-y))
 133      (setq varlist (caddr risch-ratform))
 134      (setq risch-mainvar (caadr (ratf risch-intvar)))
 135      (setq genvar (cadddr risch-ratform))
 136      (unless (some #'algpget varlist)
 137        (setq $algebraic nil)
 138        (setq $gcd (car *gcdl*)))
 139      (setq risch-var (getrischvar))
 140      (setq z (tryrisch (cdr risch-y) risch-mainvar risch-ratform risch-intvar risch-liflag risch-degree risch-var))
 141      (setf (caddr risch-ratform) varlist)
 142      (setf (cadddr risch-ratform) genvar)
 143      (return (cond ((atom (cdr z))
 144                     (disrep (car z) risch-ratform))
 145                    (t
 146                     (let (($logsimp t)
 147                           ($%e_to_numlog t))
 148                       (simplify (list* '(mplus)
 149                                        (disrep (car z) risch-ratform)
 150                                        (cdr z)))))))
 151    noun
 152      (return (list '(%integrate) exp risch-intvar))))
 153
 154 (defun rischform (l risch-intvar)
 155   (let (risch-trigint risch-hypertrigint risch-operator)
 156     (labels
 157         ((rischform-impl (l risch-intvar)
 158            (cond ((or (atom l)
 159                       (alike1 risch-intvar l)
 160                       (freeof risch-intvar l))
 161                   nil)
 162                  ((polylogp l)
 163                   (if (and (integerp (car (subfunsubs l)))
 164                            (signp g (car (subfunsubs l))))
 165                       (rischform-impl (car (subfunargs l)) risch-intvar)
 166                       (setq risch-operator t)))
 167                  ((atom (caar l))
 168                   (case (caar l)
 169                     ((%sin %cos %tan %cot %sec %csc)
 170                      (setq risch-trigint t $exponentialize t)
 171                      (rischform-impl (cadr l) risch-intvar))
 172                     ((%asin %acos %atan %acot %asec %acsc)
 173                      (setq risch-trigint t $logarc t)
 174                      (rischform-impl (cadr l) risch-intvar))
 175                     ((%sinh %cosh %tanh %coth %sech %csch)
 176                      (setq risch-hypertrigint t $exponentialize t)
 177                      (rischform-impl (cadr l) risch-intvar))
 178                     ((%asinh %acosh %atanh %acoth %asech %acsch)
 179                      (setq risch-hypertrigint t $logarc t)
 180                      (rischform-impl (cadr l) risch-intvar))
 181                     ((mtimes mplus mexpt rat %erf %log)
 182                      (mapc #'(lambda (e)
 183                                (rischform-impl e risch-intvar))
 184                            (cdr l)))
 185                     (t
 186                      (setq risch-operator (caar l)))))
 187                  (t
 188                   (setq risch-operator (caar l))))))
 189       (values (rischform-impl l risch-intvar)
 190               risch-trigint
 191               risch-hypertrigint
 192               risch-operator))))
 193
 194 (defun hypertrigint1 (exp risch-var hyperfunc)
 195   (let ((result (if hyperfunc
 196                     (sinint (resimplify exp) risch-var)
 197                     (rischint (resimplify exp) risch-var))))
 198     ;; The result can contain solvable integrals. Look for this case.
 199     (if (isinop result '%integrate)
 200         ;; Found an integral. Evaluate the result again.
 201         ;; Set the flag *in-risch-p* to make sure that we do not call
 202         ;; rischint again from the integrator. This avoids endless loops.
 203         (let ((*in-risch-p* t))
 204           (meval (list '($ev) result '$nouns)))
 205         result)))
 206
 207 (defun trigin1 (risch-*exp risch-var)
 208   (let ((yyy (hypertrigint1 risch-*exp risch-var nil)))
 209     (setq yyy (div ($expand ($num yyy))
 210                    ($expand ($denom yyy))))
 211     (let ((rischp risch-var)
 212           (rp-polylogp t)
 213           $logarc $exponentialize result)
 214       (setq result (sratsimp (if (and (freeof '$%i risch-*exp) (freeof '$li yyy))
 215                                  ($realpart yyy)
 216                                  ($rectform yyy))))
 217       ;; The result can contain solvable integrals. Look for this case.
 218       (if (isinop result '%integrate)
 219           ;; Found an integral. Evaluate the result again.
 220           ;; Set the flag *in-risch-p* to make sure that we do not call
 221           ;; rischint again from the integrator. This avoids endless loops.
 222           (let ((*in-risch-p* t))
 223             (meval (list '($ev) result '$nouns)))
 224           result))))
 225
 226 (defun tryrisch (exp risch-mainvar risch-ratform risch-intvar risch-liflag risch-degree risch-var)
 227   (prog (risch-logptdx risch-expflag risch-expstuff risch-expint risch-y)
 228      (setq risch-expstuff '(0 . 1))
 229      (cond ((eq risch-mainvar risch-var)
 230             (return (rischfprog exp risch-ratform risch-var)))
 231            ((eq (get risch-var 'leadop)
 232                 'mexpt)
 233             (setq risch-expflag t)))
 234      (multiple-value-setq (risch-y risch-logptdx risch-expint)
 235        (rischlogdprog exp risch-ratform risch-intvar risch-liflag risch-var
 236                       risch-expflag risch-mainvar risch-expint risch-degree))
 237      (dolist (rat risch-logptdx)
 238        (let (rischlogeprog-value)
 239          (setq risch-y
 240                (rischadd (multiple-value-setq (rischlogeprog-value risch-expint)
 241                            (rischlogeprog rat risch-ratform nil risch-intvar risch-expstuff
 242                                           risch-var risch-expflag risch-mainvar risch-expint))
 243                          risch-y))))
 244      (if varlist
 245          (setq risch-y (rischadd (tryrisch1 risch-expstuff risch-mainvar
 246                                             risch-ratform risch-intvar risch-liflag
 247                                             risch-degree)
 248                            risch-y)))
 249      (return (if risch-expint
 250                  (rischadd (rischexppoly risch-expint risch-var
 251                                          risch-ratform risch-intvar risch-liflag
 252                                          risch-degree risch-mainvar)
 253                            risch-y)
 254                  risch-y))))
 255
 256 (defun tryrisch1 (exp risch-mainvar risch-ratform risch-intvar risch-liflag risch-degree)
 257   (let* ((varlist (reverse (cdr (reverse varlist))))
 258          (risch-var (getrischvar)))
 259     (tryrisch exp risch-mainvar risch-ratform risch-intvar risch-liflag risch-degree risch-var)))
 260
 261 (defun rischfprog (rat risch-ratform risch-var)
 262   (multiple-value-bind (dprog-ret risch-logptdx)
 263       (dprog rat risch-ratform risch-var)
 264     (cons (cdr (ratrep* dprog-ret))
 265           (let ((varlist varlist)
 266                 (genvar (subseq genvar 0 (length varlist))))
 267             (mapcar #'(lambda (p)
 268                         (eprog p risch-ratform risch-var nil))
 269                     risch-logptdx)))))
 270
 271 (defun rischlogdprog (ratarg risch-ratform risch-intvar risch-liflag risch-var risch-expflag risch-mainvar risch-expint risch-degree)
 272   (prog (arootf deriv thebpg thetop thebot prod1 prod2 ans
 273          risch-wholepart risch-logptdx risch-parnumer risch-pardenom
 274          risch-rootfactor rischlogpoly-value)
 275      (setq ans '(0 . 1))
 276      (cond ((or (pcoefp (cdr ratarg))
 277                 (pointergp risch-var (cadr ratarg)))
 278             (multiple-value-setq (rischlogpoly-value risch-expint)
 279               (rischlogpoly ratarg risch-ratform risch-intvar risch-liflag risch-var risch-expflag risch-mainvar risch-expint risch-degree))
 280             (return (values rischlogpoly-value
 281                             risch-logptdx
 282                             risch-expint))))
 283
 284      (multiple-value-setq (risch-rootfactor risch-pardenom)
 285        (aprog (ratdenominator ratarg) risch-var))
 286      (multiple-value-setq (risch-parnumer risch-wholepart)
 287        (cprog (ratnumerator ratarg)
 288               (ratdenominator ratarg)
 289               risch-var
 290               risch-pardenom))
 291      (do ((risch-rootfactor (reverse risch-rootfactor) (cdr risch-rootfactor))
 292           (risch-parnumer (reverse risch-parnumer) (cdr risch-parnumer))
 293           (risch-klth (length risch-rootfactor) (1- risch-klth)))
 294          ((= risch-klth 1))
 295        (setq arootf (car risch-rootfactor))
 296        (cond
 297          ((pcoefp arootf))
 298          ((and (eq (get (car arootf) 'leadop) 'mexpt)
 299                (null (cdddr arootf)))
 300           (setq
 301            risch-expint
 302            (append
 303             (cond ((and (not (atom (car risch-parnumer)))
 304                         (not (atom (caar risch-parnumer)))
 305                         (eq (caaar risch-parnumer) (car arootf)))
 306                    (gennegs arootf (cdaar risch-parnumer) (cdar risch-parnumer) risch-klth))
 307                   (t (list
 308                       (list 'neg (car risch-parnumer)
 309                             (car arootf) risch-klth (cadr arootf)))))
 310             risch-expint)))
 311          ((not (zerop (pdegree arootf risch-var)))
 312           (setq deriv (spderivative arootf risch-mainvar))
 313           (setq thebpg (bprog arootf (ratnumerator deriv) risch-var))
 314           (setq thetop (car risch-parnumer))
 315           (do ((kx (1- risch-klth) (1- kx))) ((= kx 0))
 316             (setq prod1 (r* thetop (car thebpg)))
 317             (setq prod2 (r* thetop (cdr thebpg) (ratdenominator deriv)))
 318             (setq thebot (pexpt arootf kx))
 319             (setq ans (r+ ans (ratqu (r- prod2) (r* kx thebot))))
 320             (setq thetop
 321                   (r+ prod1 (ratqu (spderivative prod2 risch-mainvar) kx)))
 322             (setq thetop (cdr (ratdivide thetop thebot))))
 323           (push (ratqu thetop arootf) risch-logptdx))))
 324      (push (ratqu (car risch-parnumer) (car risch-rootfactor)) risch-logptdx)
 325      (cond ((or (pzerop ans) (pzerop (car ans)))
 326             (multiple-value-setq (rischlogpoly-value risch-expint)
 327               (rischlogpoly risch-wholepart risch-ratform risch-intvar risch-liflag risch-var risch-expflag risch-mainvar risch-expint risch-degree))
 328             (return (values rischlogpoly-value
 329                             risch-logptdx
 330                             risch-expint))))
 331      (setq thetop (cadr (pdivide (ratnumerator ans)
 332                                  (ratdenominator ans))))
 333      (multiple-value-setq (rischlogpoly-value risch-expint)
 334        (rischlogpoly risch-wholepart risch-ratform risch-intvar risch-liflag risch-var risch-expflag risch-mainvar risch-expint risch-degree))
 335      (return (values (rischadd (ncons (ratqu thetop (ratdenominator ans)))
 336                                rischlogpoly-value)
 337                      risch-logptdx
 338                      risch-expint))))
 339
 340 (defun gennegs (denom num numdenom risch-klth)
 341   (cond ((null num) nil)
 342         (t (cons (list 'neg (cadr num)
 343                        (car denom)
 344                        (- risch-klth (car num))
 345                        (r* numdenom (caddr denom) ))
 346                  (gennegs denom (cddr num) numdenom risch-klth)))))
 347
 348 (defun rischlogeprog (p risch-ratform risch-switch1 risch-intvar risch-expstuff
 349                       risch-var risch-expflag risch-mainvar risch-expint)
 350   (labels
 351       ((impl (p risch-switch1)
 352          (prog (p1e p2e p2deriv logcoef ncc dcc allcc expcoef my-divisor
 353                 risch-parnumer risch-pardenom)
 354             (if (or (pzerop p) (pzerop (car p)))
 355                 (return (rischzero)))
 356             (setq p1e (ratnumerator p))
 357             (desetq (dcc p2e) (oldcontent (ratdenominator p)))
 358             (cond ((and (not risch-switch1)
 359                         (cdr (setq risch-pardenom (intfactor p2e))))
 360                    (setq risch-parnumer nil)
 361                    (setq risch-switch1 t)
 362                    (desetq (ncc p1e) (oldcontent p1e))
 363                    (multiple-value-setq (risch-parnumer)
 364                      (cprog p1e p2e risch-var risch-pardenom))
 365                    (setq allcc (ratqu ncc dcc))
 366                    (return (do ((pnum risch-parnumer (cdr pnum))
 367                                 (pden risch-pardenom (cdr pden))
 368                                 (ans (rischzero)))
 369                                ((or (null pnum) (null pden))
 370                                 (setq risch-switch1 nil) ans)
 371                              (setq ans (rischadd
 372                                         (impl
 373                                          (r* allcc (ratqu (car pnum) (car pden)))
 374                                          risch-switch1)
 375                                         ans))))))
 376             (when (and risch-expflag (null (p-red p2e)))
 377               (push (cons 'neg p) risch-expint)
 378               (return (rischzero)))
 379             (if risch-expflag
 380                 (setq expcoef (r* (p-le p2e) (ratqu (get risch-var 'rischdiff)
 381                                                     (make-poly risch-var)))))
 382             (setq p1e (ratqu p1e (ptimes dcc (p-lc p2e)))
 383                   p2e (ratqu p2e (p-lc p2e))) ;MAKE DENOM MONIC
 384             (setq p2deriv (spderivative p2e risch-mainvar))
 385             (setq my-divisor (if risch-expflag
 386                                  (r- p2deriv (r* p2e expcoef))
 387                                  p2deriv))
 388             (when (equal my-divisor '(0 . 1))
 389               ;; (format t "HEY RISCHLOGEPROG, FOUND ZERO DIVISOR; GIVE UP.~%")
 390               (return (rischnoun p risch-ratform risch-intvar)))
 391             (setq logcoef (ratqu p1e my-divisor))
 392             (when (risch-constp logcoef risch-mainvar)
 393               (if risch-expflag
 394                   (setq risch-expstuff (r- risch-expstuff (r* expcoef logcoef))))
 395               (return
 396                 (list
 397                  '(0 . 1)
 398                  (list '(mtimes)
 399                        (disrep logcoef risch-ratform)
 400                        (logmabs (disrep p2e risch-ratform))))))
 401             (if (and risch-expflag
 402                      $liflag
 403                      *changevp*)
 404                 (let* ((newvar (gensym))
 405                        (new-int ($changevar
 406                                  `((%integrate) ,(simplify (disrep p risch-ratform)) ,risch-intvar)
 407                                  (sub newvar (get risch-var 'rischexpr))
 408                                  newvar risch-intvar))
 409                        (*changevp* nil)) ;prevents recursive changevar
 410                   (if (and (freeof risch-intvar new-int)
 411                            (freeof '%integrate
 412                                    (setq new-int (rischint (sdiff new-int newvar)
 413                                                            newvar))))
 414                       (return
 415                         (list (rzero)
 416                               (maxima-substitute (get risch-var 'rischexpr) newvar new-int))))))
 417             (return (rischnoun p risch-ratform risch-intvar)))))
 418     (values (impl p risch-switch1)
 419             risch-expint)))
 420
 421
 422 (defun findint (exp)
 423   (cond ((atom exp) nil)
 424         ((atom (car exp)) (findint (cdr exp)))
 425         ((eq (caaar exp) '%integrate) t)
 426         (t (findint (cdr exp)))))
 427
 428 (defun logequiv (fn1 fn2 risch-intvar)
 429   (freeof risch-intvar ($ratsimp (div* (remabs (leadarg fn1))
 430                                  (remabs (leadarg fn2))))))
 431
 432 (defun remabs (exp)
 433   (cond ((atom exp) exp)
 434         ((eq (caar exp) 'mabs) (cadr exp))
 435         (t exp)))
 436
 437 (defun getfnsplit (l risch-intvar)
 438   (let (coef fn)
 439     (dolist (x l (values (muln coef nil) (muln fn nil)))
 440       (if (free x risch-intvar)
 441           (push x coef)
 442           (push x fn)))))
 443
 444 (defun getfncoeff (a form risch-intvar risch-liflag risch-degree risch-cary risch-nogood risch-lians)
 445   (labels
 446       ((getfncoeff-impl (a)
 447          (cond ((null a) 0)
 448                ((equal (car a) 0)
 449                 (getfncoeff-impl (cdr a)))
 450                ((and (listp (car a))
 451                      (eq (caaar a) 'mplus)
 452                      (ratpl (getfncoeff-impl (cdar a))
 453                             (getfncoeff-impl (cdr a)))))
 454                ((and (listp (car a))
 455                      (eq (caaar a) 'mtimes))
 456                 (multiple-value-bind (coef newfn)
 457                     (getfnsplit (cdar a) risch-intvar)
 458                   ;; (car a) is a mtimes expression. We insert coef and newfn as the
 459                   ;; new arguments to the mtimes expression. This causes problems if
 460                   ;;   (1) coef is a mtimes expression too and
 461                   ;;   (2) (car a) has already a simp flag
 462                   ;; We get a nested mtimes expression, which does not sgetfncoeff-implify.
 463                   ;; We comment out the following code (DK 09/2009):
 464                   ;; (setf (cdar a) (list coef newfn))
 465
 466                   ;; Insert a complete mtimes expression without simpflag.
 467                   ;; Nested mtimes expressions sgetfncoeff-implify further.
 468                   (setf (car a) (list '(mtimes) coef newfn))
 469
 470                   (setf (cdar a) (list coef newfn))
 471                   (cond ((zerop1 coef)
 472                          (getfncoeff-impl (cdr a)))
 473                         ((and (matanp newfn)
 474                               (member '$%i varlist :test #'eq))
 475                          (let (($logarc t) ($logexpand '$all))
 476                            (rplaca a ($expand (resimplify (car a)))))
 477                          (getfncoeff-impl a))
 478                         ((and (alike1 (leadop newfn) (leadop form))
 479                               (or (alike1 (leadarg newfn) (leadarg form))
 480                                   (and (mlogp newfn)
 481                                        (logequiv form newfn risch-intvar))))
 482                          (ratpl (rform coef)
 483                                 (prog2 (rplaca a 0)
 484                                     (getfncoeff-impl (cdr a)))))
 485                         ((do ((vl varlist (cdr vl)))
 486                              ((null vl))
 487                            (and (not (atom (car vl)))
 488                                 (alike1 (leadop (car vl)) (leadop newfn))
 489                                 (if (mlogp newfn)
 490                                     (logequiv (car vl) newfn risch-intvar)
 491                                     (alike1 (car vl) newfn))
 492                                 (rplaca (cddar a) (car vl))
 493                                 (return nil))))
 494                         ((let (vlist)
 495                            (declare (special vlist))
 496                            ;; newvar1 accesses the special var vlist.
 497                            (newvar1 (car a))
 498                            (null vlist))
 499                          (setq risch-cary
 500                                (ratpl (cdr (ratrep* (car a)))
 501                                       risch-cary))
 502                          (rplaca a 0)
 503                          (getfncoeff-impl (cdr a)))
 504                         ((and risch-liflag
 505                               (mlogp form)
 506                               (mlogp newfn))
 507                          (let (res)
 508                            (multiple-value-setq (res risch-nogood)
 509                              (dilog (cons (car a) form) risch-intvar risch-degree risch-nogood))
 510                            (push res risch-lians))
 511                          (rplaca a 0)
 512                          (getfncoeff-impl (cdr a)))
 513                         ((and risch-liflag
 514                               (polylogp form)
 515                               (mlogp newfn)
 516                               (logequiv form newfn risch-intvar))
 517                          (push (mul* (cadar a) (make-li (1+ (car (subfunsubs form)))
 518                                                         (leadarg form)))
 519                                risch-lians)
 520                          (rplaca a 0)
 521                          (getfncoeff-impl (cdr a)))
 522                         (t (setq risch-nogood t) 0))))
 523                (t
 524                 (rplaca a (list '(mtimes) 1 (car a)))
 525                 (getfncoeff-impl a)))))
 526     (values (getfncoeff-impl a) risch-cary risch-nogood
 527             risch-lians)))
 528
 529
 530 (defun rischlogpoly (exp risch-ratform risch-intvar risch-liflag risch-var risch-expflag risch-mainvar risch-expint risch-degree)
 531   (let
 532       ((result
 533          (cond ((equal exp '(0 . 1))
 534                 (rischzero))
 535                (risch-expflag
 536                 (push (cons 'poly exp) risch-expint)
 537                 (rischzero))
 538                ((not (among risch-var exp))
 539                 (tryrisch1 exp risch-mainvar risch-ratform risch-intvar risch-liflag risch-degree))
 540                (t
 541                 (do ((risch-degree (pdegree (car exp) risch-var) (1- risch-degree))
 542                      (p (car exp))
 543                      (den (cdr exp))
 544                      (risch-lians ())
 545                      (sum (rzero))
 546                      (risch-cary (rzero))
 547                      (risch-y)
 548                      (z)
 549                      (ak)
 550                      (risch-nogood)
 551                      (lbkpl1))
 552                     ((minusp risch-degree)
 553                      (cons sum (append risch-lians (cdr risch-y))))
 554                   (setq ak (r- (ratqu (polcoef p risch-degree risch-var) den)
 555                                (r* (cons (1+ risch-degree) 1)
 556                                    risch-cary
 557                                    (get risch-var 'rischdiff))))
 558                   (if (not (pzerop (polcoef p risch-degree risch-var)))
 559                       (setq p (if (pcoefp p)
 560                                   (pzero)
 561                                   (psimp risch-var (p-red p)))))
 562                   (setq risch-y (tryrisch1 ak risch-mainvar
 563                                            risch-ratform risch-intvar risch-liflag risch-degree))
 564                   (setq risch-cary (car risch-y))
 565                   (and (> risch-degree 0)
 566                        (setq risch-liflag $liflag))
 567
 568                   (multiple-value-setq (z risch-cary risch-nogood risch-lians)
 569                     (getfncoeff (cdr risch-y)
 570                                 (get risch-var 'rischexpr)
 571                                 risch-intvar risch-liflag risch-degree risch-cary risch-nogood
 572                                 risch-lians))
 573                   (setq risch-liflag nil)
 574                   (cond ((and (> risch-degree 0)
 575                               (or risch-nogood
 576                                   (findint (cdr risch-y))))
 577                          (return (rischnoun sum risch-ratform
 578                                             risch-intvar
 579                                             (r+ (r* ak
 580                                                     (make-poly risch-var risch-degree 1))
 581                                                 (ratqu p den))))))
 582                   (setq lbkpl1 (ratqu z (cons (1+ risch-degree) 1)))
 583                   (setq sum (r+ (r* lbkpl1 (make-poly risch-var (1+ risch-degree) 1))
 584                                 (r* risch-cary (if (zerop risch-degree) 1
 585                                                    (make-poly risch-var risch-degree 1)))
 586                                 sum)))))))
 587     (values result risch-expint)))
 588
 589 (defun make-li (sub arg)
 590   (subfunmake '$li (ncons sub) (ncons arg)))
 591
 592 ;;integrates log(ro)^risch-degree*log(rn)' in terms of polylogs
 593 ;;finds constants c,d and integers j,k such that
 594 ;;c*ro^j+d=rn^k  If ro and rn are poly's then can assume either j=1 or k=1
 595 (defun dilog (l risch-intvar risch-degree risch-nogood)
 596   (destructuring-let* ((((nil coef nlog) . olog) l)
 597                        (narg (remabs (cadr nlog)))
 598                        (varlist varlist)
 599                        (genvar genvar)
 600                        (rn (rform narg))        ;; can add new vars to varlist
 601                        (ro (rform (cadr olog)))
 602                        (risch-var (caar ro))
 603                        ((j . k) (ratreduce (pdegree (car rn) risch-var) (pdegree (car ro) risch-var)))
 604                        (idx (gensym))
 605                        (rc) (rd))
 606     (cond ((and (= j 1) (> k 1))
 607            (setq rn (ratexpt rn k)
 608                  coef (div coef k)
 609                  narg (rdis rn)))
 610           ((and (= k 1) (> j 1))
 611            (setq ro (ratexpt ro j)
 612                  coef (div coef (f* j risch-degree))
 613                  olog (mul j olog))))
 614     (desetq (rc . rd) (ratdivide rn ro))
 615     (let
 616         ((result
 617           (cond ((and (freeof risch-intvar (rdis rc)) ;; can't use risch-constp because varlist
 618                       (freeof risch-intvar (rdis rd))) ;; is not set up with vars in correct order.
 619                  (setq narg ($ratsimp (sub 1 (div narg (rdis rd)))))
 620                  (mul* coef (power -1 (1+ risch-degree))
 621                        `((mfactorial) ,risch-degree)
 622                        (dosum (mul* (power -1 idx)
 623                                     (div* (power olog idx)
 624                                           `((mfactorial) ,idx))
 625                                     (make-li (add risch-degree (neg idx) 1) narg))
 626                               idx 0 risch-degree t)))
 627                 (t
 628                  (setq risch-nogood t)
 629                  0))))
 630       (values result risch-nogood))))
 631
 632 (defun exppolycontrol (flag f a expg n risch-ratform risch-intvar risch-liflag risch-degree risch-mainvar)
 633   (let (risch-y l risch-var (varlist varlist) (genvar genvar))
 634     (setq varlist (reverse (cdr (reverse varlist))))
 635     (setq risch-var (getrischvar))
 636     (setq risch-y (get risch-var 'leadop))
 637     (cond ((and (not (pzerop (ratnumerator f)))
 638                 (risch-constp (setq l (ratqu a f)) risch-mainvar))
 639            (cond (flag          ;; multiply in expg^n - n may be negative
 640                   (list (r* l (ratexpt (cons (list expg 1 1) 1) n))
 641                         0))
 642                  (t l)))
 643           ((eq risch-y risch-intvar)
 644            (rischexpvar flag (list f a expg n) risch-ratform risch-intvar risch-y risch-var risch-mainvar))
 645           (t
 646            (rischexplog (eq risch-y 'mexpt) flag f a
 647                         (list expg n (get risch-var 'rischarg)
 648                               risch-var (get risch-var 'rischdiff))
 649                         risch-ratform risch-intvar
 650                         risch-liflag
 651                         risch-degree
 652                         risch-y
 653                         risch-var
 654                         risch-mainvar)))))
 655
 656 (defun rischexppoly (risch-expint risch-var risch-ratform risch-intvar risch-liflag risch-degree risch-mainvar)
 657   (let (risch-y
 658         w
 659         num
 660         denom type
 661         (ans (rischzero))
 662         (expdiff (ratqu (get risch-var 'rischdiff) (list risch-var 1 1))))
 663     (do ((risch-expint risch-expint (cdr risch-expint)))
 664         ((null risch-expint)
 665          ans)
 666       (desetq (type . risch-y) (car risch-expint))
 667       (desetq (num . denom) (ratfix risch-y))
 668       (cond ((eq type 'neg)
 669              (setq w (exppolycontrol t
 670                                      (r* (- (cadr denom))
 671                                          expdiff)
 672                                      (ratqu num (caddr denom))
 673                                      risch-var
 674                                      (- (cadr denom))
 675                                      risch-ratform
 676                                      risch-intvar
 677                                      risch-liflag
 678                                      risch-degree
 679                                      risch-mainvar)))
 680             ((or (numberp num)
 681                  (not (eq (car num) risch-var)))
 682              (setq w (tryrisch1 risch-y risch-mainvar
 683                                 risch-ratform risch-intvar risch-liflag risch-degree)))
 684             (t
 685              (setq w (rischzero))
 686              (do ((num (cdr num) (cddr num)))
 687                  ((null num))
 688                (cond ((equal (car num) 0)
 689                       (setq w (rischadd
 690                                (tryrisch1 (ratqu (cadr num) denom)
 691                                           risch-mainvar
 692                                           risch-ratform risch-intvar risch-liflag risch-degree)
 693                                w)))
 694                      (t
 695                       (setq w (rischadd (exppolycontrol
 696                                          t
 697                                          (r* (car num) expdiff)
 698                                          (ratqu (cadr num) denom)
 699                                          risch-var
 700                                          (car num)
 701                                          risch-ratform
 702                                          risch-intvar
 703                                          risch-liflag
 704                                          risch-degree
 705                                          risch-mainvar)
 706                                         w)))))))
 707       (setq ans (rischadd w ans)))))
 708
 709 (defun rischexpvar (flag l risch-ratform risch-intvar risch-y risch-var risch-mainvar)
 710   (prog (lcm risch-m p risch-alphar risch-gamma delta r s
 711          tt denom k wl wv i ytemp ttemp yalpha f a expg n yn yd
 712          risch-beta)
 713      (desetq (f a expg n) l)
 714      (cond ((or (pzerop a) (pzerop (car a)))
 715             (return (cond ((null flag) (rzero))
 716                           (t (rischzero))))))
 717      (setq denom (ratdenominator f))
 718      (multiple-value-setq (p risch-alphar)
 719        (findpr (cdr (partfrac a risch-mainvar))
 720                (cdr (partfrac f risch-mainvar))
 721                risch-y
 722                risch-mainvar))
 723      (setq lcm (plcm (ratdenominator a) p))
 724      (setq risch-y (ratpl (spderivative (cons 1 p) risch-mainvar)
 725                           (ratqu f p)))
 726      (setq lcm (plcm lcm (ratdenominator risch-y)))
 727      (setq r (car (ratqu lcm p)))
 728      (setq s (car (r* lcm risch-y)))
 729      (setq tt (car (r* a lcm)))
 730      (setq risch-beta (pdegree r risch-mainvar))
 731      (setq risch-gamma (pdegree s risch-mainvar))
 732      (setq delta (pdegree tt risch-mainvar))
 733      (setq risch-alphar (max (- (1+ delta) risch-beta)
 734                              (- delta risch-gamma)))
 735      (setq risch-m 0)
 736      (cond ((equal (1- risch-beta) risch-gamma)
 737             (setq risch-y (r* -1
 738                               (ratqu (polcoef s risch-gamma risch-var)
 739                                      (polcoef r risch-beta risch-var))))
 740             (and (equal (cdr risch-y) 1)
 741                  (numberp (car risch-y))
 742                  (setq risch-m (car risch-y)))))
 743      (setq risch-alphar (max risch-alphar risch-m))
 744      (if (minusp risch-alphar)
 745          (return (if flag
 746                      (cxerfarg (rzero) expg n a risch-ratform risch-intvar risch-mainvar)
 747                      nil)))
 748      (cond ((not (and (equal risch-alphar risch-m)
 749                       (not (zerop risch-m))))
 750             (go down2)))
 751      (setq k (+ risch-alphar risch-beta -2))
 752      (setq wl nil)
 753    l2
 754      (setq wv (list (cons (polcoef tt k risch-var) 1)))
 755      (setq i risch-alphar)
 756    l1
 757      (setq wv
 758            (cons (r+ (r* (cons i 1)
 759                          (polcoef r (+ k 1 (- i)) risch-var))
 760                      (cons (polcoef s (+ k (- i)) risch-var) 1))
 761                  wv))
 762      (decf i)
 763      (cond ((> i -1)
 764             (go l1)))
 765      (setq wl (cons wv wl))
 766      (decf k)
 767      (cond ((> k -1)
 768             (go l2)))
 769      (multiple-value-setq (risch-y risch-m)
 770        (lsa wl))
 771      (if (or (eq risch-y 'singular)
 772              (eq risch-y 'inconsistent))
 773          (cond ((null flag)
 774                 (return nil))
 775                (t
 776                 (return (cxerfarg (rzero) expg n a risch-ratform risch-intvar risch-mainvar)))))
 777      (setq k 0)
 778      (setq lcm 0)
 779      (setq risch-y (cdr risch-y))
 780    l3
 781      (setq lcm
 782            (r+ (r* (car risch-y) (pexpt (list risch-mainvar 1 1) k))
 783                lcm))
 784      (incf k)
 785      (setq risch-y (cdr risch-y))
 786      (cond ((null risch-y)
 787             (return (cond ((null flag)
 788                            (ratqu lcm p))
 789                           (t
 790                            (list (r* (ratqu lcm p)
 791                                      (cons (list expg n 1) 1))
 792                                  0))))))
 793      (go l3)
 794    down2
 795      (cond ((> (1- risch-beta) risch-gamma)
 796             (setq k (+ risch-alphar (1- risch-beta)))
 797             (setq denom #'(lambda ()
 798                             (ratti risch-alphar (polcoef r risch-beta risch-var) t))))
 799            ((< (1- risch-beta) risch-gamma)
 800             (setq k (+ risch-alphar risch-gamma))
 801             (setq denom #'(lambda ()
 802                             (polcoef s risch-gamma risch-var))))
 803            (t
 804             (setq k (+ risch-alphar risch-gamma))
 805             (setq denom
 806                   #'(lambda ()
 807                       (ratpl (ratti risch-alphar (polcoef r risch-beta risch-var) t)
 808                              (polcoef s risch-gamma risch-var))))))
 809      (setq risch-y 0)
 810    loop
 811      (setq yn (polcoef (ratnumerator tt) k risch-var)
 812            yd (r* (ratdenominator tt)   ;DENOM MAY BE 0
 813                   (cond ((zerop risch-alphar)
 814                          (polcoef s risch-gamma risch-var))
 815                         (t
 816                          (funcall denom)))))
 817      (cond ((rzerop yd)
 818             (cond ((pzerop yn)
 819                    (setq k (1- k) risch-alphar (1- risch-alphar))
 820                    (go loop))           ;need more constraints?
 821                   (t
 822                    (cond
 823                      ((null flag)
 824                       (return nil))
 825                      (t
 826                       (return (cxerfarg (rzero) expg n a risch-ratform risch-intvar risch-mainvar)))))))
 827            (t
 828             (setq yalpha (ratqu yn yd))))
 829      (setq ytemp (r+ risch-y (r* yalpha
 830                                  (cons (list risch-mainvar risch-alphar 1) 1) )))
 831      (setq ttemp (r- tt (r* yalpha
 832                             (r+ (r* s (cons (list risch-mainvar risch-alphar 1) 1))
 833                                 (r* r risch-alphar
 834                                     (list risch-mainvar (1- risch-alphar) 1))))))
 835      (decf k)
 836      (decf risch-alphar)
 837      (cond ((< risch-alphar 0)
 838             (cond
 839               ((rzerop ttemp)
 840                (cond
 841                  ((null flag)
 842                   (return (ratqu ytemp p)))
 843                  (t
 844                   (return (list (ratqu (r* ytemp (cons (list expg n 1) 1))
 845                                        p)
 846                                 0)))))
 847               ((null flag)
 848                (return nil))
 849               ((and (risch-constp (setq ttemp (ratqu ttemp lcm)) risch-mainvar)
 850                     $erfflag
 851                     (equal (pdegree (car (get expg 'rischarg)) risch-mainvar) 2)
 852                     (equal (pdegree (cdr (get expg 'rischarg)) risch-mainvar) 0))
 853                (return (list (ratqu (r* ytemp (cons (list expg n 1) 1)) p)
 854                              (erfarg2 (r* n (get expg 'rischarg))
 855                                       ttemp risch-ratform risch-intvar risch-mainvar))))
 856               (t
 857                (return
 858                  (cxerfarg
 859                   (ratqu (r* risch-y (cons (list expg n 1) 1)) p)
 860                   expg
 861                   n
 862                   (ratqu tt lcm)
 863                   risch-ratform
 864                   risch-intvar
 865                   risch-mainvar))))))
 866      (setq risch-y ytemp)
 867      (setq tt ttemp)
 868      (go loop)))
 869
 870
 871 ;; *JM should be declared as an array, although it is not created
 872 ;; by this file. -- cwh
 873
 874 (defun lsa (mm)
 875   (prog (d *mosesflag m2 risch-m)
 876      (setq d (length (car mm)))
 877      ;; MTOA stands for MATRIX-TO-ARRAY.  An array is created and
 878      ;; associated functionally with the symbol *JM.  The elements
 879      ;; of the array are initialized from the matrix MM.
 880      (mtoa '*jm* (length mm) d mm)
 881      (setq risch-m (tfgeli '*jm*  (length mm) d))
 882      (cond ((or (and (null (car risch-m)) (null (cadr risch-m)))
 883                 (and (car risch-m)
 884                      (> (length (car risch-m)) (- (length mm) (1- d)))))
 885             (return (values 'singular risch-m)))
 886            ((cadr risch-m) (return (values 'inconsistent risch-m))))
 887      (setq *mosesflag t)
 888      (ptorat '*jm* (1- d) d)
 889      (setq m2 (xrutout '*jm* (1- d) d nil nil))
 890      (setq m2 (lsafix (cdr m2) (caddr risch-m)))
 891      (return (values m2 risch-m))))
 892
 893 (defun lsafix (l n)
 894   (declare (special *jm*))
 895   (do ((n n (cdr n))
 896        (l l (cdr l)))
 897       ((null l))
 898     (setf (aref *jm* 1 (car n)) (car l)))
 899   (do ((s (length l) (1- s))
 900        (ans))
 901       ((= s 0) (cons '(list) ans))
 902     (setq ans (cons (aref *jm* 1 s) ans))))
 903
 904
 905 (defun findpr (alist flist risch-y risch-mainvar &aux (p 1) fterm)
 906   (let (risch-alphar)
 907     (do ((alist alist (cdr alist))) ((null alist))
 908       (setq fterm (findflist (cadar alist) flist))
 909       (if fterm (setq flist (remove risch-y flist :count 1 :test #'eq)))
 910       (setq risch-alphar
 911             (cond ((null fterm) (caddar alist))
 912                   ((equal (caddr fterm) 1)
 913                    (fpr-dif (car flist) (caddar alist) risch-mainvar))
 914                   (t (max (- (caddar alist) (caddr fterm)) 0))))
 915       (if (not (zerop risch-alphar))
 916           (setq p (ptimes p (pexpt (cadar alist) risch-alphar)))))
 917     (do ((flist flist (cdr flist)))
 918         ((null flist))
 919       (when (equal (caddar flist) 1)
 920         (setq risch-alphar (fpr-dif (car flist) 0 risch-mainvar))
 921         (setq p (ptimes p (pexpt (cadar flist) risch-alphar)))))
 922     (values p risch-alphar)))
 923
 924 (defun fpr-dif (fterm alpha risch-mainvar)
 925   (destructuring-let* (((num den mult) fterm)
 926                        (risch-m (spderivative den risch-mainvar))
 927                        (n))
 928     (cond ((rzerop risch-m) alpha)
 929           (t (setq n (ratqu (cdr (ratdivide num den))
 930                             risch-m))
 931              (if (and (equal (cdr n) 1) (numberp (car n)))
 932                  (max (car n) alpha)
 933                  alpha)))))
 934
 935 (defun findflist (a llist)
 936   (cond ((null llist) nil)
 937         ((equal (cadar llist) a) (car llist))
 938         (t (findflist a (cdr llist)))))
 939
 940
 941 (defun rischexplog (expexpflag flag f a l
 942                     risch-ratform risch-intvar risch-liflag risch-degree risch-y risch-var risch-mainvar)
 943   (prog (lcm yy risch-m p risch-alphar risch-gamma delta
 944          mu r s tt denom ymu rbeta expg n eta logeta logdiff
 945          temp risch-cary risch-nogood vector aarray rmu rrmu rarray
 946          risch-beta
 947          risch-lians)
 948      (desetq (expg n eta logeta logdiff) l)
 949      (cond ((or (pzerop a)
 950                 (pzerop (car a)))
 951             (return (cond ((null flag)
 952                            (rzero))
 953                           (t
 954                            (rischzero))))))
 955      (multiple-value-setq (p risch-alphar)
 956        (findpr (cdr (partfrac a risch-var))
 957                (cdr (partfrac f risch-var))
 958                risch-y
 959                risch-mainvar))
 960      (setq lcm (plcm (ratdenominator a) p))
 961      (setq risch-y (ratpl (spderivative (cons 1 p) risch-mainvar)
 962                           (ratqu f p)))
 963      (setq lcm (plcm lcm (ratdenominator risch-y)))
 964      (setq r (car (ratqu lcm p)))
 965      (setq s (car (r* lcm risch-y)))
 966      (setq tt (car (r* a lcm)))
 967      (setq risch-beta (pdegree r risch-var))
 968      (setq risch-gamma (pdegree s risch-var))
 969      (setq delta (pdegree tt risch-var))
 970      (cond (expexpflag
 971             (setq mu (max (- delta risch-beta)
 972                           (- delta risch-gamma)))
 973             (go expcase)))
 974      (setq mu (max (- (1+ delta) risch-beta)
 975                    (- (1+ delta) risch-gamma)))
 976      (cond ((< risch-beta risch-gamma)
 977             (go back))
 978            ((= (1- risch-beta) risch-gamma)
 979             (go down1)))
 980      (setq risch-y (tryrisch1 (ratqu (r- (r* (polcoef r (1- risch-beta) risch-var)
 981                                              (polcoef s risch-gamma risch-var))
 982                                          (r* (polcoef r risch-beta risch-var)
 983                                              (polcoef s (1- risch-gamma) risch-var)))
 984                                      (r* (polcoef r risch-beta risch-var)
 985                                          (polcoef r risch-beta risch-var) ))
 986                               risch-mainvar risch-ratform risch-intvar risch-liflag risch-degree))
 987      (setq risch-cary (car risch-y))
 988      (multiple-value-setq (yy risch-cary risch-nogood risch-lians)
 989        (getfncoeff (cdr risch-y)
 990                    (get risch-var 'rischexpr)
 991                    risch-intvar risch-liflag risch-degree risch-cary risch-nogood
 992                    risch-lians))
 993      (cond ((and (not (findint (cdr risch-y)))
 994                  (not risch-nogood)
 995                  (not (atom yy))
 996                  (equal (cdr yy) 1)
 997                  (numberp (car yy))
 998                  (> (car yy) mu))
 999             (setq mu (car yy))))
1000      (go back)
1001    expcase
1002      (cond ((not (equal risch-beta risch-gamma))
1003             (go back)))
1004      (setq risch-y (tryrisch1 (ratqu (polcoef s risch-gamma risch-var) (polcoef r risch-beta risch-var))
1005                               risch-mainvar risch-ratform risch-intvar risch-liflag risch-degree))
1006      (cond ((findint (cdr risch-y))
1007             (go back)))
1008      (setq yy (ratqu (r* -1 (car risch-y)) eta))
1009      (cond ((and (not (atom yy))
1010                  (equal (cdr yy) 1)
1011                  (numberp (car yy))
1012                  (> (car yy) mu))
1013             (setq mu (car yy))))
1014      (go back)
1015    down1
1016      (setq risch-y (tryrisch1 (ratqu (polcoef s risch-gamma risch-var) (polcoef r risch-beta risch-var))
1017                               risch-mainvar risch-ratform risch-intvar risch-liflag risch-degree))
1018      (setq risch-cary (car risch-y))
1019      (multiple-value-setq (yy risch-cary risch-nogood risch-lians)
1020        (getfncoeff (cdr risch-y)
1021                    (get risch-var 'rischexpr)
1022                    risch-intvar risch-liflag risch-degree risch-cary risch-nogood
1023                    risch-lians))
1024      (cond ((and (not (findint (cdr risch-y)))
1025                  (not risch-nogood)
1026                  (not (atom yy))
1027                  (equal (cdr yy) 1)
1028                  (numberp (car yy))
1029                  (> (- (car yy)) mu))
1030             (setq mu (- (car yy)))))
1031    back
1032      (if (minusp mu)
1033          (return (if flag
1034                      (cxerfarg (rzero) expg n a risch-ratform risch-intvar risch-mainvar)
1035                      nil)))
1036      (cond ((> risch-beta risch-gamma)
1037             (go lsacall))
1038            ((= risch-beta risch-gamma)
1039             (go recurse)))
1040      (setq denom (polcoef s risch-gamma risch-var))
1041      (setq risch-y '(0 . 1))
1042    linearloop
1043      (setq ymu (ratqu (polcoef (ratnumerator tt) (+ mu risch-gamma) risch-var)
1044                       (r* (ratdenominator tt) denom)))
1045      (setq risch-y (r+ risch-y (setq ymu (r* ymu (pexpt (list logeta 1 1) mu) ))))
1046      (setq tt (r- tt
1047                   (r* s ymu)
1048                   (r* r (spderivative ymu risch-mainvar))))
1049      (decf mu)
1050      (cond  ((not (< mu 0))
1051              (go linearloop))
1052             ((not flag)
1053              (return (if (rzerop tt)
1054                          (ratqu risch-y p)
1055                          nil)))
1056             ((rzerop tt)
1057              (return (cons (ratqu (r* risch-y (cons (list expg n 1) 1)) p) '(0))))
1058             (t
1059              (return (cxerfarg (ratqu (r* risch-y (cons (list expg n 1) 1)) p)
1060                                expg
1061                                n
1062                                (ratqu tt lcm)
1063                                risch-ratform
1064                                risch-intvar
1065                                risch-mainvar))))
1066    recurse
1067      (setq rbeta (polcoef r risch-beta risch-var))
1068      (setq risch-y '(0 . 1))
1069    recurseloop
1070      (setq f (r+ (ratqu (polcoef s risch-gamma risch-var) rbeta)
1071                  (if expexpflag
1072                      (r* mu (spderivative eta risch-mainvar))
1073                      0)))
1074      (setq ymu (exppolycontrol nil
1075                                f
1076                                (ratqu (polcoef (ratnumerator tt)
1077                                                (+ risch-beta mu)
1078                                                risch-var)
1079                                       (r* (ratdenominator tt) rbeta))
1080                                expg
1081                                n
1082                                risch-ratform
1083                                risch-intvar
1084                                risch-liflag
1085                                risch-degree
1086                                risch-mainvar))
1087      (when (null ymu)
1088        (return (cond ((null flag)
1089                       nil)
1090                      (t
1091                       (return (cxerfarg (ratqu (r* risch-y (cons (list expg n 1) 1)) p)
1092                                         expg n (ratqu tt lcm)
1093                                         risch-ratform risch-intvar risch-mainvar))))))
1094      (setq risch-y (r+ risch-y (setq ymu (r* ymu (pexpt (list logeta 1 1) mu)))))
1095      (setq tt (r- tt
1096                   (r* s ymu)
1097                   (r* r (spderivative ymu risch-mainvar))))
1098      (decf mu)
1099      (cond
1100        ((not (< mu 0))
1101         (go recurseloop))
1102        ((not flag)
1103         (return (cond ((rzerop tt)
1104                        (ratqu risch-y p))
1105                       (t nil))))
1106        ((rzerop tt)
1107         (return (cons (ratqu (r* risch-y (cons (list expg n 1) 1)) p)
1108                       '(0))))
1109        (t
1110         (return (cxerfarg (ratqu (r* risch-y (cons (list expg n 1) 1)) p)
1111                           expg
1112                           n
1113                           (ratqu tt lcm)
1114                           risch-ratform
1115                           risch-intvar
1116                           risch-mainvar))))
1117    lsacall
1118      (setq rrmu mu)
1119    muloop
1120      (setq temp (r* (ratexpt (cons (list logeta 1 1) 1) (1- mu))
1121                     (r+ (r* s (cons (list logeta 1 1) 1))
1122                         (r* mu r logdiff ))))
1123    mu1
1124      (setq vector nil)
1125      (setq rmu (+ rrmu risch-beta))
1126    rmuloop
1127      (setq vector (cons (ratqu (polcoef (ratnumerator temp) rmu risch-var)
1128                                (ratdenominator temp)) vector))
1129      (decf rmu)
1130      (unless (< rmu 0)
1131        (go rmuloop))
1132      (decf mu)
1133      (setq aarray (append aarray (list (reverse vector))))
1134      (cond ((not (< mu 0))
1135             (go muloop))
1136            ((equal mu -2)
1137             (go skipmu)))
1138      (setq temp tt)
1139      (go mu1)
1140    skipmu
1141      (setq rarray nil)
1142    arrayloop
1143      (setq vector nil)
1144      (setq vector (mapcar 'car aarray))
1145      (setq aarray (mapcar 'cdr aarray))
1146      (setq rarray (append rarray (list vector)))
1147      (unless (null (car aarray)) (go arrayloop))
1148      (setq rmu (1+ rrmu))
1149      (setq vector nil)
1150    array1loop
1151      (setq vector (cons '(0 . 1) vector))
1152      (decf rmu)
1153      (unless (< rmu 0) (go array1loop))
1154      (setq aarray nil)
1155    array2loop
1156      (cond ((equal (car rarray) vector)
1157             nil)
1158            (t
1159             (setq aarray (cons (car rarray) aarray))))
1160      (setq rarray (cdr rarray))
1161      (when rarray (go array2loop))
1162      (setq rarray (reverse aarray))
1163      (multiple-value-setq (temp risch-m)
1164        (lsa rarray))
1165      (when (or (eq temp 'singular)
1166                (eq temp 'inconsistent))
1167        (return (if (null flag)
1168                    nil
1169                    (cxerfarg (rzero) expg n a risch-ratform risch-intvar risch-mainvar))))
1170      (setq temp (reverse  (cdr temp)))
1171      (setq rmu 0)
1172      (setq risch-y 0)
1173    l3
1174      (setq risch-y (r+ risch-y (r* (car temp) (pexpt (list logeta 1 1) rmu))))
1175      (setq temp (cdr temp))
1176      (incf rmu)
1177      (unless (> rmu rrmu)
1178        (go l3))
1179      (return (if (null flag)
1180                  (ratqu risch-y p)
1181                  (cons (r* (list expg n 1) (ratqu risch-y p)) '(0))))))
1182
1183
1184 (defun erfarg (exparg coef risch-ratform risch-mainvar)
1185   (prog (num denom erfarg)
1186      (setq exparg (r- exparg))
1187      (unless (and (setq num (pnthrootp (ratnumerator exparg) 2))
1188                   (setq denom (pnthrootp (ratdenominator exparg) 2)))
1189        (return nil))
1190      (setq erfarg (cons num denom))
1191      (if (risch-constp
1192           (setq coef (ratqu coef (spderivative erfarg risch-mainvar)))
1193           risch-mainvar)
1194          (return (simplify `((mtimes) ((rat) 1 2)
1195                              ((mexpt) $%pi ((rat) 1 2))
1196                              ,(disrep coef risch-ratform)
1197                              ((%erf) ,(disrep erfarg risch-ratform))))))))
1198
1199 (defun erfarg2 (exparg coeff risch-ratform risch-intvar risch-mainvar)
1200   (let ((risch-var risch-mainvar)
1201         a b c d)
1202     (when (and (= (pdegree (car exparg) risch-var) 2)
1203                (eq (caar exparg) risch-var)
1204                (risch-pconstp (cdr exparg) risch-mainvar)
1205                (risch-constp coeff risch-mainvar))
1206       (setq a (ratqu (r* -1 (caddar exparg))
1207                      (cdr exparg)))
1208       (setq b (disrep (ratqu (r* -1 (polcoef (car exparg) 1 risch-var))
1209                              (cdr exparg))
1210                       risch-ratform))
1211       (setq c (disrep (ratqu (r* (polcoef (car exparg) 0 risch-var))
1212                              (cdr exparg))
1213                       risch-ratform))
1214       (setq d (ratsqrt a risch-ratform))
1215       (setq a (disrep a risch-ratform))
1216       (simplify `((mtimes)
1217                   ((mtimes)
1218                    ((mexpt) $%e ((mplus) ,c
1219                                  ((mquotient) ((mexpt) ,b 2)
1220                                   ((mtimes) 4 ,a))))
1221                    ((rat) 1 2)
1222                    ,(disrep coeff risch-ratform)
1223                    ((mexpt) ,d -1)
1224                    ((mexpt) $%pi ((rat) 1 2)))
1225                   ((%erf) ((mplus)
1226                            ((mtimes) ,d ,risch-intvar)
1227                            ((mtimes) ,b ((rat) 1 2) ((mexpt) ,d -1)))))))))
1228
1229
1230 (defun cxerfarg (ans expg n numdenom risch-ratform risch-intvar risch-mainvar
1231                  &aux (arg (r* n (get expg 'rischarg)))
1232                    (fails 0))
1233   (prog (denom erfans num nerf)
1234      (desetq (num . denom) numdenom)
1235      (unless $erfflag
1236        (setq fails num)
1237        (go lose))
1238      (if (setq erfans (erfarg arg numdenom risch-ratform risch-mainvar))
1239          (return (list ans erfans)))
1240    again
1241      (when (and (not (pcoefp denom))
1242                 (null (p-red denom))
1243                 (eq (get (car denom) 'leadop) 'mexpt))
1244        (setq arg (r+ arg (r* (- (p-le denom))
1245                              (get (p-var denom) 'rischarg)))
1246              denom (p-lc denom))
1247        (go again))
1248      (loop for (coef exparg exppoly) in (explist num arg 1)
1249            do (setq coef (ratqu coef denom)
1250                     nerf (or (erfarg2 exparg coef risch-ratform risch-intvar risch-mainvar)
1251                              (erfarg exparg coef risch-ratform risch-mainvar)))
1252               (if nerf
1253                   (push nerf erfans)
1254                   (setq fails (pplus fails exppoly))))
1255    lose
1256      (return
1257        (if (pzerop fails)
1258            (cons ans erfans)
1259            (rischadd (cons ans erfans)
1260                      (rischnoun (r* (ratexpt (cons (make-poly expg) 1) n)
1261                                     (ratqu fails (cdr numdenom)))
1262                                 risch-ratform
1263                                 risch-intvar))))))
1264
1265 (defun explist (p oarg exps)
1266   (cond ((or (pcoefp p)
1267              (not (eq 'mexpt
1268                       (get (p-var p) 'leadop))))
1269          (list (list p oarg (ptimes p exps))))
1270         (t
1271          (loop with narg = (get (p-var p) 'rischarg)
1272                for (exp coef) on (p-terms p) by #'cddr
1273                nconc (explist coef
1274                               (r+ oarg (r* exp narg))
1275                               (ptimes exps
1276                                       (make-poly (p-var p) exp 1)))))))
1277
1278
1279 (defun intsetup (exp risch-*var)
1280   (prog (varlist clist $factorflag dlist genpairs old risch-y z $ratfac $keepfloat
1281          *fnewvarsw)
1282    y
1283      (setq exp (radcan1 exp risch-*var))
1284      (fnewvar exp)
1285      (setq *fnewvarsw t)
1286    a
1287      (setq clist nil)
1288      (setq dlist nil)
1289      (setq z varlist)
1290    up
1291      (setq risch-y (pop z))
1292      (cond ((freeof risch-*var risch-y)
1293             (push risch-y clist))
1294            ((eq risch-y risch-*var)
1295             nil)
1296            ((and (mexptp risch-y)
1297                  (not (eq (cadr risch-y) '$%e)))
1298             (cond ((not (freeof risch-*var (caddr risch-y)))
1299                    (setq dlist `((mexpt simp)
1300                                  $%e
1301                                  ,(mul2* (caddr risch-y)
1302                                          `((%log) ,(cadr risch-y)))))
1303                    (setq exp (maxima-substitute dlist risch-y exp))
1304                    (setq varlist nil)
1305                    (go y))
1306                   ((atom (caddr risch-y))
1307                    (cond ((numberp (caddr risch-y))
1308                           (push risch-y dlist))
1309                          (t
1310                           ;;(setq operator t)
1311                           (return (values nil t)))))
1312                   (t
1313                    (push risch-y dlist))))
1314            (t
1315             (push risch-y dlist)))
1316      (if z
1317          (go up))
1318      (if (member '$%i clist :test #'eq)
1319          (setq clist (cons '$%i (delete '$%i clist :test #'equal))))
1320      (setq varlist (append clist
1321                            (cons risch-*var
1322                                  (nreverse (sort (append dlist nil)
1323                                                  #'(lambda (a b)
1324                                                      (intgreat a b risch-*var)))))))
1325      (orderpointer varlist)
1326      (setq old varlist)
1327      (mapc #'(lambda (b)
1328                (intset1 b risch-*var))
1329            (cons risch-*var dlist))
1330      (cond ((alike old varlist)
1331             (return (values (ratrep* exp)
1332                             nil)))
1333            (t (go a)))))
1334
1335 (defun leadop (exp)
1336   (cond ((atom exp)
1337          exp)
1338         ((mqapplyp exp)
1339          (cadr exp))
1340         (t
1341          (caar exp))))
1342
1343 (defun leadarg (exp)
1344   (cond ((atom exp)
1345          0)
1346         ((and (mexptp exp) (eq (cadr exp) '$%e))
1347          (caddr exp))
1348         ((mqapplyp exp)
1349          (car (subfunargs exp)))
1350         (t
1351          (cadr exp))))
1352
1353 (defun intset1 (b risch-*var)
1354   (let (e c d)
1355     (fnewvar
1356      (setq d (if (mexptp b)             ;needed for radicals
1357                  `((mtimes simp)
1358                    ,b
1359                    ,(radcan1 (sdiff (simplify (caddr b)) risch-*var)
1360                              risch-*var))
1361                  (radcan1 (sdiff (simplify b) risch-*var)
1362                           risch-*var))))
1363     (setq d (ratrep* d))
1364     (setq c (ratrep* (leadarg b)))
1365     (setq e (cdr (assoc b (pair varlist genvar) :test #'equal)))
1366     (putprop e (leadop b) 'leadop)
1367     (putprop e b 'rischexpr)
1368     (putprop e (cdr d) 'rischdiff)
1369     (putprop e (cdr c) 'rischarg)))
1370
1371 ;; order of expressions for risch.
1372 ;; expressions containing erf and li last.
1373 ;; then order by size of expression to guarantee that
1374 ;; any subexpressions are considered smaller.
1375 ;; this relation should be transitive, since it is called by sort.
1376 (defun intgreat (a b risch-*var)
1377   (cond ((and (not (atom a))
1378               (not (atom b)))
1379          (cond ((and (not (freeof '%erf a))
1380                      (freeof '%erf b))
1381                 t)
1382                ((and (not (freeof '$li a))
1383                      (freeof '$li b))
1384                 t)
1385                ((and (freeof '$li a)
1386                      (not (freeof '$li b)))
1387                 nil)
1388                ((and (freeof '%erf a)
1389                      (not (freeof '%erf b)))
1390                 nil)
1391                ((> (conssize a) (conssize b))
1392                 t)
1393                ((< (conssize a) (conssize b))
1394                 nil)
1395                (t
1396                 (great (resimplify (fixintgreat a risch-*var))
1397                        (resimplify (fixintgreat b risch-*var))))))
1398         (t
1399          (great (resimplify (fixintgreat a risch-*var))
1400                 (resimplify (fixintgreat b risch-*var))))))
1401
1402 (defun fixintgreat (a risch-*var)
1403   (subst '/_101x risch-*var a))