;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;     The data in this file contains enhancments.                    ;;;;;
;;;                                                                    ;;;;;
;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
;;;     All rights reserved                                            ;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;     (c) Copyright 1982 Massachusetts Institute of Technology         ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(in-package "MAXIMA")
(macsyma-module trigo)

(load-macsyma-macros mrgmac)

(declare-top (genprefix tri)
	     (special varlist errorsw)
	     (flonum (tan) (cot) (sec) (csc)
		     (atan2) (atan1) (acot)
		     (sinh) (cosh) (tanh) (coth) (csch) (sech)
		     (asinh) (acsch)
		     (t//$ flonum flonum notype))
	     (*expr $bfloat teval signum1 zerop1 islinear
		    timesk addk maxima-integerp evod logarc mevenp halfangle coeff))

(declare-top (splitfile hyper))

(defmfun simp-%sinh (form y z) 
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (sinh y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs (if (zerop1 y) 0)))
	((and $%iargs (multiplep y '$%i)) (mul '$%i (cons-exp '%sin (coeff y '$%i 1))))
	((and $triginverses (not (atom y)) (if (eq '%asinh (caar y)) (cadr y))))
	((and $trigexpand (trigexpand '%sinh y)))
	($exponentialize (exponentialize '%sinh y))
	((and $halfangles (halfangle '%sinh y)))
	((and $trigsign (mminusp* y)) (neg (cons-exp '%sinh (neg y))))
	(t (eqtest (list '(%sinh) y) form))))

(defmfun simp-%cosh (form y z) 
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (cosh y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs (if (zerop1 y) 1)))
	((and $%iargs (multiplep y '$%i)) (cons-exp '%cos (coeff y '$%i 1)))
	((and $triginverses (not (atom y)) (if (eq '%acosh (caar y)) (cadr y))))
	((and $trigexpand (trigexpand '%cosh y)))
	($exponentialize (exponentialize '%cosh y))
	((and $halfangles (halfangle '%cosh y)))
	((and $trigsign (mminusp* y)) (cons-exp '%cosh (neg y)))
	(t (eqtest (list '(%cosh) y) form))))

(defmfun simp-%tanh (form y z)
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (tanh y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs (if (zerop1 y) 0)))
	((and $%iargs (multiplep y '$%i)) (mul '$%i (cons-exp '%tan (coeff y '$%i 1))))
	((and $triginverses (not (atom y)) (if (eq '%atanh (setq z (caar y))) (cadr y))))
	((and $trigexpand (trigexpand '%tanh y)))
	($exponentialize (exponentialize '%tanh y))
	((and $halfangles (halfangle '%tanh y)))
	((and $trigsign (mminusp* y)) (neg (cons-exp '%tanh (neg y))))
	(t (eqtest (list '(%tanh) y) form))))

(defmfun simp-%coth (form y z)
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (coth y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs (if (zerop1 y) (dbz-err1 'coth))))
	((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (cons-exp '%cot (coeff y '$%i 1))))
	((and $triginverses (not (atom y)) (if (eq '%acoth (caar y)) (cadr y))))
	((and $trigexpand (trigexpand '%coth y)))
	($exponentialize (exponentialize '%coth y))
	((and $halfangles (halfangle '%coth y)))
	((and $trigsign (mminusp* y)) (neg (cons-exp '%coth (neg y))))
	(t (eqtest (list '(%coth) y) form))))

(defmfun simp-%csch (form y z)
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (csch y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs (cond ((zerop1 y) (dbz-err1 'csch)))))
	((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (cons-exp '%csc (coeff y '$%i 1))))
	((and $triginverses (not (atom y)) (if (eq '%acsch (caar y)) (cadr y))))
	((and $trigexpand (trigexpand '%csch y)))
	($exponentialize (exponentialize '%csch y))
	((and $halfangles (halfangle '%csch y)))
	((and $trigsign (mminusp* y)) (neg (cons-exp '%csch (neg y))))
	(t (eqtest (list '(%csch) y) form))))

(defmfun simp-%sech (form y z)
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (sech y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs (zerop1 y)) 1)
	((and $%iargs (multiplep y '$%i)) (cons-exp '%sec (coeff y '$%i 1)))
	((and $triginverses (not (atom y)) (if (eq '%asech (caar y)) (cadr y))))
	((and $trigexpand (trigexpand '%sech y)))
	($exponentialize (exponentialize '%sech y))
	((and $halfangles (halfangle '%sech y)))
	((and $trigsign (mminusp* y)) (cons-exp '%sech (neg y)))
	(t (eqtest (list '(%sech) y) form))))

(declare-top (splitfile atrig))

(defmfun simp-%asin (form y z) 
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (asin y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs 
	      (cond ((zerop1 y) 0) ((equal 1 y) %pi//2) ((equal -1 y) (neg %pi//2))
		    ((alike1 y 1//2) (mul '((rat simp) 1 6) '$%pi)))))
	((and $%iargs (multiplep y '$%i)) (mul '$%i (cons-exp '%asinh (coeff y '$%i 1))))
	((and (eq $triginverses '$all) (not (atom y)) (if (eq '%sin (caar y)) (cadr y))))
	($logarc (logarc '%asin y))
	((and $trigsign (mminusp* y)) (neg (cons-exp '%asin (neg y))))
	(t (eqtest (list '(%asin) y) form))))

(defmfun simp-%acos (form y z)
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (acos y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs 
	      (cond ((zerop1 y) %pi//2) ((equal 1 y) 0) ((equal -1 y) '$%pi)
		    ((alike1 y 1//2) (mul '((rat simp) 1 3) '$%pi)))))
	((and (eq $triginverses '$all) (not (atom y))
	      (if (eq '%cos (caar y)) (cadr y))))
	($logarc (logarc '%acos y))
	((and $trigsign (mminusp* y)) (sub '$%pi (cons-exp '%acos (neg y))))
	(t (eqtest (list '(%acos) y) form))))

(defmfun simp-%acot (form y z)
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (acot y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs
	      (cond ((zerop1 y) %pi//2) ((equal 1 y) %pi//4) ((equal -1 y) (neg %pi//4)))))
	((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (cons-exp '%acoth (coeff y '$%i 1))))
	((and (eq $triginverses '$all) (not (atom y))
	      (if (eq '%cot (caar y)) (cadr y))))
	($logarc (logarc '%acot y))
	((and $trigsign (mminusp* y)) (neg (cons-exp '%acot (neg y))))
	(t (eqtest (list '(%acot) y) form))))

(defmfun simp-%acsc (form y z)
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (acsc y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs
	      (cond ((equal 1 y) %pi//2) ((equal -1 y) (neg %pi//2)))))
	((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (cons-exp '%acsch (coeff y '$%i 1))))
	((and (eq $triginverses '$all) (not (atom y))
	      (if (eq '%csc (caar y)) (cadr y))))
	($logarc (logarc '%acsc y))
	((and $trigsign (mminusp* y)) (neg (cons-exp '%acsc (neg y))))
	(t (eqtest (list '(%acsc) y) form))))

(defmfun simp-%asec (form y z)
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (asec y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs 
	      (cond ((equal 1 y) 0) ((equal -1 y) '$%pi))))
	((and (eq $triginverses '$all) (not (atom y))
	      (if (eq '%sec (caar y)) (cadr y))))
	($logarc (logarc '%asec y))
	((and $trigsign (mminusp* y)) (sub '$%pi (cons-exp '%asec (neg y))))
	(t (eqtest (list '(%asec) y) form))))

(declare-top (splitfile ahyper))

(defmfun simp-%asinh (form y z)
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (asinh y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs (if (zerop1 y) y)))
	((and $%iargs (multiplep y '$%i)) (mul '$%i (cons-exp '%asin (coeff y '$%i 1))))
	((and (eq $triginverses '$all) (not (atom y))
	      (if (eq '%sinh (caar y)) (cadr y))))
	($logarc (logarc '%asinh y))
	((and $trigsign (mminusp* y)) (neg (cons-exp '%asinh (neg y))))
	(t (eqtest (list '(%asinh) y) form))))

(defmfun simp-%acosh (form y z)
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (acosh y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs (if (equal y 1) 0)))
	((and (eq $triginverses '$all) (not (atom y))
	      (if (eq '%cosh (caar y)) (cadr y))))
	($logarc (logarc '%acosh y))
	(t (eqtest (list '(%acosh) y) form))))

(defmfun simp-%atanh (form y z)
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (atanh y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs (cond ((zerop1 y) 0)
			     ((or (equal y 1) (equal y -1)) (dbz-err1 'atanh)))))
	((and $%iargs (multiplep y '$%i)) (mul '$%i (cons-exp '%atan (coeff y '$%i 1))))
	((and (eq $triginverses '$all) (not (atom y))
	      (if (eq '%tanh (caar y)) (cadr y))))
	($logarc (logarc '%atanh y))
	((and $trigsign (mminusp* y)) (neg (cons-exp '%atanh (neg y))))
	(t (eqtest (list '(%atanh) y) form))))

(defmfun simp-%acoth (form y z)
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (acoth y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs (if (or (equal y 1) (equal y -1)) (dbz-err1 'acoth))))
	((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (cons-exp '%acot (coeff y '$%i 1))))
	((and (eq $triginverses '$all) (not (atom y))
	      (if (eq '%coth (caar y)) (cadr y))))
	($logarc (logarc '%acoth y))
	((and $trigsign (mminusp* y)) (neg (cons-exp '%acoth (neg y))))
	(t (eqtest (list '(%acoth) y) form))))

(defmfun simp-%acsch (form y z)
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (acsch y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs (if (zerop1 y) (dbz-err1 'acsch))))
	((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (cons-exp '%acsc (coeff y '$%i 1))))
	((and (eq $triginverses '$all) (not (atom y))
	      (if (eq '%csch (caar y)) (cadr y))))
	($logarc (logarc '%acsch y))
	((and $trigsign (mminusp* y)) (neg (cons-exp '%acsch (neg y))))
	(t (eqtest (list '(%acsch) y) form))))

(defmfun simp-%asech (form y z)
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((or (floatp y) (and $numer (integerp y))) (asech y))
	(($bfloatp y) ($bfloat form))
	((and $%piargs (cond ((equal y 1) 0)
			     ((zerop1 y) (dbz-err1 'asech)))))
	((and (eq $triginverses '$all) (not (atom y))
	      (if (eq '%sech (caar y)) (cadr y))))
	($logarc (logarc '%asech y))
	((and $trigsign (mminusp* y)) (cons-exp '%asech (neg y)))
	(t (eqtest (list '(%asech) y) form))))

(declare-top (splitfile trigex) (special $trigexpandplus $trigexpandtimes))

(defmfun $trigexpand (e)
  (cond ((atom e) e)
	((specrepp e) ($trigexpand (specdisrep e)))
	((trigexpand (caar e) (cadr e)))
	(t (recur-apply #'$trigexpand e))))

(defmfun trigexpand (op arg)
  (cond ((atom arg) nil)
	((and $trigexpandplus (eq 'mplus (caar arg)))
	 (cond ((eq '%sin op) (sin\Cos-plus (cdr arg) 1 '%sin '%cos -1))
	       ((eq '%cos op) (sin\Cos-plus (cdr arg) 0 '%sin '%cos -1))
	       ((eq '%tan op) (tan-plus (cdr arg) '%tan -1))
	       ((eq '%cot op) (cot-plus (cdr arg) '%cot -1))
	       ((eq '%csc op) (csc\Sec-plus (cdr arg) 1 '%csc '%sec -1))
	       ((eq '%sec op) (csc\Sec-plus (cdr arg) 0 '%csc '%sec -1))
	       ((eq '%sinh op) (sin\Cos-plus (cdr arg) 1 '%sinh '%cosh 1))
	       ((eq '%cosh op) (sin\Cos-plus (cdr arg) 0 '%sinh '%cosh 1))
	       ((eq '%tanh op) (tan-plus (cdr arg) '%tanh 1))
	       ((eq '%coth op) (cot-plus (cdr arg) '%coth 1))
	       ((eq '%csch op) (csc\Sec-plus (cdr arg) 1 '%csch '%sech 1))
	       ((eq '%sech op) (csc\Sec-plus (cdr arg) 0 '%csch '%sech 1))))
	((and $trigexpandtimes (eq 'mtimes (caar arg)) (eq (ml-typep (cadr arg)) 'fixnum))
	 (cond ((eq '%sin op) (sin\Cos-times (cddr arg) 1 (cadr arg) '%sin '%cos -1))
	       ((eq '%cos op) (sin\Cos-times (cddr arg) 0 (cadr arg) '%sin '%cos -1))
	       ((eq '%tan op) (tan-times (cddr arg) (cadr arg) '%tan -1))
	       ((eq '%cot op) (cot-times (cddr arg) (cadr arg) '%cot -1))
	       ((eq '%csc op) (csc\Sec-times (cddr arg) 1 (cadr arg) '%csc '%sec -1))
	       ((eq '%sec op) (csc\Sec-times (cddr arg) 0 (cadr arg) '%csc '%sec -1))
	       ((eq '%sinh op) (sin\Cos-times (cddr arg) 1 (cadr arg) '%sinh '%cosh 1))
	       ((eq '%cosh op) (sin\Cos-times (cddr arg) 0 (cadr arg) '%sinh '%cosh 1))
	       ((eq '%tanh op) (tan-times (cddr arg) (cadr arg) '%tanh 1))
	       ((eq '%coth op) (cot-times (cddr arg) (cadr arg) '%coth 1))
	       ((eq '%csch op) (csc\Sec-times (cddr arg) 1 (cadr arg) '%csch '%sech 1))
	       ((eq '%sech op) (csc\Sec-times (cddr arg) 0 (cadr arg) '%csch '%sech 1))))))


(defun sin\Cos-plus (l n f1 f2 flag)
  (do ((i n (f+ 2 i)) (len (length l)) (sign 1 (f* flag sign)) (result))
      ((> i len) (simplify (cons '(mplus) result)))
    (setq result (mpc (cond ((minusp sign) '(-1 (mtimes))) (t '((mtimes)))) l result f1 f2 len i))))

(defun tan-plus (l f flag) 
  (do ((i 1 (f+ 2 i)) (sign 1 (f* flag sign)) (len (length l)) (num) (den (list 1)))
      ((> i len) (div* (cons '(mplus) num) (cons '(mplus) den)))
    (setq num (mpc1 (list sign '(mtimes)) l num f len i)
	  den (cond ((= len i) den)
		    (t (mpc1 (list (f* flag sign) '(mtimes)) l den f len (f1+ i)))))))

(defun cot-plus (l f flag)
  (do ((i (length l) (f- i 2)) (len (length l)) (sign 1 (f* flag sign)) (num) (den))
      ((< i 0) (div* (cons '(mplus) num) (cons '(mplus) den)))
    (setq num (mpc1 (list sign '(mtimes)) l num f len i)
	  den (cond ((= 0 i) den)
		    (t (mpc1 (list sign '(mtimes)) l den f len (f1- i)))))))

(defun csc\Sec-plus (l n f1 f2 flag)
  (div* (do ((l l (cdr l)) (result)) ((null l) (cons '(mtimes) result))
	  (setq result (cons (cons-exp f1 (car l)) (cons (cons-exp f2 (car l)) result))))
	(sin\Cos-plus l n f1 f2 flag)))

(defun sin\Cos-times (l m n f1 f2 flag)
  ;; Assume m,n < 2^17, but Binom may become big
  ;; Flag is 1 or -1
  (setq f1 (cons-exp f1 (cons '(mtimes) l)) f2 (cons-exp f2 (cons '(mtimes) l)))
  (do ((i m (f+ 2 i)) (end (abs n)) (result)
       (binom (cond ((= 0 m) 1) (t (abs n))) (quotient (times (f* flag (f- end i 1) (f- end i)) binom) (f* (f+ 2 i) (f1+ i)))))
      ((> i end) (setq result (simplify (cons '(mplus) result)))
       (cond ((and (= 1 m) (minusp n)) (neg result)) (t result)))
    (setq result (cons (mul binom (power f1 i) (power f2 (f- end i))) result))))

(defun tan-times (l n f flag)
  (setq f (cons-exp f (cons '(mtimes) l)))
  (do ((i 1 (f+ 2 i)) (end (abs n)) (num) (den (list 1))
       (binom (abs n) (quotient (times (f- end i 1) binom) (f+ 2 i))))
      ((> i end) (setq num (div* (cons '(mplus) num) (cons '(mplus) den)))
       (cond ((minusp n) (neg num)) (t num)))
    (setq num (cons (mul binom (power f i)) num) 
	  den (cond ((= end i) den)
		    (t (cons (mul (setq binom (// (f* flag (f- end i) binom) (f1+ i)))
				  (power f (f1+ i)))
			     den))))))

(defun cot-times (l n f flag)
  (setq f (cons-exp f (cons '(mtimes) l)))
  (do ((i (abs n) (f- i 2)) (end (abs n)) (num) (den)
       (binom 1 (// (f* flag (f1- i) binom) (f- end i -2))))
      ((< i 0) (setq num (div* (cons '(mplus) num) (cons '(mplus) den)))
       (if (minusp n) (neg num) num))
    (setq num (cons (mul binom (power f i)) num)
	  den (if (= 0 i) den
		  (cons (mul (setq binom (// (f* i binom) (f- end i -1))) (power f (f1- i))) den)))))

(defun csc\Sec-times (l m n f1 f2 flag)
  (div* (mul (power (cons-exp f1 (cons '(mtimes) l)) (abs n))
	     (power (cons-exp f2 (cons '(mtimes) l)) (abs n)))
	(sin\Cos-times l m n f1 f2 flag)))

(defun mpc (dl ul result f1 f2 di ui) 
  (cond ((= 0 ui)
	 (cons (reconc dl (mapcar #'(lambda (l) (cons-exp f2 l)) ul))
	       result))
	((= di ui)
	 (cons (reconc dl (mapcar #'(lambda (l) (cons-exp f1 l)) ul))
	       result))
	(t (mpc (cons (cons-exp f1 (car ul)) dl) (cdr ul)
		(mpc (cons (cons-exp f2 (car ul)) dl)
		     (cdr ul) result f1 f2 (f1- di) ui) f1 f2
		(f1- di) (f1- ui)))))

(defun mpc1 (dl ul result f di ui) 
  (cond ((= 0 ui) (cons (reverse dl) result))
	((= di ui)
	 (cons (reconc dl (mapcar #'(lambda (l) (cons-exp f l)) ul)) result))
	(t (mpc1 (cons (cons-exp f (car ul)) dl) (cdr ul)
		 (mpc1 dl (cdr ul) result f (f1- di) ui) f
		 (f1- di) (f1- ui)))))

;; Local Modes:
;; Mode: LISP
;; Comment Col: 40
;; End:
