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 (setf (get 'conjugate-superscript-star
'tex-lbp
) (tex-lbp 'mexpt
))
19 (setf (get 'conjugate-superscript-star
'tex-rbp
) (tex-rbp 'mexpt
))
21 (defun tex-conjugate (x l r
)
22 ;; Punt to TeX output for MEXPT, but defeat the special case for
23 ;; powers of trig functions, e.g. sin(x)^2 which is set as sin^2 x,
24 ;; by supplying an expression which has a different operator
25 ;; (namely CONJUGATE-SUPERSCRIPT-STAR) instead of MEXPT.
26 (tex-mexpt `((conjugate-superscript-star) ,(second x
) "\\ast") l r
))
28 (defprop $conjugate tex-conjugate tex
)
29 (defprop $conjugate
121. tex-lbp
)
30 (defprop $conjugate
120. tex-rbp
)
32 (defprop $conjugate simp-conjugate operators
)
34 ;; Maybe $conjugate should have a msimpind property. But with some Maxima versions,
35 ;; kill(conjugate) eliminates the msimpind property; after that, conjugate gives rubbish.
36 ;; Until this is resolved, $conjugate doesn't have a msimpind property.
39 (:load-toplevel
:execute
)
40 (let (($context
'$global
) (context '$global
))
41 (meval '(($declare
) $conjugate $complex
))
42 ;; Let's remove built-in symbols from list for user-defined properties.
43 (setq $props
(remove '$conjugate $props
))))
45 ;; When a function commutes with the conjugate, give the function the
46 ;; commutes-with-conjugate property. The log function commutes with
47 ;; the conjugate on all of C except on the negative real axis. Thus
48 ;; log does not get the commutes-with-conjugate property. Instead,
49 ;; log gets the conjugate-function property.
51 ;; What important functions have I missed?
53 ;; (1) Arithmetic operators
55 (setf (get 'mplus
'commutes-with-conjugate
) t
)
56 (setf (get 'mtimes
'commutes-with-conjugate
) t
)
57 ;(setf (get 'mnctimes 'commutes-with-conjugate) t) ;; generally I think users will want this
58 (setf (get '%signum
'commutes-with-conjugate
) t
) ;; x=/=0, conjugate(signum(x)) = conjugate(x/abs(x)) = signum(conjugate(x))
59 ;; Trig-like functions and other such functions
61 (setf (get '%cosh
'commutes-with-conjugate
) t
)
62 (setf (get '%sinh
'commutes-with-conjugate
) t
)
63 (setf (get '%tanh
'commutes-with-conjugate
) t
)
64 (setf (get '%sech
'commutes-with-conjugate
) t
)
65 (setf (get '%csch
'commutes-with-conjugate
) t
)
66 (setf (get '%coth
'commutes-with-conjugate
) t
)
67 (setf (get '%cos
'commutes-with-conjugate
) t
)
68 (setf (get '%sin
'commutes-with-conjugate
) t
)
69 (setf (get '%tan
'commutes-with-conjugate
) t
)
70 (setf (get '%sec
'commutes-with-conjugate
) t
)
71 (setf (get '%csc
'commutes-with-conjugate
) t
)
72 (setf (get '%cot
'commutes-with-conjugate
) t
)
73 (setf (get '%atan2
'commutes-with-conjugate
) t
)
75 (setf (get '%jacobi_cn
'commutes-with-conjugate
) t
)
76 (setf (get '%jacobi_sn
'commutes-with-conjugate
) t
)
77 (setf (get '%jacobi_dn
'commutes-with-conjugate
) t
)
79 (setf (get '%gamma
'commutes-with-conjugate
) t
)
80 (setf (get '$pochhammer
'commutes-with-conjugate
) t
)
84 (setf (get '$matrix
'commutes-with-conjugate
) t
)
85 (setf (get 'mlist
'commutes-with-conjugate
) t
)
86 (setf (get '$set
'commutes-with-conjugate
) t
)
90 (setf (get 'mequal
'commutes-with-conjugate
) t
)
91 (setf (get 'mnotequal
'commutes-with-conjugate
) t
)
92 (setf (get '%transpose
'commutes-with-conjugate
) t
)
96 (setf (get '$max
'commutes-with-conjugate
) t
)
97 (setf (get '$min
'commutes-with-conjugate
) t
)
99 ;; When a function has the conjugate-function property, use a non-generic function to conjugate it.
100 ;; The argument to a conjugate function for an operator op is the CL list of arguments to op. For
101 ;; example, the conjugate function for log gets the argument for log, not the expression log(x).
102 ;; It would be a bit more efficient if a conjugate function received the full expression--that
103 ;; way for a pure nounform return (for example, return conjugate(log(x))), a conjugate function
104 ;; would not not need to apply the operator to the argument to the conjugate function, instead it
105 ;; could simply paste ($conjugate simp) onto the expression.
107 ;; Not done: conjugate-functions for all the inverse trigonometric functions.
109 ;; Trig like and hypergeometric like functions
111 (setf (get '%log
'conjugate-function
) 'conjugate-log
)
112 (setf (get '%plog
'conjugate-function
) 'conjugate-plog
)
113 (setf (get 'mexpt
'conjugate-function
) 'conjugate-mexpt
)
114 (setf (get '%asin
'conjugate-function
) 'conjugate-asin
)
115 (setf (get '%acos
'conjugate-function
) 'conjugate-acos
)
116 (setf (get '%atan
'conjugate-function
) 'conjugate-atan
)
117 (setf (get '%atanh
'conjugate-function
) 'conjugate-atanh
)
118 (setf (get '%asec
'conjugate-function
) 'conjugate-asec
)
119 (setf (get '%acsc
'conjugate-function
) 'conjugate-acsc
)
121 (setf (get '%bessel_j
'conjugate-function
) 'conjugate-bessel-j
)
122 (setf (get '%bessel_y
'conjugate-function
) 'conjugate-bessel-y
)
123 (setf (get '%bessel_i
'conjugate-function
) 'conjugate-bessel-i
)
124 (setf (get '%bessel_k
'conjugate-function
) 'conjugate-bessel-k
)
126 (setf (get '%hankel_1
'conjugate-function
) 'conjugate-hankel-1
)
127 (setf (get '%hankel_2
'conjugate-function
) 'conjugate-hankel-2
)
128 (setf (get '%log_gamma
'conjugate-function
) 'conjugate-log-gamma
)
130 ;; conjugate of polylogarithm li & psi
131 (setf (get '$li
'conjugate-function
) 'conjugate-li
)
132 (setf (get '$psi
'conjugate-function
) 'conjugate-psi
)
135 (setf (get '%sum
'conjugate-function
) 'conjugate-sum
)
136 (setf (get '%product
'conjugate-function
) 'conjugate-product
)
138 ;; Return true iff Maxima can prove that z is not on the
139 ;; negative real axis.
141 (defun off-negative-real-axisp (z)
142 (setq z
(trisplit z
)) ; split into real and imaginary
143 (or (eql t
(mnqp (cdr z
) 0)) ; y # 0
144 (eql t
(mgqp (car z
) 0)))) ; x >= 0
146 (defun on-negative-real-axisp (z)
147 (setq z
(trisplit z
))
148 (and (eql t
(meqp (cdr z
) 0))
149 (eql t
(mgrp 0 (car z
)))))
151 (defun off-negative-one-to-onep (z)
152 (setq z
(trisplit z
)) ; split z into real and imaginary parts
154 (eq t
(mnqp (cdr z
) 0)) ; y # 0
155 (eq t
(mgrp (car z
) 1)) ; x > 1
156 (eq t
(mgrp -
1 (car z
))))) ; -1 > x
158 (defun in-domain-of-asin (z)
159 (setq z
(trisplit z
)) ; split z into real and imaginary parts
160 (let ((x (car z
)) (y (cdr z
))) ;z = x+%i*y
162 (eq t
(mnqp y
0)) ; y # 0
164 (eq t
(mgrp x -
1)) ; x > -1
165 (eq t
(mgrp 1 x
)))))) ; x < 1
167 ;; Return conjugate(log(x)). Actually, x is a lisp list (x).
169 (defun conjugate-log (x)
171 (cond ((off-negative-real-axisp x
)
172 (take '(%log
) (take '($conjugate
) x
)))
173 ((on-negative-real-axisp x
)
174 (add (take '(%log
) (neg x
)) (mul -
1 '$%i
'$%pi
)))
175 (t (list '($conjugate simp
) (take '(%log
) x
)))))
178 ;; Return conjugate(plog(x)); again, x is the CL list (x).
179 (defun conjugate-plog (x)
181 (cond ((off-negative-real-axisp x
)
182 (take '(%plog
) (take '($conjugate
) x
)))
183 ((on-negative-real-axisp x
)
184 (add (take '(%plog
) (neg x
)) (mul -
1 '$%i
'$%pi
)))
185 (t (list '($conjugate simp
) (take '(%plog
) x
)))))
187 ;; Return conjugate(x^p), where e = (x, p). Suppose x isn't on the negative real axis.
188 ;; Then conjugate(x^p) == conjugate(exp(p * log(x))) == exp(conjugate(p) * conjugate(log(x)))
189 ;; == exp(conjugate(p) * log(conjugate(x)) = conjugate(x)^conjugate(p). Thus, when
190 ;; x is off the negative real axis, commute the conjugate with ^. Also if p is an integer
191 ;; ^ commutes with the conjugate.
193 ;; We don't need to call $ratdisrep before checking if p is a declared integer--the
194 ;; simpcheck at the top level of simp-conjugate does that for us. So we can call
195 ;; maxima-integerp on p instead of using $featurep.
197 ;; The rule that is commented out is, I think, correct, but I'm not sure how useful it is and
198 ;; the testsuite plus the share testsuite never use this rule. For now, let's keep
201 ;; Running the testsuite plus the share testsuite calls conjugate-mexpt 63,441
202 ;; times. This is far more times than all the other conjugate functions. Of these
203 ;; calls, the exponent is an integer 63,374 times. So for efficiency, we check
204 ;; (maxima-integerp p) first.
206 ;; The case of a nounform return only happens 9 times. For the nounform return, the power has
207 ;; been simplified at the higher level. So at least for running the testsuite, we shouldn't
208 ;; worry all that much about re-simplifying the power for the nounform return.
210 (defun conjugate-mexpt (e)
211 (let ((x (first e
)) (p (second e
)))
212 (cond ((or (maxima-integerp p
) (off-negative-real-axisp x
))
213 (power (take '($conjugate
) x
) (take '($conjugate
) p
)))
214 ;((on-negative-real-axisp x) ;conjugate(x^p) = exp(-%i %pi conjugate(p)) (-x)^p
215 ; (setq p (take '($conjugate) p))
216 ; (mul (power '$%e (mul -1 '$%i '$%pi p)) (power (mul -1 x) p)))
218 (list '($conjugate simp
) (power x p
))))))
220 (defun conjugate-sum (e)
221 (if (and ($featurep
(third e
) '$real
) ($featurep
(fourth e
) '$real
))
222 (take '(%sum
) (take '($conjugate
) (first e
)) (second e
) (third e
) (fourth e
))
223 (list '($conjugate simp
) (simplifya (cons '(%sum
) e
) t
))))
225 (defun conjugate-product (e)
226 (if (and ($featurep
(third e
) '$real
) ($featurep
(fourth e
) '$real
))
227 (take '(%product
) (take '($conjugate
) (first e
)) (second e
) (third e
) (fourth e
))
228 (list '($conjugate simp
) (simplifya (cons '(%product
) e
) t
))))
230 (defun conjugate-asin (x)
232 (if (in-domain-of-asin x
) (take '(%asin
) (take '($conjugate
) x
))
233 (list '($conjugate simp
) (take '(%asin
) x
))))
235 (defun conjugate-acos (x)
237 (if (in-domain-of-asin x
) (take '(%acos
) (take '($conjugate
) x
))
238 (list '($conjugate simp
) (take '(%acos
) x
))))
240 (defun conjugate-acsc (x)
242 (if (off-negative-one-to-onep x
) (take '(%acsc
) (take '($conjugate
) x
))
243 (list '($conjugate simp
) (take '(%acsc
) x
))))
245 (defun conjugate-asec (x)
247 (if (off-negative-one-to-onep x
) (take '(%asec
) (take '($conjugate
) x
))
248 (list '($conjugate simp
) (take '(%asec
) x
))))
250 (defun conjugate-atan (x)
253 (setq xx
(mul '$%i x
))
254 (if (in-domain-of-asin xx
)
255 (take '(%atan
) (take '($conjugate
) x
))
256 (list '($conjugate simp
) (take '(%atan
) x
)))))
258 ;; atanh and asin are entire on the same set; DLMF http://dlmf.nist.gov/4.37.F1 and
259 ;; http://dlmf.nist.gov/4.23.F1
261 (defun conjugate-atanh (x)
263 (if (in-domain-of-asin x
) (take '(%atanh
) (take '($conjugate
) x
))
264 (list '($conjugate simp
) (take '(%atanh
) x
))))
266 ;; Integer order Bessel functions are entire; thus they commute with the
267 ;; conjugate (Schwartz refection principle). But non-integer order Bessel
268 ;; functions are not analytic along the negative real axis. Notice that DLMF
269 ;; http://dlmf.nist.gov/10.11.E9 isn't correct; we have, for example
270 ;; conjugate(bessel_j(1/2,-1)) =/= bessel_j(1/2,conjugate(-1))
272 (defun conjugate-bessel-j (z)
273 (let ((n (first z
)) (x (second z
)))
274 (if (or ($featurep n
'$integer
) (off-negative-real-axisp x
))
275 (take '(%bessel_j
) (take '($conjugate
) n
) (take '($conjugate
) x
))
276 (list '($conjugate simp
) (simplifya (cons '(%bessel_j
) z
) t
)))))
278 (defun conjugate-bessel-y (z)
279 (let ((n (first z
)) (x (second z
)))
280 (if (off-negative-real-axisp x
)
281 (take '(%bessel_y
) (take '($conjugate
) n
) (take '($conjugate
) x
))
282 (list '($conjugate simp
) (simplifya (cons '(%bessel_y
) z
) t
)))))
284 (defun conjugate-bessel-i (z)
285 (let ((n (first z
)) (x (second z
)))
286 (if (or ($featurep n
'$integer
) (off-negative-real-axisp x
))
287 (take '(%bessel_i
) (take '($conjugate
) n
) (take '($conjugate
) x
))
288 (list '($conjugate simp
) (simplifya (cons '(%bessel_i
) z
) t
)))))
290 (defun conjugate-bessel-k (z)
291 (let ((n (first z
)) (x (second z
)))
292 (if (off-negative-real-axisp x
)
293 (take '(%bessel_k
) (take '($conjugate
) n
) (take '($conjugate
) x
))
294 (list '($conjugate simp
) (simplifya (cons '(%bessel_k
) z
) t
)))))
296 (defun conjugate-hankel-1 (z)
297 (let ((n (first z
)) (x (second z
)))
298 (if (off-negative-real-axisp x
)
299 (take '(%hankel_2
) (take '($conjugate
) n
) (take '($conjugate
) x
))
300 (list '($conjugate simp
) (simplifya (cons '(%hankel_1
) z
) t
)))))
302 (defun conjugate-hankel-2 (z)
303 (let ((n (first z
)) (x (second z
)))
304 (if (off-negative-real-axisp x
)
305 (take '(%hankel_1
) (take '($conjugate
) n
) (take '($conjugate
) x
))
306 (list '($conjugate simp
) (simplifya (cons '(%hankel_2
) z
) t
)))))
308 (defun conjugate-log-gamma (z)
310 (if (off-negative-real-axisp z
)
311 (take '(%log_gamma
) (take '($conjugate
) z
))
312 (list '($conjugate simp
) (take '(%log_gamma
) z
))))
314 ;; conjugate of polylogarithm li[s](x), where z = (s,x). We have li[s](x) = x+x^2/2^s+x^3/3^s+...
315 ;; Since for all integers k, we have conjugate(x^k/k^s) = conjugate(x)^k/k^conjugate(s), we
316 ;; commute conjugate with li.
317 (defun conjugate-li (z)
318 (let ((s (take '($conjugate
) (first z
))) (x (take '($conjugate
) (second z
))))
319 (take '(mqapply) `(($li array
) ,s
) x
)))
321 (defun conjugate-psi (z)
322 (let ((s (take '($conjugate
) (first z
))) (x (take '($conjugate
) (second z
))))
323 (take '(mqapply) `(($psi array
) ,s
) x
)))
325 ;; When all derivative variables & orders are real, commute the derivative with
327 (defun conjugate-derivative (z)
328 (cond ((every #'manifestly-real-p
(cdr z
))
329 (setq z
(cons (take '($conjugate
) (car z
)) (cdr z
)))
330 (simplifya (cons (list '%derivative
) z
) t
))
332 (list '($conjugate simp
) (simplifya (cons (list '%derivative
) z
) t
)))))
334 (setf (get '%derivative
'conjugate-function
) 'conjugate-derivative
)
336 ;; When a function maps "everything" into the reals, put real-valued on the
337 ;; property list of the function name. This duplicates some knowledge that
338 ;; $rectform has. So it goes.
340 (setf (get '%imagpart
'real-valued
) t
)
341 (setf (get 'mabs
'real-valued
) t
)
342 (setf (get '%realpart
'real-valued
) t
)
343 (setf (get '%carg
'real-valued
) t
)
344 (setf (get '$ceiling
'real-valued
) t
)
345 (setf (get '$floor
'real-valued
) t
)
346 (setf (get '$mod
'real-valued
) t
)
347 (setf (get '$unit_step
'real-valued
) t
)
348 (setf (get '$charfun
'real-valued
) t
)
351 ;; The function manifestly-real-p makes some effort to determine if its input is
354 ;; manifestly-real-p isn't a great name, but it's OK. Since (manifestly-real-p '$inf) --> true
355 ;; it might be called manifestly-extended-real-p. A nonscalar isn't real.
357 ;; There might be some advantage to requiring that the subscripts to a $subvarp
358 ;; all be real. Why? Well li[n] maps reals to reals when n is real, but li[n] does
359 ;; not map the reals to reals when n is nonreal.
361 (defun manifestly-real-p (e)
366 (not (manifestly-pure-imaginary-p e
))
367 (not (manifestly-complex-p e
))
368 (not (manifestly-nonreal-p e
)))
369 (and (consp e
) (consp (car e
)) (get (caar e
) 'real-valued
)) ;F(xxx), where F is declared real-valued
370 (and ($subvarp e
) (manifestly-real-p ($op e
)))))) ;F[n], where F is declared real-valued
372 ;; The function manifestly-pure-imaginary-p makes some effort to determine if its input is
375 (defun manifestly-pure-imaginary-p (e)
381 (and (symbolp e
) (kindp e
'$imaginary
) (not ($nonscalarp e
)))
382 (and ($subvarp e
) (manifestly-pure-imaginary-p ($op e
)))))
383 ;; For now, let's use $csign on constant expressions only; once $csign improves,
384 ;; the ban on nonconstant expressions can be removed.
385 (and ($constantp e
) (not (eq '$und e
)) (not (eq '$ind e
)) (eq '$imaginary
($csign e
))))))
387 ;; Don't use (kindp e '$complex)!
389 (defun manifestly-complex-p (e)
391 (or (and (symbolp e
) (decl-complexp e
) (not ($nonscalarp e
)))
393 (and ($subvarp e
) (manifestly-complex-p ($op e
)) (not ($nonscalarp e
))))))
395 (defun manifestly-nonreal-p (e)
396 (and (symbolp e
) (or (member e
`($und $ind t nil
)) ($nonscalarp e
))))
398 ;; We could make commutes_with_conjugate and maps_to_reals features. But I
399 ;; doubt it would get much use.
401 (defun simp-conjugate (e f z
)
403 (setq e
(simpcheck (cadr e
) z
)) ; simp and disrep if necessary
405 (cond ((complexp e
) (conjugate e
)) ; never happens, but might someday.
406 ((manifestly-real-p e
) e
)
407 ((manifestly-pure-imaginary-p e
) (mul -
1 e
))
408 ((or (manifestly-nonreal-p e
) ($mapatom e
))
409 (list '($conjugate simp
) e
))
411 ((op-equalp e
'$conjugate
) (car (margs e
)))
413 ((and (symbolp (mop e
)) (get (mop e
) 'real-valued
)) e
)
415 ((and (symbolp (mop e
)) (get (mop e
) 'commutes-with-conjugate
))
416 (simplify (cons (list (mop e
)) (mapcar #'(lambda (s) (take '($conjugate
) s
)) (margs e
)))))
418 ((setq f
(and (symbolp (mop e
)) (get (mop e
) 'conjugate-function
)))
419 (funcall f
(margs e
)))
421 ;;subscripted functions
422 ((setq f
(and ($subvarp
(mop e
)) (get (caar (mop e
)) 'conjugate-function
)))
423 (funcall f
(append (margs (mop e
)) (margs e
))))
426 (list '($conjugate simp
) e
))))
