src/trigo.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 trigo)
  14
  15 (load-macsyma-macros mrgmac)
  16
  17 (def-simplifier sinh (y)
  18   (cond ((flonum-eval (mop form) y))
  19         ((big-float-eval (mop form) y))
  20         ((taylorize (mop form) (second form)))
  21         ((and $%piargs (if (zerop1 y) 0)))
  22         ((and $%iargs (multiplep y '$%i)) (mul '$%i (ftake* '%sin (coeff y '$%i 1))))
  23         ((and $triginverses (not (atom y))
  24               (let ((fcn (caar y))
  25                     (arg (cadr y)))
  26                 (cond ((eq '%asinh fcn)
  27                        arg)
  28                       ((eq '%acosh fcn)
  29                        ;; ratsimp(logarc(exponentialize(sinh(acosh(x))))),algebraic;
  30                        ;; -> sqrt(x-1)*sqrt(x+1)
  31                        (mul (power (sub arg 1) 1//2)
  32                             (power (add arg 1) 1//2)))
  33                       ((eq '%atanh fcn)
  34                        ;; radcan(logarc(exponentialize(sinh(atanh(x)))));
  35                        ;; -> x/(sqrt(1-x)*sqrt(1+x))
  36                        (div arg
  37                             (mul (power (sub 1 arg) 1//2)
  38                                  (power (add 1 arg) 1//2))))))))
  39         ((and $trigexpand (trigexpand '%sinh y)))
  40         ($exponentialize (exponentialize '%sinh y))
  41         ((and $halfangles (halfangle '%sinh y)))
  42         ((apply-reflection-simp (mop form) y $trigsign))
  43         ;((and $trigsign (mminusp* y)) (neg (ftake* '%sinh (neg y))))
  44         (t (give-up))))
  45
  46 (def-simplifier cosh (y)
  47   (cond ((flonum-eval (mop form) y))
  48         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
  49         ((taylorize (mop form) (second form)))
  50         ((and $%piargs (if (zerop1 y) 1)))
  51         ((and $%iargs (multiplep y '$%i)) (ftake* '%cos (coeff y '$%i 1)))
  52         ((and $triginverses (not (atom y))
  53               (let ((fcn (caar y))
  54                     (arg (cadr y)))
  55                 (cond ((eq '%acosh fcn)
  56                        arg)
  57                       ((eq '%asinh fcn)
  58                        ;; ex: cosh(asinh(x));
  59                        ;; ex,exponentialize,logarc;
  60                        ;; ratsimp(%),algebraic
  61                        ;; -> sqrt(x^2+1)
  62                        ;;
  63                        (sqrt1+x^2 arg))
  64                       ((eq '%atanh fcn)
  65                        ;; ex: cosh(atanh(x))
  66                        ;; radcan(logarc(exponentialize(ex)))
  67                        ;; -> 1/sqrt(1-x)/sqrt(1+x)
  68                        (div 1
  69                             (mul (power (sub 1 arg) 1//2)
  70                                  (power (add 1 arg) 1//2))))))))
  71         ((and $trigexpand (trigexpand '%cosh y)))
  72         ($exponentialize (exponentialize '%cosh y))
  73         ((and $halfangles (halfangle '%cosh y)))
  74         ((apply-reflection-simp (mop form) y $trigsign))
  75         ;((and $trigsign (mminusp* y)) (ftake* '%cosh (neg y)))
  76         (t (give-up))))
  77
  78 (def-simplifier tanh (y)
  79   (cond ((flonum-eval (mop form) y))
  80         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
  81         ((taylorize (mop form) (second form)))
  82         ((and $%piargs (if (zerop1 y) 0)))
  83         ((and $%iargs (multiplep y '$%i)) (mul '$%i (ftake* '%tan (coeff y '$%i 1))))
  84         ((and $triginverses (not (atom y))
  85               (let ((fcn (caar y))
  86                     (arg (cadr y)))
  87                 (cond ((eq '%atanh fcn)
  88                        arg)
  89                       ((eq '%asinh fcn)
  90                        ;; ratsimp(logarc(exponentialize(tanh(asinh(x))))),algebraic;
  91                        ;; --> x/sqrt(1+x^2)
  92                        (div arg (sqrt1+x^2 arg)))
  93                       ((eq '%acosh fcn)
  94                        ;; ratsimp(logarc(exponentialize(tanh(acosh(x))))),algebraic;
  95                        ;; sqrt(x-1)*sqrt(x+1)/x
  96                        (div (mul (power (sub arg 1) 1//2)
  97                                  (power (add arg 1) 1//2))
  98                             arg))))))
  99         ((and $trigexpand (trigexpand '%tanh y)))
 100         ($exponentialize (exponentialize '%tanh y))
 101         ((and $halfangles (halfangle '%tanh y)))
 102         ((apply-reflection-simp (mop form) y $trigsign))
 103         ;((and $trigsign (mminusp* y)) (neg (ftake* '%tanh (neg y))))
 104         (t (give-up))))
 105
 106 (def-simplifier coth (y)
 107   (cond ((flonum-eval (mop form) y))
 108         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 109         ((taylorize (mop form) (second form)))
 110         ((and $%piargs (if (zerop1 y) (domain-error y 'coth))))
 111         ((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (ftake* '%cot (coeff y '$%i 1))))
 112         ((and $triginverses (not (atom y)) (if (eq '%acoth (caar y)) (cadr y))))
 113         ((and $trigexpand (trigexpand '%coth y)))
 114         ($exponentialize (exponentialize '%coth y))
 115         ((and $halfangles (halfangle '%coth y)))
 116         ((apply-reflection-simp (mop form) y $trigsign))
 117         ;((and $trigsign (mminusp* y)) (neg (ftake* '%coth (neg y))))
 118         (t (give-up))))
 119
 120 (def-simplifier csch (y)
 121   (cond ((flonum-eval (mop form) y))
 122         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 123         ((taylorize (mop form) (second form)))
 124         ((and $%piargs (cond ((zerop1 y) (domain-error y 'csch)))))
 125         ((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (ftake* '%csc (coeff y '$%i 1))))
 126         ((and $triginverses (not (atom y)) (if (eq '%acsch (caar y)) (cadr y))))
 127         ((and $trigexpand (trigexpand '%csch y)))
 128         ($exponentialize (exponentialize '%csch y))
 129         ((and $halfangles (halfangle '%csch y)))
 130         ((apply-reflection-simp (mop form) y $trigsign))
 131         ;((and $trigsign (mminusp* y)) (neg (ftake* '%csch (neg y))))
 132         (t (give-up))))
 133
 134 (def-simplifier sech (y)
 135   (cond ((flonum-eval (mop form) y))
 136         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 137         ((taylorize (mop form) (second form)))
 138         ((and $%piargs (zerop1 y)) 1)
 139         ((and $%iargs (multiplep y '$%i)) (ftake* '%sec (coeff y '$%i 1)))
 140         ((and $triginverses (not (atom y)) (if (eq '%asech (caar y)) (cadr y))))
 141         ((and $trigexpand (trigexpand '%sech y)))
 142         ($exponentialize (exponentialize '%sech y))
 143         ((and $halfangles (halfangle '%sech y)))
 144         ((apply-reflection-simp (mop form) y $trigsign))
 145         ;((and $trigsign (mminusp* y)) (ftake* '%sech (neg y)))
 146         (t (give-up))))
 147
 148 (def-simplifier asin (y)
 149   (cond ((flonum-eval (mop form) y))
 150         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 151         ((taylorize (mop form) (second form)))
 152         ((and $%piargs
 153               ;; Recognize some special values
 154               (cond ((zerop1 y)
 155                      0)
 156                     ((equal 1 y)
 157                      (div '$%pi 2))
 158                     ((equal -1 y)
 159                      (div '$%pi -2))
 160                     ((alike1 y 1//2)
 161                      (div '$%pi 6))
 162                     ((alike1 y -1//2)
 163                      (div '$%pi -6))
 164                     ;; 1/sqrt(2)
 165                     ((alike1 y (power* 2 -1//2))
 166                      (div '$%pi 4))
 167                     ;; sqrt(3)/2
 168                     ((alike1 y (div (power* 3 1//2) 2))
 169                      (div '$%pi 3))
 170                     ;; -sqrt(3)/2
 171                     ((alike1 y (div (power* 3 1//2) -2))
 172                      (div '$%pi -3)))))
 173         ((and $%iargs (multiplep y '$%i)) (mul '$%i (ftake* '%asinh (coeff y '$%i 1))))
 174         ((and (not (atom y)) (member (caar y) '(%cos %sin))
 175               (if ($constantp (cadr y))
 176                   (let ((y-val (mfuncall '$mod
 177                                          (if (eq (caar y) '%sin) (cadr y) (m- %pi//2 (cadr y)))
 178                                          (m* 2 '$%pi))))
 179                     (cond ((eq (mlsp y-val %pi//2) t) y-val)
 180                           ((eq (mlsp y-val (m* 3 %pi//2)) t) (m- '$%pi y-val))
 181                           ((eq (mlsp y-val (m* 2 '$%pi)) t) (m- y-val (m* 2 '$%pi))))))))
 182         ((and (eq $triginverses t) (not (atom y)) (eq (caar y) '%sin)
 183               (if (and (member (csign (m- (cadr y) %pi//2)) '($nz $neg) :test #'eq)
 184                        (member (csign (m+ (cadr y) %pi//2)) '($pz $pos) :test #'eq))
 185                   (cadr y))))
 186         ((and (eq $triginverses '$all) (not (atom y))
 187               (if (eq '%sin (caar y)) (cadr y))))
 188         ($logarc (logarc '%asin y))
 189         ((apply-reflection-simp (mop form) y $trigsign))
 190         ;((and $trigsign (mminusp* y)) (neg (ftake* '%asin (neg y))))
 191         (t (give-up))))
 192
 193 (def-simplifier acos (y)
 194   (cond ((flonum-eval (mop form) y))
 195         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 196         ((taylorize (mop form) (second form)))
 197         ((and $%piargs
 198               ;; Recognize some special values
 199               (cond ((zerop1 y)
 200                      (div '$%pi 2))
 201                     ((equal 1 y)
 202                      0)
 203                     ((equal -1 y)
 204                      '$%pi)
 205                     ((alike1 y 1//2)
 206                      (div '$%pi 3))
 207                     ((alike1 y -1//2)
 208                      (mul '$%pi (div 2 3)))
 209                     ;; 1/sqrt(2)
 210                     ((alike1 y (power* 2 -1//2))
 211                      (div '$%pi 4))
 212                     ;; sqrt(3)/2
 213                     ((alike1 y (div (power* 3 1//2) 2))
 214                      (div '$%pi 6))
 215                     ;; -sqrt(3)/2
 216                     ((alike1 y (div (power* 3 1//2) -2))
 217                      (mul '$%pi (div 5 6))))))
 218         ((and (not (atom y)) (member (caar y) '(%cos %sin))
 219               (if ($constantp (cadr y))
 220                   (let ((y-val (mfuncall '$mod
 221                                          (if (eq (caar y) '%cos) (cadr y) (m- %pi//2 (cadr y)))
 222                                          (m* 2 '$%pi))))
 223                     (cond ((eq (mlsp y-val '$%pi) t) y-val)
 224                           ((eq (mlsp y-val (m* 2 '$%pi)) t) (m- (m* 2 '$%pi) y-val)))))))
 225         ((and (eq $triginverses '$all) (not (atom y))
 226               (if (eq '%cos (caar y)) (cadr y))))
 227         ((and (eq $triginverses t) (not (atom y)) (eq (caar y) '%cos)
 228               (if (and (member (csign (m- (cadr y) '$%pi)) '($nz $neg) :test #'eq)
 229                        (member (csign (cadr y)) '($pz $pos) :test #'eq))
 230                   (cadr y))))
 231         ($logarc (logarc '%acos y))
 232         ((apply-reflection-simp (mop form) y $trigsign))
 233         ;((and $trigsign (mminusp* y)) (sub '$%pi (ftake* '%acos (neg y))))
 234         (t (give-up))))
 235
 236 (def-simplifier acot (y)
 237   (cond ((flonum-eval (mop form) y))
 238         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 239         ((taylorize (mop form) (second form)))
 240         ((and $%piargs
 241               (cond ((zerop1 y) (div '$%pi 2))
 242                     ((equal 1 y) (div '$%pi 4))
 243                     ((equal -1 y) (div '$%pi -4))
 244                     ;; 1/sqrt(3)
 245                     ((alike1 y '((mexpt) 3 ((rat) -1 2))) (div '$%pi 3))
 246                     ;; sqrt(3)
 247                     ((alike1 y '((mexpt) 3 ((rat) 1 2))) (div '$%pi 6)))))
 248         ((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (ftake* '%acoth (coeff y '$%i 1))))
 249         ((and (not (atom y)) (member (caar y) '(%cot %tan))
 250               (if ($constantp (cadr y))
 251                   (let ((y-val (mfuncall '$mod
 252                                          (if (eq (caar y) '%cot) (cadr y) (m- %pi//2 (cadr y)))
 253                                          '$%pi)))
 254                     (cond ((eq (mlsp y-val %pi//2) t) y-val)
 255                           ((eq (mlsp y-val '$%pi) t) (m- y-val '$%pi)))))))
 256         ((and (eq $triginverses '$all) (not (atom y))
 257               (if (eq '%cot (caar y)) (cadr y))))
 258         ((and (eq $triginverses t) (not (atom y)) (eq (caar y) '%cot)
 259               (if (and (member (csign (m- (cadr y) %pi//2)) '($nz $neg) :test #'eq)
 260                        (member (csign (m+ (cadr y) %pi//2)) '($pz $pos) :test #'eq))
 261                   (cadr y))))
 262         ($logarc (logarc '%acot y))
 263         ((apply-reflection-simp (mop form) y $trigsign))
 264         ;((and $trigsign (mminusp* y)) (neg (ftake* '%acot (neg y))))
 265         (t (give-up))))
 266
 267 (def-simplifier acsc (y)
 268   (cond ((flonum-eval (mop form) y))
 269         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 270         ((taylorize (mop form) (second form)))
 271         ((and $%piargs
 272               (cond ((equal 1 y) (div '$%pi 2))
 273                     ((equal -1 y) (div '$%pi -2))
 274                     ((equal y 2) (div '$%pi 6))
 275                     ;; sqrt(2)
 276                     ((alike1 y '((mexpt) 2 ((rat) 1 2))) (div '$%pi 4))
 277                     ;; 2*sqrt(3)/3
 278                     ((alike1 y '((mtimes) 2 ((mexpt) 3 ((rat) -1 2))))
 279                      (div '$%pi 3)))))
 280         ((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (ftake* '%acsch (coeff y '$%i 1))))
 281         ((and (not (atom y)) (eq '%csc (caar y))
 282               (if ($constantp (cadr y))
 283                   (let ((y-val (mfuncall '$mod (cadr y) (m* 2 '$%pi))))
 284                     (cond ((eq (mlsp y-val %pi//2) t) y-val)
 285                           ((eq (mlsp y-val (m* 3 %pi//2)) t) (m- '$%pi y-val))
 286                           ((eq (mlsp y-val (m* 2 '$%pi)) t) (m- y-val (m* 2 '$%pi))))))))
 287         ((and (eq $triginverses '$all) (not (atom y))
 288               (if (eq '%csc (caar y)) (cadr y))))
 289         ($logarc (logarc '%acsc y))
 290         ((apply-reflection-simp (mop form) y $trigsign))
 291         ;((and $trigsign (mminusp* y)) (neg (ftake* '%acsc (neg y))))
 292         (t (give-up))))
 293
 294 (def-simplifier asec (y)
 295   (cond ((flonum-eval (mop form) y))
 296         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 297         ((taylorize (mop form) (second form)))
 298         ((and $%piargs
 299               (cond ((equal 1 y) 0)
 300                     ((equal -1 y) '$%pi)
 301                     ((equal 2 y) (div '$%pi 3))
 302                     ;; sqrt(2)
 303                     ((alike1 y '((mexpt) 2 ((rat) 1 2))) (div '$%pi 4))
 304                     ;; 2/sqrt(3)
 305                     ((alike1 y '((mtimes) 2 ((mexpt) 3 ((rat) -1 2))))
 306                      (div '$%pi 6)))))
 307         ((and (not (atom y)) (eq '%sec (caar y))
 308               (if ($constantp (cadr y))
 309                   (let ((y-val (mfuncall '$mod (cadr y) (m* 2 '$%pi))))
 310                     (cond ((eq (mlsp y-val '$%pi) t) y-val)
 311                           ((eq (mlsp y-val (m* 2 '$%pi)) t) (m- (m* 2 '$%pi) y-val)))))))
 312         ((and (eq $triginverses '$all) (not (atom y))
 313               (if (eq '%sec (caar y)) (cadr y))))
 314         ($logarc (logarc '%asec y))
 315         ((apply-reflection-simp (mop form) y $trigsign))
 316         ;;((and $trigsign (mminusp* y)) (sub '$%pi (ftake* '%asec (neg y))))
 317         (t (give-up))))
 318
 319 (def-simplifier asinh (y)
 320   (cond ((flonum-eval (mop form) y))
 321         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 322         ((taylorize (mop form) (second form)))
 323         ((and $%piargs (if (zerop1 y) y)))
 324         ((and $%iargs (multiplep y '$%i)) (mul '$%i (ftake* '%asin (coeff y '$%i 1))))
 325         ((and (eq $triginverses '$all) (not (atom y))
 326               (if (eq '%sinh (caar y)) (cadr y))))
 327         ($logarc (logarc '%asinh y))
 328         ((apply-reflection-simp (mop form) y $trigsign))
 329         ;((and $trigsign (mminusp* y)) (neg (ftake* '%asinh (neg y))))
 330         (t (give-up))))
 331
 332 (def-simplifier acosh (y)
 333   (cond ((flonum-eval (mop form) y))
 334         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 335         ((taylorize (mop form) (second form)))
 336         ((and $%piargs (if (equal y 1) 0)))
 337         ((and (eq $triginverses '$all) (not (atom y))
 338               (if (eq '%cosh (caar y)) (cadr y))))
 339         ($logarc (logarc '%acosh y))
 340         (t (give-up))))
 341
 342 (def-simplifier atanh (y)
 343   (cond ((flonum-eval (mop form) y))
 344         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 345         ((taylorize (mop form) (second form)))
 346         ((and $%piargs (cond ((zerop1 y) 0)
 347                              ((or (equal y 1) (equal y -1)) (domain-error y 'atanh)))))
 348         ((and $%iargs (multiplep y '$%i)) (mul '$%i (ftake* '%atan (coeff y '$%i 1))))
 349         ((and (eq $triginverses '$all) (not (atom y))
 350               (if (eq '%tanh (caar y)) (cadr y))))
 351         ($logarc (logarc '%atanh y))
 352         ((apply-reflection-simp (mop form) y $trigsign))
 353         ;((and $trigsign (mminusp* y)) (neg (ftake* '%atanh (neg y))))
 354         (t (give-up))))
 355
 356 (def-simplifier acoth (y)
 357   (cond ((flonum-eval (mop form) y))
 358         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 359         ((taylorize (mop form) (second form)))
 360         ((and $%piargs (if (or (equal y 1) (equal y -1)) (domain-error y 'acoth))))
 361         ((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (ftake* '%acot (coeff y '$%i 1))))
 362         ((and (eq $triginverses '$all) (not (atom y))
 363               (if (eq '%coth (caar y)) (cadr y))))
 364         ($logarc (logarc '%acoth y))
 365         ((apply-reflection-simp (mop form) y $trigsign))
 366         ;((and $trigsign (mminusp* y)) (neg (ftake* '%acoth (neg y))))
 367         (t (give-up))))
 368
 369 (def-simplifier acsch (y)
 370   (cond ((flonum-eval (mop form) y))
 371         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 372         ((taylorize (mop form) (second form)))
 373         ((and $%piargs (if (zerop1 y) (domain-error y 'acsch))))
 374         ((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (ftake* '%acsc (coeff y '$%i 1))))
 375         ((and (eq $triginverses '$all) (not (atom y))
 376               (if (eq '%csch (caar y)) (cadr y))))
 377         ($logarc (logarc '%acsch y))
 378         ((apply-reflection-simp (mop form) y $trigsign))
 379         ;((and $trigsign (mminusp* y)) (neg (ftake* '%acsch (neg y))))
 380         (t (give-up))))
 381
 382 (def-simplifier asech (y)
 383   (cond ((flonum-eval (mop form) y))
 384         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 385         ((taylorize (mop form) (second form)))
 386         ((and $%piargs (cond ((equal y 1) 0)
 387                              ((zerop1 y) (domain-error y 'asech)))))
 388         ((and (eq $triginverses '$all) (not (atom y))
 389               (if (eq '%sech (caar y)) (cadr y))))
 390         ($logarc (logarc '%asech y))
 391         ((apply-reflection-simp (mop form) y $trigsign))
 392         ;((and $trigsign (mminusp* y)) (ftake* '%asech (neg y)))
 393         (t (give-up))))
 394
 395 (declare-top (special $trigexpandplus $trigexpandtimes))
 396
 397 (defmfun ($trigexpand :properties ((evfun t))) (e)
 398   (cond ((atom e) e)
 399         ((specrepp e) ($trigexpand (specdisrep e)))
 400         ((trigexpand (caar e) (cadr e)))
 401         (t (recur-apply #'$trigexpand e))))
 402
 403 (defun trigexpand (op arg)
 404   (cond ((atom arg) nil)
 405         ((and $trigexpandplus (eq 'mplus (caar arg)))
 406          (cond ((eq '%sin op) (sin/cos-plus (cdr arg) 1 '%sin '%cos -1))
 407                ((eq '%cos op) (sin/cos-plus (cdr arg) 0 '%sin '%cos -1))
 408                ((eq '%tan op) (tan-plus (cdr arg) '%tan -1))
 409                ((eq '%cot op) (cot-plus (cdr arg) '%cot -1))
 410                ((eq '%csc op) (csc/sec-plus (cdr arg) 1 '%csc '%sec -1))
 411                ((eq '%sec op) (csc/sec-plus (cdr arg) 0 '%csc '%sec -1))
 412                ((eq '%sinh op) (sin/cos-plus (cdr arg) 1 '%sinh '%cosh 1))
 413                ((eq '%cosh op) (sin/cos-plus (cdr arg) 0 '%sinh '%cosh 1))
 414                ((eq '%tanh op) (tan-plus (cdr arg) '%tanh 1))
 415                ((eq '%coth op) (cot-plus (cdr arg) '%coth 1))
 416                ((eq '%csch op) (csc/sec-plus (cdr arg) 1 '%csch '%sech 1))
 417                ((eq '%sech op) (csc/sec-plus (cdr arg) 0 '%csch '%sech 1))))
 418         ((and $trigexpandtimes (eq 'mtimes (caar arg)) (fixnump (cadr arg)))
 419          (cond ((eq '%sin op) (sin/cos-times (cddr arg) 1 (cadr arg) '%sin '%cos -1))
 420                ((eq '%cos op) (sin/cos-times (cddr arg) 0 (cadr arg) '%sin '%cos -1))
 421                ((eq '%tan op) (tan-times (cddr arg) (cadr arg) '%tan -1))
 422                ((eq '%cot op) (cot-times (cddr arg) (cadr arg) '%cot -1))
 423                ((eq '%csc op) (csc/sec-times (cddr arg) 1 (cadr arg) '%csc '%sec -1))
 424                ((eq '%sec op) (csc/sec-times (cddr arg) 0 (cadr arg) '%csc '%sec -1))
 425                ((eq '%sinh op) (sin/cos-times (cddr arg) 1 (cadr arg) '%sinh '%cosh 1))
 426                ((eq '%cosh op) (sin/cos-times (cddr arg) 0 (cadr arg) '%sinh '%cosh 1))
 427                ((eq '%tanh op) (tan-times (cddr arg) (cadr arg) '%tanh 1))
 428                ((eq '%coth op) (cot-times (cddr arg) (cadr arg) '%coth 1))
 429                ((eq '%csch op) (csc/sec-times (cddr arg) 1 (cadr arg) '%csch '%sech 1))
 430                ((eq '%sech op) (csc/sec-times (cddr arg) 0 (cadr arg) '%csch '%sech 1))))))
 431
 432 (defun sin/cos-plus (l n f1 f2 flag)
 433   (do ((i n (+ 2 i))
 434        (len (length l))
 435        (sign 1 (* flag sign))
 436        (result))
 437       ((> i len) (simplify (cons '(mplus) result)))
 438     (setq result (mpc (cond ((minusp sign) '(-1 (mtimes)))
 439                             (t '((mtimes)))) l result f1 f2 len i))))
 440
 441 (defun tan-plus (l f flag)
 442   (do ((i 1 (+ 2 i))
 443        (sign 1 (* flag sign))
 444        (len (length l))
 445        (num)
 446        (den (list 1)))
 447       ((> i len) (div* (cons '(mplus) num) (cons '(mplus) den)))
 448     (setq num (mpc1 (list sign '(mtimes)) l num f len i)
 449           den (cond ((= len i) den)
 450                     (t (mpc1 (list (* flag sign) '(mtimes)) l den f len (1+ i)))))))
 451
 452 (defun cot-plus (l f flag)
 453   (do ((i (length l) (- i 2)) (len (length l)) (sign 1 (* flag sign)) (num) (den))
 454       ((< i 0) (div* (cons '(mplus) num) (cons '(mplus) den)))
 455     (setq num (mpc1 (list sign '(mtimes)) l num f len i)
 456           den (cond ((= 0 i) den)
 457                     (t (mpc1 (list sign '(mtimes)) l den f len (1- i)))))))
 458
 459 (defun csc/sec-plus (l n f1 f2 flag)
 460   (div* (do ((l l (cdr l))
 461              (result))
 462             ((null l) (cons '(mtimes) result))
 463           (setq result (cons (ftake* f1 (car l)) (cons (ftake* f2 (car l)) result))))
 464         (sin/cos-plus l n f2 f1 flag)))
 465
 466 (defun sin/cos-times (l m n f1 f2 flag)
 467   ;; Assume m,n < 2^17, but Binom may become big
 468   ;; Flag is 1 or -1
 469   (setq f1 (ftake* f1 (cons '(mtimes) l)) f2 (ftake* f2 (cons '(mtimes) l)))
 470   (do ((i m (+ 2 i))
 471        (end (abs n))
 472        (result)
 473        (binom (cond ((= 0 m) 1)
 474                     (t (abs n)))
 475               (quotient (* flag (- end i 1) (- end i) binom) (* (+ 2 i) (1+ i)))))
 476       ((> i end) (setq result (simplify (cons '(mplus) result)))
 477        (cond ((and (= 1 m) (minusp n)) (neg result)) (t result)))
 478     (setq result (cons (mul binom (power f1 i) (power f2 (- end i))) result))))
 479
 480 (defun tan-times (l n f flag)
 481   (setq f (ftake* f (cons '(mtimes) l)))
 482   (do ((i 1 (+ 2 i))
 483        (end (abs n))
 484        (num)
 485        (den (list 1))
 486        (binom (abs n) (quotient (* (- end i 1) binom) (+ 2 i))))
 487       ((> i end) (setq num (div* (cons '(mplus) num) (cons '(mplus) den)))
 488        (cond ((minusp n) (neg num))
 489              (t num)))
 490     (setq num (cons (mul binom (power f i)) num)
 491           den (cond ((= end i) den)
 492                     (t (cons (mul (setq binom (truncate (* flag (- end i) binom) (1+ i)))
 493                                   (power f (1+ i)))
 494                              den))))))
 495
 496 (defun cot-times (l n f flag)
 497   (setq f (ftake* f (cons '(mtimes) l)))
 498   (do ((i (abs n) (- i 2))
 499        (end (abs n))
 500        (num)
 501        (den)
 502        (binom 1 (truncate (* flag (1- i) binom) (- end i -2))))
 503       ((< i 0) (setq num (div* (cons '(mplus) num) (cons '(mplus) den)))
 504        (if (minusp n) (neg num) num))
 505     (setq num (cons (mul binom (power f i)) num)
 506           den (if (= 0 i)
 507                   den
 508                   (cons (mul (setq binom (truncate (* i binom) (- end i -1))) (power f (1- i))) den)))))
 509
 510 (defun csc/sec-times (l m n f1 f2 flag)
 511   (div* (mul (power (ftake* f1 (cons '(mtimes) l)) (abs n))
 512              (power (ftake* f2 (cons '(mtimes) l)) (abs n)))
 513         (sin/cos-times l m n f2 f1 flag)))
 514
 515 (defun mpc (dl ul result f1 f2 di ui)
 516   (cond ((= 0 ui)
 517          (cons (revappend dl (mapcar #'(lambda (l) (ftake* f2 l)) ul)) result))
 518         ((= di ui)
 519          (cons (revappend dl (mapcar #'(lambda (l) (ftake* f1 l)) ul)) result))
 520         (t (mpc (cons (ftake* f1 (car ul)) dl) (cdr ul)
 521                 (mpc (cons (ftake* f2 (car ul)) dl)
 522                      (cdr ul) result f1 f2 (1- di) ui) f1 f2
 523                 (1- di) (1- ui)))))
 524
 525 (defun mpc1 (dl ul result f di ui)
 526   (cond ((= 0 ui) (cons (reverse dl) result))
 527         ((= di ui)
 528          (cons (revappend dl (mapcar #'(lambda (l) (ftake* f l)) ul)) result))
 529         (t (mpc1 (cons (ftake* f (car ul)) dl) (cdr ul)
 530                  (mpc1 dl (cdr ul) result f (1- di) ui) f
 531                  (1- di) (1- ui)))))
 532
 533 ;; Local Modes:
 534 ;; Mode: LISP
 535 ;; Comment Col: 40
 536 ;; End: