1;;  Copyright 2005, 2006 by Barton Willis
2
3;;  This is free software; you can redistribute it and/or
4;;  modify it under the terms of the GNU General Public License,
5;;  http://www.gnu.org/copyleft/gpl.html.
6
7;; This software has NO WARRANTY, not even the implied warranty of
8;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
9
10(in-package :maxima)
11
12(macsyma-module conjugate)
13
14($put '$conjugate 1 '$version)
15;; Let's remove built-in symbols from list for user-defined properties.
16(setq $props (remove '$conjugate $props))
17
18(defprop $conjugate tex-postfix tex)
19(defprop $conjugate ("^\\star") texsym)
20(defprop $conjugate 160. tex-lbp)
21(defprop $conjugate simp-conjugate operators)
22
23(eval-when
24    #+gcl (load eval)
25    #-gcl (:load-toplevel :execute)
26    (let (($context '$global) (context '$global))
27      (meval '(($declare) $conjugate $complex))
28      ;; Let's remove built-in symbols from list for user-defined properties.
29      (setq $props (remove '$conjugate $props))))
30
31;; When a function commutes with the conjugate, give the function the
32;; commutes-with-conjugate property. The log function commutes with
33;; the conjugate on all of C except on the negative real axis. Thus
34;; log does not get the commutes-with-conjugate property.  Instead,
35;; log gets the conjugate-function property.
36
37;; What important functions have I missed?
38
39;; (1) Arithmetic operators
40
41(setf (get 'mplus 'commutes-with-conjugate) t)
42(setf (get 'mtimes 'commutes-with-conjugate) t)
43;(setf (get 'mnctimes 'commutes-with-conjugate) t) ;; generally I think users will want this
44(setf (get '%signum 'commutes-with-conjugate) t) ;; x=/=0, conjugate(signum(x)) = conjugate(x/abs(x)) = signum(conjugate(x))
45;; Trig-like functions and other such functions
46
47(setf (get '%cosh 'commutes-with-conjugate) t)
48(setf (get '%sinh 'commutes-with-conjugate) t)
49(setf (get '%tanh 'commutes-with-conjugate) t)
50(setf (get '%sech 'commutes-with-conjugate) t)
51(setf (get '%csch 'commutes-with-conjugate) t)
52(setf (get '%coth 'commutes-with-conjugate) t)
53(setf (get '%cos 'commutes-with-conjugate) t)
54(setf (get '%sin 'commutes-with-conjugate) t)
55(setf (get '%tan 'commutes-with-conjugate) t)
56(setf (get '%sec 'commutes-with-conjugate) t)
57(setf (get '%csc 'commutes-with-conjugate) t)
58(setf (get '%cot 'commutes-with-conjugate) t)
59(setf (get '$atan2 'commutes-with-conjugate) t)
60
61(setf (get '%jacobi_cn 'commutes-with-conjugate) t)
62(setf (get '%jacobi_sn 'commutes-with-conjugate) t)
63(setf (get '%jacobi_dn 'commutes-with-conjugate) t)
64
65(setf (get '%gamma 'commutes-with-conjugate) t)
66(setf (get '$pochhammer 'commutes-with-conjugate) t)
67
68;; Collections
69
70(setf (get '$matrix 'commutes-with-conjugate) t)
71(setf (get 'mlist 'commutes-with-conjugate) t)
72(setf (get '$set 'commutes-with-conjugate) t)
73
74;; Relations
75
76(setf (get 'mequal 'commutes-with-conjugate) t)
77(setf (get 'mnotequal 'commutes-with-conjugate) t)
78(setf (get '%transpose 'commutes-with-conjugate) t)
79
80;; Oddball functions
81
82(setf (get '$max 'commutes-with-conjugate) t)
83(setf (get '$min 'commutes-with-conjugate) t)
84
85;; When a function has the conjugate-function property,
86;; use a non-generic function to conjugate it. Not done:
87;; conjugate-functions for all the inverse trigonometric
88;; functions.
89
90;; Trig like and hypergeometric like functions
91
92(setf (get '%log 'conjugate-function) 'conjugate-log)
93(setf (get 'mexpt 'conjugate-function) 'conjugate-mexpt)
94(setf (get '%asin 'conjugate-function) 'conjugate-asin)
95(setf (get '%acos 'conjugate-function) 'conjugate-acos)
96(setf (get '%atan 'conjugate-function) 'conjugate-atan)
97(setf (get '%atanh 'conjugate-function) 'conjugate-atanh)
98
99;;(setf (get '$asec 'conjugate-function) 'conjugate-asec)
100;;(setf (get '$acsc 'conjugate-function) 'conjugate-acsc)
101(setf (get '%bessel_j 'conjugate-function) 'conjugate-bessel-j)
102(setf (get '%bessel_y 'conjugate-function) 'conjugate-bessel-y)
103(setf (get '%bessel_i 'conjugate-function) 'conjugate-bessel-i)
104(setf (get '%bessel_k 'conjugate-function) 'conjugate-bessel-k)
105
106(setf (get '%hankel_1 'conjugate-function) 'conjugate-hankel-1)
107(setf (get '%hankel_2 'conjugate-function) 'conjugate-hankel-2)
108
109;; Other things:
110
111(setf (get '%sum 'conjugate-function) 'conjugate-sum)
112(setf (get '%product 'conjugate-function) 'conjugate-product)
113
114;; Return true iff Maxima can prove that z is not on the
115;; negative real axis.
116
117(defun off-negative-real-axisp (z)
118  (setq z (trisplit z))	          ; split into real and imaginary
119  (or (eq t (mnqp (cdr z) 0))     ; y #  0
120      (eq t (mgqp (car z) 0))))   ; x >= 0
121
122(defun on-negative-real-axisp (z)
123  (setq z (trisplit z))
124  (and (eq t (meqp (cdr z) 0))
125       (eq t (mgrp 0 (car z)))))
126
127(defun in-domain-of-asin (z)
128  (setq z (trisplit z))
129  (let ((x (car z)) (y (cdr z)))
130    (or (eq t (mgrp y 0))
131	(eq t (mgrp 0 y))
132	(and
133	 (eq t (mgrp x -1))
134	 (eq t (mgrp 1 x))))))
135
136;; Return conjugate(log(x)). Actually, x is a lisp list (x).
137
138(defun conjugate-log (x)
139  (setq x (car x))
140  (cond ((off-negative-real-axisp x)
141	 (take '(%log) (take '($conjugate) x)))
142	((on-negative-real-axisp x)
143	 (add (take '(%log) (neg x)) (mul -1 '$%i '$%pi)))
144	(t `(($conjugate simp) ((%log simp) ,x)))))
145
146;; Return conjugate(x^p), where e = (x, p). Suppose x isn't on the negative real axis.
147;; Then conjugate(x^p) == conjugate(exp(p * log(x))) == exp(conjugate(p) * conjugate(log(x)))
148;; == exp(conjugate(p) * log(conjugate(x)) = conjugate(x)^conjugate(p). Thus, when
149;; x is off the negative real axis, commute the conjugate with ^. Also if p is an integer
150;; ^ commutes with the conjugate.
151
152(defun conjugate-mexpt (e)
153  (let ((x (first e)) (p (second e)))
154    (if (or (off-negative-real-axisp x) ($featurep p '$integer))
155	(power (take '($conjugate) x) (take '($conjugate) p))
156      `(($conjugate simp) ,(power x p)))))
157
158(defun conjugate-sum (e)
159  (take '(%sum) (take '($conjugate) (first e)) (second e) (third e) (fourth e)))
160
161(defun conjugate-product (e)
162  (take '(%product) (take '($conjugate) (first e)) (second e) (third e) (fourth e)))
163
164(defun conjugate-asin (x)
165  (setq x (car x))
166  (if (in-domain-of-asin x) (take '(%asin) (take '($conjugate) x))
167    `(($conjugate simp) ((%asin) ,x))))
168
169(defun conjugate-acos (x)
170  (setq x (car x))
171  (if (in-domain-of-asin x) (take '(%acos) (take '($conjugate) x))
172    `(($conjugate simp) ((%acos) ,x))))
173
174(defun conjugate-atan (x)
175  (let ((xx))
176    (setq x (car x))
177    (setq xx (mul '$%i x))
178    (if (in-domain-of-asin xx)
179        (take '(%atan) (take '($conjugate) x))
180        `(($conjugate simp) ((%atan) ,x)))))
181
182;; atanh and asin are entire on the same set; see A&S Fig. 4.4 and 4.7.
183
184(defun conjugate-atanh (x)
185  (setq x (car x))
186  (if (in-domain-of-asin x) (take '(%atanh) (take '($conjugate) x))
187    `(($conjugate simp) ((%atanh) ,x))))
188
189;; Integer order Bessel functions are entire; thus they commute with the
190;; conjugate (Schwartz refection principle).  But non-integer order Bessel
191;; functions are not analytic along the negative real axis. Notice that A&S
192;; 9.1.40 isn't correct -- it says that the real order Bessel functions
193;; commute with the conjugate. Not true.
194
195(defun conjugate-bessel-j (z)
196  (let ((n (first z)) (x (second z)))
197    (if (or ($featurep n '$integer) (off-negative-real-axisp x))
198        (take '(%bessel_j) (take '($conjugate) n) (take '($conjugate) x))
199       `(($conjugate simp) ((%bessel_j simp) ,@z)))))
200
201(defun conjugate-bessel-y (z)
202  (let ((n (first z)) (x (second z)))
203    (if (off-negative-real-axisp x)
204        (take '(%bessel_y) (take '($conjugate) n) (take '($conjugate) x))
205       `(($conjugate simp) ((%bessel_y simp) ,@z)))))
206
207(defun conjugate-bessel-i (z)
208  (let ((n (first z)) (x (second z)))
209    (if (or ($featurep n '$integer) (off-negative-real-axisp x))
210        (take '(%bessel_i) (take '($conjugate) n) (take '($conjugate) x))
211       `(($conjugate simp) ((%bessel_i simp) ,@z)))))
212
213(defun conjugate-bessel-k (z)
214  (let ((n (first z)) (x (second z)))
215    (if (off-negative-real-axisp x)
216        (take '(%bessel_k) (take '($conjugate) n) (take '($conjugate) x))
217       `(($conjugate simp) ((%bessel_k simp) ,@z)))))
218
219(defun conjugate-hankel-1 (z)
220  (let ((n (first z)) (x (second z)))
221    (if (off-negative-real-axisp x)
222        (take '(%hankel_2) (take '($conjugate) n) (take '($conjugate) x))
223       `(($conjugate simp) ((%hankel_1 simp) ,@z)))))
224
225(defun conjugate-hankel-2 (z)
226  (let ((n (first z)) (x (second z)))
227    (if (off-negative-real-axisp x)
228        (take '(%hankel_1) (take '($conjugate) n) (take '($conjugate) x))
229       `(($conjugate simp) ((%hankel_2 simp) ,@z)))))
230
231;; When a function maps "everything" into the reals, put real-valued on the
232;; property list of the function name. This duplicates some knowledge that
233;; $rectform has. So it goes. The functions floor and ceiling also aren't
234;; defined off the real-axis. I suppose these functions could be given the
235;; real-valued property.
236
237(setf (get '%imagpart 'real-valued) t)
238(setf (get 'mabs 'real-valued) t)
239(setf (get '%realpart 'real-valued) t)
240(setf (get '%carg 'real-valued) t)
241
242;; manifestly-real-p isn't a great name, but it's OK. Since (manifestly-real-p '$inf) --> true
243;; it might be called manifestly-extended-real-p. A nonscalar isn't real.
244
245;; There might be some advantage to requiring that the subscripts to a $subvarp
246;; all be real.  Why? Well li[n] maps reals to reals when n is real, but li[n] does
247;; not map the reals to reals when n is nonreal.
248
249(defun manifestly-real-p (e)
250  (let (($inflag t))
251    (and ($mapatom e)
252	 (not (manifestly-pure-imaginary-p e))
253	 (not (manifestly-complex-p e))
254	 (not (manifestly-nonreal-p e))
255	 (or
256	  ($numberp e)
257	  (symbolp e)
258	  (and ($subvarp e) (manifestly-real-p ($op e)))))))
259
260(defun manifestly-pure-imaginary-p (e)
261  (let (($inflag t))
262    (or
263     (and ($mapatom e)
264	  (or
265	   (eq e '$%i)
266	   (and (symbolp e) (kindp e '$imaginary) (not ($nonscalarp e)))
267	   (and ($subvarp e) (manifestly-pure-imaginary-p ($op e)))))
268     ;; For now, let's use $csign on constant expressions only; once $csign improves,
269     ;; the ban on nonconstant expressions can be removed
270     (and ($constantp e) (eq '$imaginary ($csign e))))))
271
272;; Don't use (kindp e '$complex)!
273
274(defun manifestly-complex-p (e)
275  (let (($inflag t))
276    (or (and (symbolp e) (decl-complexp e) (not ($nonscalarp e)))
277	(eq e '$infinity)
278	(and ($subvarp e) (manifestly-complex-p ($op e))
279	     (not ($nonscalarp e))))))
280
281(defun manifestly-nonreal-p (e)
282  (and (symbolp e) (or (member e `($und $ind $zeroa $zerob t nil)) ($nonscalarp e))))
283
284;; For a subscripted function, conjugate always returns the conjugate noun-form.
285;; This could be repaired. For now, we don't have a scheme for conjugate(li[m](x)).
286
287;; We could make commutes_with_conjugate and maps_to_reals features. But I
288;; doubt it would get much use.
289
290(defun simp-conjugate (e f z)
291  (oneargcheck e)
292  (setq e (simpcheck (cadr e) z))	; simp and disrep if necessary
293  (cond ((complexp e) (conjugate e))    ; never happens, but might someday.
294	((manifestly-real-p e) e)
295	((manifestly-pure-imaginary-p e) (mul -1 e))
296	((manifestly-nonreal-p e) `(($conjugate simp) ,e))
297	(($mapatom e) `(($conjugate simp) ,e))
298	((op-equalp e '$conjugate) (car (margs e)))
299
300	((and (symbolp (mop e)) (get (mop e) 'real-valued)) e)
301
302	((and (symbolp (mop e)) (get (mop e) 'commutes-with-conjugate))
303	 (simplify (cons (list (mop e)) (mapcar #'(lambda (s) (take '($conjugate) s)) (margs e)))))
304
305	((setq f (and (symbolp (mop e)) (get (mop e) 'conjugate-function)))
306	 (funcall f (margs e)))
307
308	(t `(($conjugate simp) ,e))))
309