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