src/trgred.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 1981 Massachusetts Institute of Technology         ;;;
   9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10
  11 (in-package :maxima)
  12
  13 (macsyma-module trgred)
  14
  15 (declare-top (special var *schatc-ans* *trigred *noexpand *lin))
  16
  17 (load-macsyma-macros rzmac)
  18
  19 (defvar *trans-list-plus*
  20   '((((mplus) ((coeffpt) (c true) ((mexpt) ((%tan) (x true)) 2))
  21        (var* (uvar) c))
  22       ((mtimes) c ((mexpt) ((%sec) x) 2)))
  23     (((mplus) ((coeffpt) (c true) ((mexpt) ((%cot) (x true)) 2))
  24       (var* (uvar) c))
  25      ((mtimes) c ((mexpt) ((%csc) x) 2)))
  26     (((mplus) ((coeffpt) (c true) ((mexpt) ((%tanh) (x true)) 2))
  27       ((mtimes) -1 (var* (uvar) c)))
  28      ((mtimes) -1 c ((mexpt) ((%sech) x) 2)))
  29     (((mplus) ((coeffpt) (c true) ((mexpt) ((%coth) (x true)) 2))
  30       ((mtimes) -1 (var* (uvar) c)))
  31      ((mtimes) c ((mexpt) ((%csch) x) 2)))))
  32
  33 (defvar *triglaws*
  34       '(* %sin (* %cot %cos %sec %tan) %cos (* %tan %sin %csc %cot)
  35           %tan (* %cos %sin %csc %sec) %cot (* %sin %cos %sec %csc)
  36           %sec (* %sin %tan %cot %csc) %csc (* %cos %cot %tan %sec)))
  37
  38 (defvar *hyperlaws*
  39       '(* %sinh (* %coth %cosh %sech %tanh) %cosh (* %tanh %sinh %csch %coth)
  40           %tanh (* %cosh %sinh %csch %sech) %coth (* %sinh %cosh %sech %csch)
  41           %sech (* %sinh %tanh %coth %csch) %csch (* %cosh %coth %tanh %sech)))
  42
  43 (defvar *sc^ndisp* '((%sin . sin^n) (%cos . cos^n) (%sinh . sinh^n) (%cosh . cosh^n)))
  44
  45 (defvar *laws*)
  46 (defvar *trigbuckets*)
  47 (defvar *hyperbuckets*)
  48
  49 ;;The Trigreduce file contains a group of routines which can be used to
  50 ;;make trigonometric simplifications of expressions.  The bulk of the
  51 ;;routines here involve the reductions of products of sin's and cos's.
  52 ;;
  53 ;;      *TRIGRED        indicates that the special simplifications for
  54 ;;                      $TRIGREDUCE are to be used.
  55 ;;      *NOEXPAND       indicates that trig functions of sums of
  56 ;;                      angles are not to be used.
  57
  58 ;; Without the check for the value of $simp, trigreduce can enter an
  59 ;; endless loop--for example trigreduce((-1)*x). There is at least one
  60 ;; test in the testsuite that sets simp to false and then needs to compute
  61 ;; a limit. But the limit code calls trigreduce and the test fails. The
  62 ;; top level check for simp seems like the best way to prevent infinite
  63 ;; loops when simp is false.
  64 (defmfun ($trigreduce :properties ((evfun t))) (exp &optional (var '*novar))
  65   (let ((*trigred t)
  66         (*noexpand t)
  67         $trigexpand $verbose $ratprint)
  68     (declare (special $trigexpand *trigred))
  69     (if $simp (gcdred (sp1 exp)) exp)))
  70
  71 ;; The first pass in power series expansion (used by $powerseries in
  72 ;; series.lisp), but also used by $trigreduce. Expands / reduces the expression
  73 ;; E as a function VAR, controlled by the *TRIGRED and *NOEXPAND flags.
  74 (defun sp1 (e)
  75   (cond
  76     ((atom e) e)
  77
  78     ;; We recognise some special patterns that are expressed as sums, such as
  79     ;; 1+tan(x)^2 => sec(x)^2. Rewrite using the first matching pattern (and
  80     ;; recurse to try to simplify further). If no pattern matches, expand each
  81     ;; term in the sum.
  82     ((eq (caar e) 'mplus)
  83      (or (dolist (pair *trans-list-plus*)
  84            (let ((match (m2 e (first pair))))
  85              (when match
  86                (return (sp1 (sch-replace match (second pair)))))))
  87          (m+l (mapcar #'sp1 (cdr e)))))
  88
  89     ((eq (caar e) 'mtimes)
  90      (sp1times e))
  91
  92     ((eq (caar e) 'mexpt)
  93      (sp1expt (sp1 (cadr e)) (sp1 (caddr e))))
  94
  95     ((eq (caar e) '%log)
  96      (sp1log (sp1 (cadr e))))
  97
  98     ((member (caar e) '(%cos %sin %tan %cot %sec %csc
  99                         %cosh %sinh %tanh %coth %sech %csch) :test #'eq)
 100      (sp1trig (list (car e) (let* ((*noexpand t)) (sp1 (cadr e))))))
 101
 102     ((member (caar e) '(%asin %acos %atan %acot %asec %acsc
 103                         %asinh %acosh %atanh %acoth %asech %acsch) :test #'eq)
 104      (sp1atrig (caar e) (let* ((*noexpand t)) (sp1 (cadr e)))))
 105
 106     ((eq (caar e) 'mrat) (sp1 (ratdisrep e)))
 107
 108     ((mbagp e)
 109      (cons (list (caar e)) (mapcar #'(lambda (u)
 110                                        (gcdred (sp1 u)))
 111                                    (cdr e))))
 112     ((eq (caar e) '%integrate)
 113      (list* '(%integrate) (sp1 (cadr e)) (cddr e)))
 114
 115     (t (recur-apply #'sp1 e))))
 116
 117 (defun trigfp (e)
 118   (or (and (not (atom e)) (trigp (caar e))) (equal e 1)))
 119
 120 (defun gcdred (e)
 121   (cond ((atom e) e)
 122         ((eq (caar e) 'mplus) (m+l (mapcar #'gcdred (cdr e))))
 123         ((eq (caar e) 'mtimes)
 124          (let* ((nn '(1))
 125                 (nd '(1))
 126                 (gcd nil))
 127            (do ((e (cdr e) (cdr e)))
 128                ((null e)
 129                 (setq nn (m*l nn) nd (m*l nd)))
 130              (cond ((and (mexptp (car e))
 131                          (or (signp l (caddar e))
 132                              (and (mtimesp (caddar e))
 133                                   (signp l (cadr (caddar e))))))
 134                     (setq nd (cons (m^ (cadar e) (m- (caddar e))) nd)))
 135                    ((ratnump (car e))
 136                     (setq nn (cons (cadar e) nn)
 137                           nd (cons (caddar e) nd)))
 138                    ((setq nn (cons (car e) nn)))))
 139            (cond ((equal nd 1) nn)
 140                  ((equal (setq gcd ($gcd nn nd)) 1) e)
 141                  ((div* (cadr ($divide nn gcd))
 142                         (cadr ($divide nd gcd)))))))
 143         (t (recur-apply #'gcdred e))))
 144
 145 (defun sp1times (e)
 146   (let* ((fr nil)
 147          (g '(1))
 148          (*trigbuckets* nil)
 149          (*hyperbuckets* nil)
 150          (tr nil)
 151          (hyp nil)
 152          (*lin '(0)))
 153     (dolist (factor (cdr e))
 154       (cond ((or (mnump factor)
 155                  (and (not (eq var '*novar)) (free factor var)))
 156              (push factor fr))
 157             ((atom factor) (push factor g))
 158             ((or (trigfp factor)
 159                  (and (eq (caar factor) 'mexpt)
 160                       (trigfp (cadr factor))))
 161              (sp1add factor))
 162             (t
 163              (push factor g))))
 164     (setq g (mapcar #'sp1 g))
 165
 166     (mapc #'(lambda (q)  (sp1sincos q t)) *trigbuckets*)
 167     (mapc #'(lambda (q) (sp1sincos q nil)) *hyperbuckets*)
 168     (setq fr (cons (m^ 1//2 (m+l *lin)) fr)
 169           *lin nil)
 170     (setq tr (cons '* (mapcan #'sp1untrep *trigbuckets*)))
 171     (setq g (nconc (sp1tlin tr t) (sp1tplus *lin t) g)
 172           *lin nil)
 173     (setq hyp (cons '* (mapcan #'sp1untrep *hyperbuckets*)))
 174     (setq g (nconc (sp1tlin hyp nil) (sp1tplus *lin nil) g))
 175     (setq g ($expand (let* (($keepfloat t))
 176                        (declare (special $keepfloat))
 177                        ($ratsimp (cons '(mtimes) g)))))
 178     (if (mtimesp g)
 179         (setq g (mapcar #'sp1 (cdr g)))
 180         (setq g (list (sp1 g))))
 181     (m*l (cons 1 (nconc g fr (cdr tr) (cdr hyp))))))
 182
 183 (defun sp1tlin (l trig)
 184   (cond ((null (cdr l)) nil)
 185         ((and (eq (caaadr l) 'mexpt)
 186               (integerp (caddr (cadr l)))
 187               (member (caaadr (cadr l))
 188                     (if trig '(%sin %cos) '(%sinh %cosh)) :test #'eq))
 189          (cons (funcall (cdr (assoc (caaadr (cadr l)) *sc^ndisp* :test #'eq))
 190                         (caddr (cadr l)) (cadadr (cadr l)))
 191                (sp1tlin (rplacd l (cddr l)) trig)))
 192         ((member (caaadr l) (if trig '(%sin %cos) '(%sinh %cosh)) :test #'eq)
 193          (push (cadr l) *lin)
 194          (sp1tlin (rplacd l (cddr l)) trig))
 195         ((sp1tlin (cdr l) trig))))
 196
 197 ;; Rewrite a product of sines and cosines as a sum
 198 ;;
 199 ;; L is a list of sines and cosines. For example, consider the list
 200 ;;
 201 ;;  sin(x), sin(2*x), sin(3*x)
 202 ;;
 203 ;; This represents the product sin(x)*sin(2*x)*sin(3*x).
 204 ;;
 205 ;; ANS starts as sin(x). Then for each term in the rest of the list, we multiply
 206 ;; the answer that we have found so far by that term. The result will be a sum
 207 ;; of sines. In this example, sin(x)*sin(2*x) gives us
 208 ;;
 209 ;;  1/2 * (cos(x) - cos(3*x))
 210 ;;
 211 ;; In fact we don't calculate the 1/2 coefficient in sp1sintcos: you always get
 212 ;; a factor of 2^(k-1), where k is the length of the list, so this is calculated
 213 ;; at the bottom of sp1tplus. Anyway, next we calculate cos(x)*sin(3*x) and
 214 ;; -cos(3*x)*sin(3*x) and sum the answers. Note that -cos(3*x) will crop up
 215 ;; represented as ((mtimes) -1 ((%cos) ((mtimes) 3 $x))). See note in the let
 216 ;; form for info on what form ANS must take.
 217 (defun sp1tplus (l trig)
 218   (if (or (null l) (null (cdr l)))
 219       l
 220       ;; ANS is a list containing the terms in a sum for the expanded
 221       ;; expression. Each element in this list is either of the form sc(x),
 222       ;; where sc is sin or cos, or of the form ((mtimes) coeff ((sc) $x)),
 223       ;; where coeff is some coefficient.
 224       ;;
 225       ;; multiply-sc-terms rewrites a*sc as a sum of sines and cosines. The
 226       ;; result is a list containing the terms in a sum which is
 227       ;; mathematically equal to a*sc. Assuming that term is of one of the
 228       ;; forms described for ANS below and that SC is of the form sc(x), the
 229       ;; elements of the resulting list will all be of suitable form for
 230       ;; inclusion into ANS.
 231       (flet ((multiply-sc-terms (term sc)
 232                (let* ((coefficient (when (mtimesp term) (cadr term)))
 233                       (term-sc (if (mtimesp term) (caddr term) term))
 234                       (expanded (sp1sintcos term-sc sc trig)))
 235                  ;; expanded will now either be sin(foo) or cos(foo) OR it
 236                  ;; will be a sum of such terms.
 237                  (cond
 238                    ((not coefficient) (list expanded))
 239                    ((or (atom expanded)
 240                         (member (caar expanded) '(%sin %cos %sinh %cosh)
 241                                 :test 'eq))
 242                     (list (m* coefficient expanded)))
 243                    ((mplusp expanded)
 244                     (mapcar (lambda (summand) (m* coefficient summand))
 245                             (cdr expanded)))
 246                    (t
 247                     ;; SP1SINTCOS can also return numbers and constant expressions.
 248                     ;; Assume that's the case here.
 249                     (list (m* coefficient expanded))))))
 250              ;; Treat EXPR as a sum and return a list of its terms
 251              (terms-of-sum (expr)
 252                (if (mplusp expr) (cdr expr) (ncons expr))))
 253
 254         (let ((ans (list (first l))))
 255           (dolist (sc (rest l))
 256             (setq ans (terms-of-sum
 257                        (m+l (mapcan (lambda (q)
 258                                       (multiply-sc-terms q sc)) ans)))))
 259           (list (list '(rat) 1 (expt 2 (1- (length l))))
 260                 (m+l ans))))))
 261
 262 ;; The core of trigreduce. Performs transformations like sin(x)*cos(x) =>
 263 ;; sin(2*x)
 264 ;;
 265 ;; This function only does something non-trivial if both a and b have one of
 266 ;; sin, cos, sinh and cosh as top-level operators. (Note the first term in the
 267 ;; cond: we assume that if a,b are non-atomic and not both of them are
 268 ;; hyperbolic/trigonometric then we can just multiply the two terms)
 269 (defun sp1sintcos (a b trig)
 270   (let* ((x nil)
 271          (y nil))
 272     (cond ((or (atom a) (atom b)
 273                (not (member (caar a) '(%sin %cos %sinh %cosh) :test #'eq))
 274                (not (member (caar b) '(%sin %cos %sinh %cosh) :test #'eq)))
 275            (mul3 2 a b))
 276           ((prog2 (setq x (m+ (cadr a) (cadr b)) y (m- (cadr a) (cadr b)))
 277                (null (eq (caar a) (caar b))))
 278            (setq b (if trig '(%sin) '(%sinh)))
 279            (or (eq (caar a) '%sin) (eq (caar a) '%sinh)
 280                (setq y (m- y)))
 281            (m+ (list b x) (list b y)))
 282           ((member (caar a) '(%cos %cosh) :test #'eq)
 283            (m+ (list (list (caar a)) x)
 284                (list (list (caar a)) y)))
 285           (trig
 286            (m- (list '(%cos) y) (list '(%cos) x)))
 287           ((m- (list '(%cosh) x) (list '(%cosh) y))))))
 288
 289 ;; For COS(X)^2, TRIGBUCKET is (X (1 (COS . 2))) or, more generally,
 290 ;; (arg (numfactor-of-arg (operator . exponent)))
 291
 292 (defun sp1add (e)
 293   (let* ((n (cond ((eq (caar e) 'mexpt)
 294                    (cond ((= (signum1 (caddr e)) -1)
 295                           (prog1 (m- (caddr e))
 296                             (setq e (cons (list (zl-get (caaadr e) 'recip)) (cdadr e)))))
 297                          ((prog1 (caddr e) (setq e (cadr e))))))
 298                   ( 1 )))
 299          (arg (sp1kget (cadr e)))
 300          (buc nil)
 301          (*laws* *hyperlaws*))
 302     (cond ((member (caar e) '(%sin %cos %tan %cot %sec %csc) :test #'eq)
 303            (cond ((setq buc (assoc (cdr arg) *trigbuckets* :test #'equal))
 304                   (setq *laws* *triglaws*)
 305                   (sp1addbuc (caar e) (car arg) n buc))
 306                  ((setq *trigbuckets*
 307                         (cons (list (cdr arg) (list (car arg) (cons (caar e) n)))
 308                               *trigbuckets*)))))
 309           ((setq buc (assoc (cdr arg) *hyperbuckets* :test #'equal))
 310            (sp1addbuc (caar e) (car arg) n buc))
 311           ((setq *hyperbuckets*
 312                  (cons (list (cdr arg) (list (car arg) (cons (caar e) n)))
 313                        *hyperbuckets*))))))
 314
 315 (defun sp1addbuc (f arg n b) ;FUNCTION, ARGUMENT, EXPONENT, BUCKET LIST
 316   (cond ((and (cdr b) (alike1 arg (caadr b))) ;GOES IN THIS BUCKET
 317          (sp1putbuc f n (cadr b)))
 318         ((or (null (cdr b)) (great (caadr b) arg))
 319          (rplacd b (cons (list arg (cons f n)) (cdr b))))
 320         ((sp1addbuc f arg n (cdr b)))))
 321
 322 (defun sp1putbuc (f n *buc)             ;PUT IT IN THERE
 323   (do ((buc *buc (cdr buc)))
 324       ((null (cdr buc))
 325        (rplacd buc (list (cons f n))))
 326     (cond ((eq f (caadr buc))           ;SAME FUNCTION
 327            (return
 328              (rplacd (cadr buc) (m+ n (cdadr buc))))) ;SO BOOST EXPONENT
 329           ((eq (caadr buc) (zl-get f 'recip)) ;RECIPROCAL FUNCTIONS
 330            (setq n (m- (cdadr buc) n))
 331            (return
 332              (cond ((signp e n) (rplacd buc (cddr buc)))
 333                    ((= (signum1 n) -1)
 334                     (rplaca (cadr buc) f)
 335                     (rplacd (cadr buc) (neg n)))
 336                    (t (rplacd (cadr buc) n)))))
 337           (t (let* ((nf (zl-get (zl-get *laws* (caadr buc)) f))
 338                     (m nil))
 339                (cond ((null nf))        ;NO SIMPLIFICATIONS HERE
 340                      ((equal n (cdadr buc)) ;EXPONENTS MATCH
 341                       (rplacd buc (cddr buc))
 342                       (return
 343                         (sp1putbuc1 nf n *buc))) ;TO MAKE SURE IT DOESN'T OCCUR TWICE
 344                      ((eq (setq m (sp1great n (cdadr buc))) 'nomatch))
 345                      (m (setq m (cdadr buc))
 346                         (rplacd buc (cddr buc))
 347                         (sp1putbuc1 nf m *buc)
 348                         (sp1putbuc1 f (m- n m) *buc)
 349                         (return t))
 350                      (t (rplacd (cadr buc) (m- (cdadr buc) n))
 351                         (return (sp1putbuc1 nf n *buc)))))))))
 352
 353 (defun sp1putbuc1 (f n buc)
 354   (cond ((null (cdr buc))
 355          (rplacd buc (list (cons f n))))
 356         ((eq f (caadr buc))
 357          (rplacd (cadr buc) (m+ n (cdadr buc))))
 358         ((sp1putbuc1 f n (cdr buc)))))
 359
 360 (defun sp1great (x y)
 361   (let* ((a nil)
 362          (b nil))
 363     (cond ((mnump x)
 364            (cond ((mnump y) (great x y)) (t 'nomatch)))
 365           ((or (atom x) (atom y)) 'nomatch)
 366           ((and (eq (caar x) (caar y))
 367                 (alike (cond ((mnump (cadr x))
 368                               (setq a (cadr x)) (cddr x))
 369                              (t (setq a 1) (cdr x)))
 370                        (cond ((mnump (cadr y))
 371                               (setq b (cadr y)) (cddr y))
 372                              (t (setq b 1) (cdr y)))))
 373            (great a b))
 374           (t 'nomatch))))
 375
 376 (defun sp1untrep (b)
 377   (mapcan
 378    #'(lambda (buc)
 379        (mapcar #'(lambda (term)
 380                    (let* ((bas (simplifya (list (list (car term))
 381                                                 (m* (car b) (car buc)))
 382                                           t)))
 383                      (cond ((equal (cdr term) 1) bas)
 384                            ((m^ bas (cdr term))))))
 385                (cdr buc)))
 386    (cdr b)))
 387
 388 (defun sp1kget (e)                      ;FINDS NUMERIC COEFFICIENTS
 389   (or (and (mtimesp e) (numberp (cadr e))
 390            (cons (cadr e) (m*l (cddr e))))
 391       (cons 1 e)))
 392
 393 (defun sp1sincos (l trig)
 394   (mapcar #'(lambda (q) (sp1sincos2 (m* (car l) (car q)) q trig)) (cdr l)))
 395
 396 (defun sp1sincos2 (arg l trig)
 397   (let* ((a nil))
 398     (cond ((null (cdr l)))
 399           ((and
 400             (setq a (member (caadr l)
 401                             (if (null trig)
 402                                 '(%sinh %cosh %sinh %csch %sech %csch)
 403                                 '(%sin %cos %sin %csc %sec %csc)) :test #'eq))
 404             (cddr l))            ;THERE MUST BE SOMETHING TO MATCH TO.
 405            (sp1sincos1 (cadr a) l arg trig))
 406           ((sp1sincos2 arg (cdr l) trig)))))
 407
 408 (defun sp1sincos1 (s l arg trig)
 409   (let* ((g nil)
 410          (e 1))
 411     (do ((ll (cdr l) (cdr ll)))
 412         ((null (cdr ll)) t)
 413       (cond ((eq s (caadr ll))
 414              (setq arg (m* 2 arg))
 415              (cond (trig
 416                     (cond ((member s '(%sin %cos) :test #'eq)
 417                            (setq s '%sin))
 418                           ((setq s '%csc e -1))))
 419                    (t
 420                     (cond ((member s '(%sinh %cosh) :test #'eq)
 421                            (setq s '%sinh))
 422                           ((setq s '%csch e -1)))))
 423              (cond ((alike1 (cdadr ll) (cdadr l))
 424                     (sp1addto s arg (cdadr l))
 425                     (setq *lin (cons (m* e (cdadr l)) *lin))
 426                     (rplacd ll (cddr ll))   ;;;MUST BE IN THIS ORDER!!
 427                     (rplacd l (cddr l))
 428                     (return t))
 429                    ((eq (setq g (sp1great (cdadr l) (cdadr ll))) 'nomatch))
 430                    ((null g)
 431                     (rplacd (cadr ll) (m- (cdadr ll) (cdadr l)))
 432                     (sp1addto s arg (cdadr l))
 433                     (setq *lin (cons (m* e (cdadr l)) *lin))
 434                     (rplacd l (cddr l))
 435                     (return t))
 436                    (t
 437                     (rplacd (cadr l) (m- (cdadr l) (cdadr ll)))
 438                     (sp1addto s arg (cdadr ll))
 439                     (push (m* e (cdadr ll)) *lin)
 440                     (rplacd ll (cddr ll))
 441                     (return t))))))))
 442
 443 (defun sp1addto (fn arg exp)
 444   (setq arg (list (list fn) arg))
 445   (sp1add (if (equal exp 1) arg (m^ arg exp))))
 446
 447 (defun sp1expt (b e)
 448   (cond ((mexptp b)
 449          (power (sp1 b) (sp1 e)))
 450         ((and (null (trigfp b)) (free e var))
 451          (m^ b e))
 452         ((equal b '$%e)
 453          (sp1expt2 e))
 454         ((and (null (eq var '*novar)) (free b var))
 455          (sp1expt2 (m* (list '(%log) b) e)))
 456         ((and (consp b) (consp (car b)) (member (caar b) '(%sin %cos %tan %cot %sec %csc
 457                                                                 %sinh %cosh %tanh %coth %sech %csch) :test #'eq))
 458          (cond ((= (signum1 e) -1)
 459                 (sp1expt (list (list (zl-get (caar b) 'recip)) (cadr b))
 460                          (neg e)))
 461                ((and (signp g e)
 462                      (member (caar b) '(%sin %cos %sinh %cosh) :test #'eq))
 463                 (funcall (cdr (assoc (caar b) *sc^ndisp* :test #'eq)) e (cadr b)))
 464                ((m^ b e))))
 465         ((m^ b e))))
 466
 467 (defun sp1expt2 (e)
 468   (let* ((*schatc-ans* (m2 e '((mplus) ((coeffpp) (fr freevar)) ((coeffpp) (exp true)))))
 469          (fr (cdr (assoc 'fr *schatc-ans* :test #'eq)))
 470          (exp (cdr (assoc 'exp *schatc-ans* :test #'eq))))
 471     (cond ((equal fr 0)
 472            (m^ '$%e exp))
 473           ((m* (m^ '$%e fr) (m^ '$%e exp))))))
 474
 475 ;; Split TERMS into (VALUES NON-NEG OTHER) where NON-NEG and OTHER are a
 476 ;; partition of the elements of TERMS. Expressions that are known not to be
 477 ;; negative are placed in NON-NEG and all others end up in OTHER.
 478 ;;
 479 ;; This function is used to safely split products when expanding logarithms to
 480 ;; avoid accidentally ending up with something like
 481 ;;
 482 ;;   log(1 - x) => log(-1) + log(x-1).
 483 ;;
 484 ;; Note that we don't check a term is strictly positive: if it was actually
 485 ;; zero, the logarithm was bogus in the first place.
 486 (defun non-negative-split (terms)
 487   (let ((non-neg) (other))
 488     (dolist (term terms)
 489       (if (memq ($sign term) '($pos $pz $zero))
 490           (push term non-neg)
 491           (push term other)))
 492     (values non-neg other)))
 493
 494 ;; Try to expand a logarithm for use in a power series in VAR by splitting up
 495 ;; products.
 496 (defun sp1log (e &optional no-recurse)
 497   (cond
 498     ;; If E is free of VAR, is an atom, or we're supposed to be reducing rather
 499     ;; than expanding, then just return E.
 500     ((or *trigred (atom e) (free e var))
 501      (list '(%log) e))
 502
 503     ;; The logarithm of a sum doesn't simplify very nicely, but call $factor to
 504     ;; see if we can pull out one or more terms and then recurse (setting
 505     ;; NO-RECURSE to make sure we don't end up in a loop)
 506     ((eq (caar e) 'mplus)
 507      (let* ((exp (m1- e)) *a *n)
 508        (declare (special *n *a))
 509        (cond
 510          ((smono exp var)
 511           (list '(%log) e))
 512          ((not no-recurse)
 513           (sp1log ($factor e) t))
 514          (t (sp1log2 e)))))
 515
 516     ;; A product is much more promising. Do the transformation log(ab) =>
 517     ;; log(a)+log(b) and pass it to SP1 for further simplification.
 518     ;;
 519     ;; We need to be a little careful here because eg. factor(1-x) gives
 520     ;; -(x-1). We don't want to end up with a log(-1) term! So check the sign of
 521     ;; terms and only pull out the terms we know to be non-negative. If the
 522     ;; argument was a negative real in the first place then we'd already got
 523     ;; rubbish, but otherwise we won't pull out anything we don't want.
 524     ((eq (caar e) 'mtimes)
 525      (multiple-value-bind (non-neg other) (non-negative-split (cdr e))
 526        (cond
 527          ((null non-neg) (sp1log2 e))
 528          (t
 529           (sp1 (m+l (mapcar #'sp1log (append other non-neg))))))))
 530
 531     ;; Similarly, transform log(a^b) => b log(a) and pass back to SP1.
 532     ((eq (caar e) 'mexpt)
 533      (sp1 (m* (caddr e) (list '(%log) (cadr e)))))
 534
 535     ;; If we can't find any other expansions, pass the result to SP1LOG2, which
 536     ;; tries again after expressing E as integrate(diff(e)/e).
 537     ((sp1log2 e))))
 538
 539 ;; We didn't manage to expand the expression, so make use of the fact that
 540 ;; diff(log(f(x)), x) = f'(x)/f(x) and return integrate(f'(x)/f(x), x), hoping
 541 ;; that a later stage will be able to do something useful with it.
 542 ;;
 543 ;; We have to be a little bit careful because an indefinite integral might have
 544 ;; the wrong constant term. Instead, rewrite as
 545 ;;
 546 ;;  log(f(x0+h)) = log(f(x0+h)) - log(f(x0)) + log(f(x0))
 547 ;;               = integrate(diff(log(f(x0+k)), k), k, 0, h) + log(f(x0))
 548 ;;               = integrate(diff(f(x0+k))/f(x0+k), k, 0, h) + log(f(x0))
 549 ;;
 550 ;; The "x0" about which we expand is always zero (see the code in $powerseries)
 551 (defun sp1log2 (e)
 552   (when $verbose
 553     (mtell (intl:gettext "trigreduce: failed to expand.~%~%"))
 554     (show-exp (list '(%log) e))
 555     (mtell (intl:gettext "trigreduce: try again after applying rule:~2%~M~%~%")
 556            (list '(mlabel) nil
 557                  (out-of
 558                   `((mequal)
 559                     ((%log) ,e)
 560                     ((%integrate)
 561                      ((mquotient) ((%derivative) ,e ,var 1) ,e) ,var))))))
 562   (let* ((dummy-sym ($gensym)))
 563     (m+ (list '(%log) ($limit e var 0))
 564         (list '(%integrate)
 565               (maxima-substitute dummy-sym var
 566                                  (sp1 (m// (sdiff e var) e)))
 567               dummy-sym 0 var))))
 568
 569 (defun sp1trig (e)
 570   (cond ((atom (cadr e)) (simplify e))
 571         ((eq (caaadr e) (zl-get (caar e) '$inverse)) (sp1 (cadadr e)))
 572         ((eq (caaadr e) (zl-get (zl-get (caar e) 'recip) '$inverse))
 573          (sp1 (m// (cadadr e))))
 574         ((and (null *trigred) (null *noexpand) (eq (caaadr e) 'mplus))
 575          (sp1trigex e))
 576         ( e )))
 577
 578 ;; Return the expansion of ((trigfun) ((mplus) a b)). For example sin(a+b) =
 579 ;; sin(a)cos(b) + cos(a)sin(b).
 580 (defun expand-trig-of-sum (trigfun a b)
 581   (flet ((expand-it (op f1 f2 f3 f4)
 582            (funcall op
 583                     (m* (sp1trig (list f1 a)) (sp1trig (list f2 b)))
 584                     (m* (sp1trig (list f3 a)) (sp1trig (list f4 b))))))
 585     (ecase trigfun
 586       (%sin  (expand-it #'add2* '(%sin)  '(%cos)  '(%cos)  '(%sin)))
 587       (%cos  (expand-it #'sub*  '(%cos)  '(%cos)  '(%sin)  '(%sin)))
 588       (%sinh (expand-it #'add2* '(%sinh) '(%cosh) '(%cosh) '(%sinh)))
 589       (%cosh (expand-it #'sub*  '(%cosh) '(%cosh) '(%sinh) '(%sinh))))))
 590
 591 ;; Try to expand f(a+b) where f is sin, cos, sinh or cosh.
 592 (defun sp1trigex (e)
 593   (schatchen-cond w
 594     ;; Ideally, we'd like to split the argument of the trig function into terms
 595     ;; that involve VAR and those that are free of it.
 596     ((m2 (cadr e) '((mplus) ((coeffpp) (a freevar)) ((coeffpp) (b true))))
 597      (a b)
 598
 599      ;; Make sure that if B is zero then so is A (to simplify the cond)
 600      (when (signp e b) (rotatef a b))
 601
 602      ;; Assuming we didn't just swap them, A will be free of VAR and B will
 603      ;; contain any other terms. If A is zero (because the argument of trig
 604      ;; function is a sum of terms, all of which involve VAR), then fall back on
 605      ;; a different splitting, by terms of taking the first term of B.
 606      (cond
 607        ((and (signp e a)
 608              (not (atom b))
 609              (eq (caar b) 'mplus))
 610         (expand-trig-of-sum (caar e)
 611                             (cadr b)
 612                             (if (cdddr b)
 613                                 (cons (car b) (cddr b))
 614                                 (caddr b))))
 615
 616        ;; For some weird reason, B isn't a sum. Give up.
 617        ((signp e a) e)
 618
 619        ;; Do the splitting we intended in the first place.
 620        (t
 621         (expand-trig-of-sum (caar e) a b))))
 622
 623     ;; E doesn't match f(a+b). Return it unmodified.
 624     (t nil e)))
 625
 626 (defun sp1atrig (fn exp)
 627   (cond ((atom exp)
 628          (sp1atrig2 fn exp))
 629         ((eq fn (zl-get (caar exp) '$inverse))
 630          (sp1 (cadr exp)))
 631         (t (sp1atrig2 fn exp))))
 632
 633 (defun sp1atrig2 (fn exp)
 634   (cond ((member fn '(%cot %sec %csc %coth %sech %csch) :test #'eq)
 635          (setq exp (sp1 (m// exp))
 636                fn (cdr (assoc fn '((%acot . %atan) (%asec . %acos) (%acsc . %asin)
 637                                   (%acoth . %atanh) (%asech . %acosh) (%acsch . %asinh)) :test #'eq)))))
 638   (cond ((and (null *trigred)
 639               (member fn '(%acos %acosh) :test #'eq))
 640          (m+ half%pi (list
 641                       (list (cdr (assoc fn '((%acos . %asin) (%acosh . %asinh)) :test #'eq)))
 642                       exp)))
 643         ((list (list fn) exp))))
 644
 645 (defun sin^n (%n v)
 646   (sc^n %n v (if (oddp %n) '(%sin) '(%cos)) (not (oddp %n))
 647         (m^ -1 (m+ (ash %n -1) 'k))))
 648
 649 (defun sinh^n (%n v)
 650   (if (oddp %n)
 651       (sc^n %n v '(%sinh) nil (m^ -1 'k))
 652       (if (zerop (mod %n 4))
 653           (sc^n %n v '(%cosh) t (m^ -1 'k))
 654           (m- (sc^n %n v '(%cosh) t (m- (m^ -1 'k)))))))
 655
 656 (defun cos^n (%n v)
 657   (sc^n %n v '(%cos) (not (oddp %n)) 1))
 658
 659 (defun cosh^n (%n v)
 660   (sc^n %n v '(%cosh) (not (oddp %n)) 1))
 661
 662 (defun sc^n (%n v fn fl coef)
 663   (when (minusp %n)
 664     (merror "trigreduce: internal error; %N must be nonnegative, found: ~M") %n)
 665   (m* (list '(rat) 1 (expt 2 %n))
 666       (m+ (if fl
 667               (list '(%binomial) %n (ash %n -1))
 668               0)
 669           (maxima-substitute v 'trig-var
 670                              (dosum (m+ (m* 2
 671                                             (list '(%binomial) %n 'k)
 672                                             coef
 673                                             (list fn (m* 'trig-var
 674                                                          (m+ %n (m* -2 'k))))))
 675                                     'k 0 (ash (1- %n) -1) t)))))