1 ;; Copyright 2005, 2006, 2020, 2021 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
)
23 ;; Maybe $conjugate should have a msimpind property. But with some Maxima versions,
24 ;; kill(conjugate) eliminates the msimpind property; after that, conjugate gives rubbish.
25 ;; Until this is resolved, $conjugate doesn't have a msimpind property.
28 (:load-toplevel
:execute
)
29 (let (($context
'$global
) (context '$global
))
30 (meval '(($declare
) $conjugate $complex
))
31 ;; Let's remove built-in symbols from list for user-defined properties.
32 (setq $props
(remove '$conjugate $props
))))
34 ;; When a function commutes with the conjugate, give the function the
35 ;; commutes-with-conjugate property. The log function commutes with
36 ;; the conjugate on all of C except on the negative real axis. Thus
37 ;; log does not get the commutes-with-conjugate property. Instead,
38 ;; log gets the conjugate-function property.
40 ;; What important functions have I missed?
42 ;; (1) Arithmetic operators
44 (setf (get 'mplus
'commutes-with-conjugate
) t
)
45 (setf (get 'mtimes
'commutes-with-conjugate
) t
)
46 ;(setf (get 'mnctimes 'commutes-with-conjugate) t) ;; generally I think users will want this
47 (setf (get '%signum
'commutes-with-conjugate
) t
) ;; x=/=0, conjugate(signum(x)) = conjugate(x/abs(x)) = signum(conjugate(x))
48 ;; Trig-like functions and other such functions
50 (setf (get '%cosh
'commutes-with-conjugate
) t
)
51 (setf (get '%sinh
'commutes-with-conjugate
) t
)
52 (setf (get '%tanh
'commutes-with-conjugate
) t
)
53 (setf (get '%sech
'commutes-with-conjugate
) t
)
54 (setf (get '%csch
'commutes-with-conjugate
) t
)
55 (setf (get '%coth
'commutes-with-conjugate
) t
)
56 (setf (get '%cos
'commutes-with-conjugate
) t
)
57 (setf (get '%sin
'commutes-with-conjugate
) t
)
58 (setf (get '%tan
'commutes-with-conjugate
) t
)
59 (setf (get '%sec
'commutes-with-conjugate
) t
)
60 (setf (get '%csc
'commutes-with-conjugate
) t
)
61 (setf (get '%cot
'commutes-with-conjugate
) t
)
62 (setf (get '$atan2
'commutes-with-conjugate
) t
)
64 (setf (get '%jacobi_cn
'commutes-with-conjugate
) t
)
65 (setf (get '%jacobi_sn
'commutes-with-conjugate
) t
)
66 (setf (get '%jacobi_dn
'commutes-with-conjugate
) t
)
68 (setf (get '%gamma
'commutes-with-conjugate
) t
)
69 (setf (get '$pochhammer
'commutes-with-conjugate
) t
)
73 (setf (get '$matrix
'commutes-with-conjugate
) t
)
74 (setf (get 'mlist
'commutes-with-conjugate
) t
)
75 (setf (get '$set
'commutes-with-conjugate
) t
)
79 (setf (get 'mequal
'commutes-with-conjugate
) t
)
80 (setf (get 'mnotequal
'commutes-with-conjugate
) t
)
81 (setf (get '%transpose
'commutes-with-conjugate
) t
)
85 (setf (get '$max
'commutes-with-conjugate
) t
)
86 (setf (get '$min
'commutes-with-conjugate
) t
)
88 ;; When a function has the conjugate-function property, use a non-generic function to conjugate it.
89 ;; The argument to a conjugate function for an operator op is the CL list of arguments to op. For
90 ;; example, the conjugate function for log gets the argument for log, not the expression log(x).
91 ;; It would be a bit more efficient if a conjugate function received the full expression--that
92 ;; way for a pure nounform return (for example, return conjugate(log(x))), a conjugate function
93 ;; would not not need to apply the operator to the argument to the conjugate function, instead it
94 ;; could simply paste ($conjugate simp) onto the expression.
96 ;; Not done: conjugate-functions for all the inverse trigonometric functions.
98 ;; Trig like and hypergeometric like functions
100 (setf (get '%log
'conjugate-function
) 'conjugate-log
)
101 (setf (get '%plog
'conjugate-function
) 'conjugate-plog
)
102 (setf (get 'mexpt
'conjugate-function
) 'conjugate-mexpt
)
103 (setf (get '%asin
'conjugate-function
) 'conjugate-asin
)
104 (setf (get '%acos
'conjugate-function
) 'conjugate-acos
)
105 (setf (get '%atan
'conjugate-function
) 'conjugate-atan
)
106 (setf (get '%atanh
'conjugate-function
) 'conjugate-atanh
)
107 (setf (get '%asec
'conjugate-function
) 'conjugate-asec
)
108 (setf (get '%acsc
'conjugate-function
) 'conjugate-acsc
)
110 (setf (get '%bessel_j
'conjugate-function
) 'conjugate-bessel-j
)
111 (setf (get '%bessel_y
'conjugate-function
) 'conjugate-bessel-y
)
112 (setf (get '%bessel_i
'conjugate-function
) 'conjugate-bessel-i
)
113 (setf (get '%bessel_k
'conjugate-function
) 'conjugate-bessel-k
)
115 (setf (get '%hankel_1
'conjugate-function
) 'conjugate-hankel-1
)
116 (setf (get '%hankel_2
'conjugate-function
) 'conjugate-hankel-2
)
117 (setf (get '%log_gamma
'conjugate-function
) 'conjugate-log-gamma
)
119 ;; conjugate of polylogarithm li & psi
120 (setf (get '$li
'conjugate-function
) 'conjugate-li
)
121 (setf (get '$psi
'conjugate-function
) 'conjugate-psi
)
124 (setf (get '%sum
'conjugate-function
) 'conjugate-sum
)
125 (setf (get '%product
'conjugate-function
) 'conjugate-product
)
127 ;; Return true iff Maxima can prove that z is not on the
128 ;; negative real axis.
130 (defun off-negative-real-axisp (z)
131 (setq z
(trisplit z
)) ; split into real and imaginary
132 (or (eql t
(mnqp (cdr z
) 0)) ; y # 0
133 (eql t
(mgqp (car z
) 0)))) ; x >= 0
135 (defun on-negative-real-axisp (z)
136 (setq z
(trisplit z
))
137 (and (eql t
(meqp (cdr z
) 0))
138 (eql t
(mgrp 0 (car z
)))))
140 (defun off-negative-one-to-onep (z)
141 (setq z
(trisplit z
)) ; split z into real and imaginary parts
143 (eq t
(mnqp (cdr z
) 0)) ; y # 0
144 (eq t
(mgrp (car z
) 1)) ; x > 1
145 (eq t
(mgrp -
1 (car z
))))) ; -1 > x
147 (defun in-domain-of-asin (z)
148 (setq z
(trisplit z
)) ; split z into real and imaginary parts
149 (let ((x (car z
)) (y (cdr z
))) ;z = x+%i*y
151 (eq t
(mnqp y
0)) ; y # 0
153 (eq t
(mgrp x -
1)) ; x > -1
154 (eq t
(mgrp 1 x
)))))) ; x < 1
156 ;; Return conjugate(log(x)). Actually, x is a lisp list (x).
158 (defun conjugate-log (x)
160 (cond ((off-negative-real-axisp x
)
161 (take '(%log
) (take '($conjugate
) x
)))
162 ((on-negative-real-axisp x
)
163 (add (take '(%log
) (neg x
)) (mul -
1 '$%i
'$%pi
)))
164 (t (list '($conjugate simp
) (take '(%log
) x
)))))
167 ;; Return conjugate(plog(x)); again, x is the CL list (x).
168 (defun conjugate-plog (x)
170 (cond ((off-negative-real-axisp x
)
171 (take '(%plog
) (take '($conjugate
) x
)))
172 ((on-negative-real-axisp x
)
173 (add (take '(%plog
) (neg x
)) (mul -
1 '$%i
'$%pi
)))
174 (t (list '($conjugate simp
) (take '(%plog
) x
)))))
176 ;; Return conjugate(x^p), where e = (x, p). Suppose x isn't on the negative real axis.
177 ;; Then conjugate(x^p) == conjugate(exp(p * log(x))) == exp(conjugate(p) * conjugate(log(x)))
178 ;; == exp(conjugate(p) * log(conjugate(x)) = conjugate(x)^conjugate(p). Thus, when
179 ;; x is off the negative real axis, commute the conjugate with ^. Also if p is an integer
180 ;; ^ commutes with the conjugate.
182 ;; We don't need to call $ratdisrep before checking if p is a declared integer--the
183 ;; simpcheck at the top level of simp-conjugate does that for us. So we can call
184 ;; maxima-integerp on p instead of using $featurep.
186 ;; The rule that is commented out is, I think, correct, but I'm not sure how useful it is and
187 ;; the testsuite plus the share testsuite never use this rule. For now, let's keep
190 ;; Running the testsuite plus the share testsuite calls conjugate-mexpt 63,441
191 ;; times. This is far more times than all the other conjugate functions. Of these
192 ;; calls, the exponent is an integer 63,374 times. So for efficiency, we check
193 ;; (maxima-integerp p) first.
195 ;; The case of a nounform return only happens 9 times. For the nounform return, the power has
196 ;; been simplified at the higher level. So at least for running the testsuite, we shouldn't
197 ;; worry all that much about re-simplifying the power for the nounform return.
199 (defun conjugate-mexpt (e)
200 (let ((x (first e
)) (p (second e
)))
201 (cond ((or (maxima-integerp p
) (off-negative-real-axisp x
))
202 (power (take '($conjugate
) x
) (take '($conjugate
) p
)))
203 ;((on-negative-real-axisp x) ;conjugate(x^p) = exp(-%i %pi conjugate(p)) (-x)^p
204 ; (setq p (take '($conjugate) p))
205 ; (mul (power '$%e (mul -1 '$%i '$%pi p)) (power (mul -1 x) p)))
207 (list '($conjugate simp
) (power x p
))))))
209 (defun conjugate-sum (e)
210 (if (and ($featurep
(third e
) '$real
) ($featurep
(fourth e
) '$real
))
211 (take '(%sum
) (take '($conjugate
) (first e
)) (second e
) (third e
) (fourth e
))
212 (list '($conjugate simp
) (simplifya (cons '(%sum
) e
) t
))))
214 (defun conjugate-product (e)
215 (if (and ($featurep
(third e
) '$real
) ($featurep
(fourth e
) '$real
))
216 (take '(%product
) (take '($conjugate
) (first e
)) (second e
) (third e
) (fourth e
))
217 (list '($conjugate simp
) (simplifya (cons '(%product
) e
) t
))))
219 (defun conjugate-asin (x)
221 (if (in-domain-of-asin x
) (take '(%asin
) (take '($conjugate
) x
))
222 (list '($conjugate simp
) (take '(%asin
) x
))))
224 (defun conjugate-acos (x)
226 (if (in-domain-of-asin x
) (take '(%acos
) (take '($conjugate
) x
))
227 (list '($conjugate simp
) (take '(%acos
) x
))))
229 (defun conjugate-acsc (x)
231 (if (off-negative-one-to-onep x
) (take '(%acsc
) (take '($conjugate
) x
))
232 (list '($conjugate simp
) (take '(%acsc
) x
))))
234 (defun conjugate-asec (x)
236 (if (off-negative-one-to-onep x
) (take '(%asec
) (take '($conjugate
) x
))
237 (list '($conjugate simp
) (take '(%asec
) x
))))
239 (defun conjugate-atan (x)
242 (setq xx
(mul '$%i x
))
243 (if (in-domain-of-asin xx
)
244 (take '(%atan
) (take '($conjugate
) x
))
245 (list '($conjugate simp
) (take '(%atan
) x
)))))
247 ;; atanh and asin are entire on the same set; DLMF http://dlmf.nist.gov/4.37.F1 and
248 ;; http://dlmf.nist.gov/4.23.F1
250 (defun conjugate-atanh (x)
252 (if (in-domain-of-asin x
) (take '(%atanh
) (take '($conjugate
) x
))
253 (list '($conjugate simp
) (take '(%atanh
) x
))))
255 ;; Integer order Bessel functions are entire; thus they commute with the
256 ;; conjugate (Schwartz refection principle). But non-integer order Bessel
257 ;; functions are not analytic along the negative real axis. Notice that DLMF
258 ;; http://dlmf.nist.gov/10.11.E9 isn't correct; we have, for example
259 ;; conjugate(bessel_j(1/2,-1)) =/= bessel_j(1/2,conjugate(-1))
261 (defun conjugate-bessel-j (z)
262 (let ((n (first z
)) (x (second z
)))
263 (if (or ($featurep n
'$integer
) (off-negative-real-axisp x
))
264 (take '(%bessel_j
) (take '($conjugate
) n
) (take '($conjugate
) x
))
265 (list '($conjugate simp
) (simplifya (cons '(%bessel_j
) z
) t
)))))
267 (defun conjugate-bessel-y (z)
268 (let ((n (first z
)) (x (second z
)))
269 (if (off-negative-real-axisp x
)
270 (take '(%bessel_y
) (take '($conjugate
) n
) (take '($conjugate
) x
))
271 (list '($conjugate simp
) (simplifya (cons '(%bessel_y
) z
) t
)))))
273 (defun conjugate-bessel-i (z)
274 (let ((n (first z
)) (x (second z
)))
275 (if (or ($featurep n
'$integer
) (off-negative-real-axisp x
))
276 (take '(%bessel_i
) (take '($conjugate
) n
) (take '($conjugate
) x
))
277 (list '($conjugate simp
) (simplifya (cons '(%bessel_i
) z
) t
)))))
279 (defun conjugate-bessel-k (z)
280 (let ((n (first z
)) (x (second z
)))
281 (if (off-negative-real-axisp x
)
282 (take '(%bessel_k
) (take '($conjugate
) n
) (take '($conjugate
) x
))
283 (list '($conjugate simp
) (simplifya (cons '(%bessel_k
) z
) t
)))))
285 (defun conjugate-hankel-1 (z)
286 (let ((n (first z
)) (x (second z
)))
287 (if (off-negative-real-axisp x
)
288 (take '(%hankel_2
) (take '($conjugate
) n
) (take '($conjugate
) x
))
289 (list '($conjugate simp
) (simplifya (cons '(%hankel_1
) z
) t
)))))
291 (defun conjugate-hankel-2 (z)
292 (let ((n (first z
)) (x (second z
)))
293 (if (off-negative-real-axisp x
)
294 (take '(%hankel_1
) (take '($conjugate
) n
) (take '($conjugate
) x
))
295 (list '($conjugate simp
) (simplifya (cons '(%hankel_2
) z
) t
)))))
297 (defun conjugate-log-gamma (z)
299 (if (off-negative-real-axisp z
)
300 (take '(%log_gamma
) (take '($conjugate
) z
))
301 (list '($conjugate simp
) (take '(%log_gamma
) z
))))
303 ;; conjugate of polylogarithm li[s](x), where z = (s,x). We have li[s](x) = x+x^2/2^s+x^3/3^s+...
304 ;; Since for all integers k, we have conjugate(x^k/k^s) = conjugate(x)^k/k^conjugate(s), we
305 ;; commute conjugate with li.
306 (defun conjugate-li (z)
307 (let ((s (take '($conjugate
) (first z
))) (x (take '($conjugate
) (second z
))))
308 (take '(mqapply) `(($li array
) ,s
) x
)))
310 (defun conjugate-psi (z)
311 (let ((s (take '($conjugate
) (first z
))) (x (take '($conjugate
) (second z
))))
312 (take '(mqapply) `(($psi array
) ,s
) x
)))
314 ;; When all derivative variables & orders are real, commute the derivative with
316 (defun conjugate-derivative (z)
317 (cond ((every #'manifestly-real-p
(cdr z
))
318 (setq z
(cons (take '($conjugate
) (car z
)) (cdr z
)))
319 (simplifya (cons (list '%derivative
) z
) t
))
321 (list '($conjugate simp
) (simplifya (cons (list '%derivative
) z
) t
)))))
323 (setf (get '%derivative
'conjugate-function
) 'conjugate-derivative
)
325 ;; When a function maps "everything" into the reals, put real-valued on the
326 ;; property list of the function name. This duplicates some knowledge that
327 ;; $rectform has. So it goes.
329 (setf (get '%imagpart
'real-valued
) t
)
330 (setf (get 'mabs
'real-valued
) t
)
331 (setf (get '%realpart
'real-valued
) t
)
332 (setf (get '%carg
'real-valued
) t
)
333 (setf (get '$ceiling
'real-valued
) t
)
334 (setf (get '$floor
'real-valued
) t
)
335 (setf (get '$mod
'real-valued
) t
)
336 (setf (get '$unit_step
'real-valued
) t
)
337 (setf (get '$charfun
'real-valued
) t
)
340 ;; The function manifestly-real-p makes some effort to determine if its input is
343 ;; manifestly-real-p isn't a great name, but it's OK. Since (manifestly-real-p '$inf) --> true
344 ;; it might be called manifestly-extended-real-p. A nonscalar isn't real.
346 ;; There might be some advantage to requiring that the subscripts to a $subvarp
347 ;; all be real. Why? Well li[n] maps reals to reals when n is real, but li[n] does
348 ;; not map the reals to reals when n is nonreal.
350 (defun manifestly-real-p (e)
355 (not (manifestly-pure-imaginary-p e
))
356 (not (manifestly-complex-p e
))
357 (not (manifestly-nonreal-p e
)))
358 (and (consp e
) (consp (car e
)) (get (caar e
) 'real-valued
)) ;F(xxx), where F is declared real-valued
359 (and ($subvarp e
) (manifestly-real-p ($op e
)))))) ;F[n], where F is declared real-valued
361 ;; The function manifestly-pure-imaginary-p makes some effort to determine if its input is
364 (defun manifestly-pure-imaginary-p (e)
370 (and (symbolp e
) (kindp e
'$imaginary
) (not ($nonscalarp e
)))
371 (and ($subvarp e
) (manifestly-pure-imaginary-p ($op e
)))))
372 ;; For now, let's use $csign on constant expressions only; once $csign improves,
373 ;; the ban on nonconstant expressions can be removed.
374 (and ($constantp e
) (not (eq '$und e
)) (not (eq '$ind e
)) (eq '$imaginary
($csign e
))))))
376 ;; Don't use (kindp e '$complex)!
378 (defun manifestly-complex-p (e)
380 (or (and (symbolp e
) (decl-complexp e
) (not ($nonscalarp e
)))
382 (and ($subvarp e
) (manifestly-complex-p ($op e
)) (not ($nonscalarp e
))))))
384 (defun manifestly-nonreal-p (e)
385 (and (symbolp e
) (or (member e
`($und $ind t nil
)) ($nonscalarp e
))))
387 ;; We could make commutes_with_conjugate and maps_to_reals features. But I
388 ;; doubt it would get much use.
390 (defun simp-conjugate (e f z
)
392 (setq e
(simpcheck (cadr e
) z
)) ; simp and disrep if necessary
394 (cond ((complexp e
) (conjugate e
)) ; never happens, but might someday.
395 ((manifestly-real-p e
) e
)
396 ((manifestly-pure-imaginary-p e
) (mul -
1 e
))
397 ((or (manifestly-nonreal-p e
) ($mapatom e
))
398 (list '($conjugate simp
) e
))
400 ((op-equalp e
'$conjugate
) (car (margs e
)))
402 ((and (symbolp (mop e
)) (get (mop e
) 'real-valued
)) e
)
404 ((and (symbolp (mop e
)) (get (mop e
) 'commutes-with-conjugate
))
405 (simplify (cons (list (mop e
)) (mapcar #'(lambda (s) (take '($conjugate
) s
)) (margs e
)))))
407 ((setq f
(and (symbolp (mop e
)) (get (mop e
) 'conjugate-function
)))
408 (funcall f
(margs e
)))
410 ;;subscripted functions
411 ((setq f
(and ($subvarp
(mop e
)) (get (caar (mop e
)) 'conjugate-function
)))
412 (funcall f
(append (margs (mop e
)) (margs e
))))
415 (list '($conjugate simp
) e
))))
