1 ;; Copyright 2005, 2006, 2020 by Barton Willis
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.
7 ;; This software has NO WARRANTY, not even the implied warranty of
8 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
12 (macsyma-module conjugate
)
14 ($put
'$conjugate
1 '$version
)
15 ;; Let's remove built-in symbols from list for user-defined properties.
16 (setq $props
(remove '$conjugate $props
))
18 (defprop $conjugate tex-postfix tex
)
19 (defprop $conjugate
("^\\star") texsym
)
20 (defprop $conjugate
160. tex-lbp
)
21 (defprop $conjugate simp-conjugate operators
)
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
))))
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.
37 ;; What important functions have I missed?
39 ;; (1) Arithmetic operators
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
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
)
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
)
65 (setf (get '%gamma
'commutes-with-conjugate
) t
)
66 (setf (get '$pochhammer
'commutes-with-conjugate
) t
)
70 (setf (get '$matrix
'commutes-with-conjugate
) t
)
71 (setf (get 'mlist
'commutes-with-conjugate
) t
)
72 (setf (get '$set
'commutes-with-conjugate
) t
)
76 (setf (get 'mequal
'commutes-with-conjugate
) t
)
77 (setf (get 'mnotequal
'commutes-with-conjugate
) t
)
78 (setf (get '%transpose
'commutes-with-conjugate
) t
)
82 (setf (get '$max
'commutes-with-conjugate
) t
)
83 (setf (get '$min
'commutes-with-conjugate
) t
)
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
90 ;; Trig like and hypergeometric like functions
92 (setf (get '%log
'conjugate-function
) 'conjugate-log
)
93 (setf (get '%plog
'conjugate-function
) 'conjugate-log
)
94 (setf (get 'mexpt
'conjugate-function
) 'conjugate-mexpt
)
95 (setf (get '%asin
'conjugate-function
) 'conjugate-asin
)
96 (setf (get '%acos
'conjugate-function
) 'conjugate-acos
)
97 (setf (get '%atan
'conjugate-function
) 'conjugate-atan
)
98 (setf (get '%atanh
'conjugate-function
) 'conjugate-atanh
)
100 ;;(setf (get '$asec 'conjugate-function) 'conjugate-asec)
101 ;;(setf (get '$acsc 'conjugate-function) 'conjugate-acsc)
102 (setf (get '%bessel_j
'conjugate-function
) 'conjugate-bessel-j
)
103 (setf (get '%bessel_y
'conjugate-function
) 'conjugate-bessel-y
)
104 (setf (get '%bessel_i
'conjugate-function
) 'conjugate-bessel-i
)
105 (setf (get '%bessel_k
'conjugate-function
) 'conjugate-bessel-k
)
107 (setf (get '%hankel_1
'conjugate-function
) 'conjugate-hankel-1
)
108 (setf (get '%hankel_2
'conjugate-function
) 'conjugate-hankel-2
)
109 (setf (get '%log_gamma
'conjugate-function
) 'conjugate-log-gamma
)
111 ;; conjugate of polylogarithm li & psi
112 (setf (get '$li
'conjugate-function
) 'conjugate-li
)
113 (setf (get '$psi
'conjugate-function
) 'conjugate-psi
)
116 (setf (get '%sum
'conjugate-function
) 'conjugate-sum
)
117 (setf (get '%product
'conjugate-function
) 'conjugate-product
)
119 ;; Return true iff Maxima can prove that z is not on the
120 ;; negative real axis.
122 (defun off-negative-real-axisp (z)
123 (setq z
(trisplit z
)) ; split into real and imaginary
124 (or (eql t
(mnqp (cdr z
) 0)) ; y # 0
125 (eql t
(mgqp (car z
) 0)))) ; x >= 0
127 (defun on-negative-real-axisp (z)
128 (setq z
(trisplit z
))
129 (and (eql t
(meqp (cdr z
) 0))
130 (eql t
(mgrp 0 (car z
)))))
132 (defun in-domain-of-asin (z)
133 (setq z
(trisplit z
))
134 (let ((x (car z
)) (y (cdr z
)))
135 (or (eql t
(mgrp y
0))
139 (eql t
(mgrp 1 x
))))))
141 ;; Return conjugate(log(x)). Actually, x is a lisp list (x).
143 (defun conjugate-log (x)
145 (cond ((off-negative-real-axisp x
)
146 (take '(%log
) (take '($conjugate
) x
)))
147 ((on-negative-real-axisp x
)
148 (add (take '(%log
) (neg x
)) (mul -
1 '$%i
'$%pi
)))
149 (t `(($conjugate simp
) ((%log simp
) ,x
)))))
151 ;; Return conjugate(x^p), where e = (x, p). Suppose x isn't on the negative real axis.
152 ;; Then conjugate(x^p) == conjugate(exp(p * log(x))) == exp(conjugate(p) * conjugate(log(x)))
153 ;; == exp(conjugate(p) * log(conjugate(x)) = conjugate(x)^conjugate(p). Thus, when
154 ;; x is off the negative real axis, commute the conjugate with ^. Also if p is an integer
155 ;; ^ commutes with the conjugate.
157 (defun conjugate-mexpt (e)
158 (let ((x (first e
)) (p (second e
)))
159 (if (or (off-negative-real-axisp x
) ($featurep p
'$integer
))
160 (power (take '($conjugate
) x
) (take '($conjugate
) p
))
161 `(($conjugate simp
) ,(power x p
)))))
163 (defun conjugate-sum (e)
164 (take '(%sum
) (take '($conjugate
) (first e
)) (second e
) (third e
) (fourth e
)))
166 (defun conjugate-product (e)
167 (take '(%product
) (take '($conjugate
) (first e
)) (second e
) (third e
) (fourth e
)))
169 (defun conjugate-asin (x)
171 (if (in-domain-of-asin x
) (take '(%asin
) (take '($conjugate
) x
))
172 `(($conjugate simp
) ((%asin
) ,x
))))
174 (defun conjugate-acos (x)
176 (if (in-domain-of-asin x
) (take '(%acos
) (take '($conjugate
) x
))
177 `(($conjugate simp
) ((%acos
) ,x
))))
179 (defun conjugate-atan (x)
182 (setq xx
(mul '$%i x
))
183 (if (in-domain-of-asin xx
)
184 (take '(%atan
) (take '($conjugate
) x
))
185 `(($conjugate simp
) ((%atan
) ,x
)))))
187 ;; atanh and asin are entire on the same set; see A&S Fig. 4.4 and 4.7.
189 (defun conjugate-atanh (x)
191 (if (in-domain-of-asin x
) (take '(%atanh
) (take '($conjugate
) x
))
192 `(($conjugate simp
) ((%atanh
) ,x
))))
194 ;; Integer order Bessel functions are entire; thus they commute with the
195 ;; conjugate (Schwartz refection principle). But non-integer order Bessel
196 ;; functions are not analytic along the negative real axis. Notice that A&S
197 ;; 9.1.40 isn't correct -- it says that the real order Bessel functions
198 ;; commute with the conjugate. Not true.
200 (defun conjugate-bessel-j (z)
201 (let ((n (first z
)) (x (second z
)))
202 (if (or ($featurep n
'$integer
) (off-negative-real-axisp x
))
203 (take '(%bessel_j
) (take '($conjugate
) n
) (take '($conjugate
) x
))
204 `(($conjugate simp
) ((%bessel_j simp
) ,@z
)))))
206 (defun conjugate-bessel-y (z)
207 (let ((n (first z
)) (x (second z
)))
208 (if (off-negative-real-axisp x
)
209 (take '(%bessel_y
) (take '($conjugate
) n
) (take '($conjugate
) x
))
210 `(($conjugate simp
) ((%bessel_y simp
) ,@z
)))))
212 (defun conjugate-bessel-i (z)
213 (let ((n (first z
)) (x (second z
)))
214 (if (or ($featurep n
'$integer
) (off-negative-real-axisp x
))
215 (take '(%bessel_i
) (take '($conjugate
) n
) (take '($conjugate
) x
))
216 `(($conjugate simp
) ((%bessel_i simp
) ,@z
)))))
218 (defun conjugate-bessel-k (z)
219 (let ((n (first z
)) (x (second z
)))
220 (if (off-negative-real-axisp x
)
221 (take '(%bessel_k
) (take '($conjugate
) n
) (take '($conjugate
) x
))
222 `(($conjugate simp
) ((%bessel_k simp
) ,@z
)))))
224 (defun conjugate-hankel-1 (z)
225 (let ((n (first z
)) (x (second z
)))
226 (if (off-negative-real-axisp x
)
227 (take '(%hankel_2
) (take '($conjugate
) n
) (take '($conjugate
) x
))
228 `(($conjugate simp
) ((%hankel_1 simp
) ,@z
)))))
230 (defun conjugate-hankel-2 (z)
231 (let ((n (first z
)) (x (second z
)))
232 (if (off-negative-real-axisp x
)
233 (take '(%hankel_1
) (take '($conjugate
) n
) (take '($conjugate
) x
))
234 `(($conjugate simp
) ((%hankel_2 simp
) ,@z
)))))
236 (defun conjugate-log-gamma (z)
238 (if (off-negative-real-axisp z
)
239 (take '(%log_gamma
) (take '($conjugate
) z
))
240 `(($conjugate simp
) ((%log_gamma simp
) ,z
))))
242 ;; conjugate of polylogarithm li[s](x), where z = (s,x). We have li[s](x) = x+x^2/2^s+x^3/3^s+...
243 ;; Since for all integers k, we have conjugate(x^k/k^s) = conjugate(x)^k/k^conjugate(s), we
244 ;; commute conjugate with li.
245 (defun conjugate-li (z)
246 (let ((s (take '($conjugate
) (first z
))) (x (take '($conjugate
) (second z
))))
247 (take '(mqapply) `(($li array
) ,s
) x
)))
249 (defun conjugate-psi (z)
250 (let ((s (take '($conjugate
) (first z
))) (x (take '($conjugate
) (second z
))))
251 (take '(mqapply) `(($psi array
) ,s
) x
)))
253 ;; When a function maps "everything" into the reals, put real-valued on the
254 ;; property list of the function name. This duplicates some knowledge that
255 ;; $rectform has. So it goes.
257 (setf (get '%imagpart
'real-valued
) t
)
258 (setf (get 'mabs
'real-valued
) t
)
259 (setf (get '%realpart
'real-valued
) t
)
260 (setf (get '%carg
'real-valued
) t
)
261 (setf (get '$ceiling
'real-valued
) t
)
262 (setf (get '$floor
'real-valued
) t
)
263 (setf (get '$mod
'real-valued
) t
)
264 (setf (get '$unit_step
'real-valued
) t
)
265 (setf (get '$charfun
'real-valued
) t
)
268 ;; The function manifestly-real-p makes some effort to determine if its input is
271 ;; manifestly-real-p isn't a great name, but it's OK. Since (manifestly-real-p '$inf) --> true
272 ;; it might be called manifestly-extended-real-p. A nonscalar isn't real.
274 ;; There might be some advantage to requiring that the subscripts to a $subvarp
275 ;; all be real. Why? Well li[n] maps reals to reals when n is real, but li[n] does
276 ;; not map the reals to reals when n is nonreal.
278 (defun manifestly-real-p (e)
283 (not (manifestly-pure-imaginary-p e
))
284 (not (manifestly-complex-p e
))
285 (not (manifestly-nonreal-p e
)))
286 (and (consp e
) (consp (car e
)) (get (caar e
) 'real-valued
)) ;F(xxx), where F is declared real-valued
287 (and ($subvarp e
) (manifestly-real-p ($op e
)))))) ;F[n], where F is declared real-valued
289 ;; The function manifestly-pure-imaginary-p makes some effort to determine if its input is
292 (defun manifestly-pure-imaginary-p (e)
298 (and (symbolp e
) (kindp e
'$imaginary
) (not ($nonscalarp e
)))
299 (and ($subvarp e
) (manifestly-pure-imaginary-p ($op e
)))))
300 ;; For now, let's use $csign on constant expressions only; once $csign improves,
301 ;; the ban on nonconstant expressions can be removed
302 (and ($constantp e
) (eql '$imaginary
($csign e
))))))
304 ;; Don't use (kindp e '$complex)!
306 (defun manifestly-complex-p (e)
308 (or (and (symbolp e
) (decl-complexp e
) (not ($nonscalarp e
)))
310 (and ($subvarp e
) (manifestly-complex-p ($op e
)) (not ($nonscalarp e
))))))
312 (defun manifestly-nonreal-p (e)
313 (and (symbolp e
) (or (member e
`($und $ind $zeroa $zerob t nil
)) ($nonscalarp e
))))
315 ;; We could make commutes_with_conjugate and maps_to_reals features. But I
316 ;; doubt it would get much use.
318 (defun simp-conjugate (e f z
)
320 (setq e
(simpcheck (cadr e
) z
)) ; simp and disrep if necessary
321 (cond ((complexp e
) (conjugate e
)) ; never happens, but might someday.
322 ((manifestly-real-p e
) e
)
323 ((manifestly-pure-imaginary-p e
) (mul -
1 e
))
324 ((manifestly-nonreal-p e
) `(($conjugate simp
) ,e
))
325 (($mapatom e
) `(($conjugate simp
) ,e
))
326 ((op-equalp e
'$conjugate
) (car (margs e
)))
328 ((and (symbolp (mop e
)) (get (mop e
) 'real-valued
)) e
)
330 ((and (symbolp (mop e
)) (get (mop e
) 'commutes-with-conjugate
))
331 (simplify (cons (list (mop e
)) (mapcar #'(lambda (s) (take '($conjugate
) s
)) (margs e
)))))
333 ((setq f
(and (symbolp (mop e
)) (get (mop e
) 'conjugate-function
)))
334 (funcall f
(margs e
)))
336 ;;subscripted functions
337 ((setq f
(and ($subvarp
(mop e
)) (get (caar (mop e
)) 'conjugate-function
)))
339 (funcall f
(append (margs (mop e
)) (margs e
))))
341 (t `(($conjugate simp
) ,e
))))
