src/trgred.lisp

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