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.
29 #-gcl
(:load-toplevel
:execute
)
30 (let (($context
'$global
) (context '$global
))
31 (meval '(($declare
) $conjugate $complex
))
32 ;; Let's remove built-in symbols from list for user-defined properties.
33 (setq $props
(remove '$conjugate $props
))))
35 ;; When a function commutes with the conjugate, give the function the
36 ;; commutes-with-conjugate property. The log function commutes with
37 ;; the conjugate on all of C except on the negative real axis. Thus
38 ;; log does not get the commutes-with-conjugate property. Instead,
39 ;; log gets the conjugate-function property.
41 ;; What important functions have I missed?
43 ;; (1) Arithmetic operators
45 (setf (get 'mplus
'commutes-with-conjugate
) t
)
46 (setf (get 'mtimes
'commutes-with-conjugate
) t
)
47 ;(setf (get 'mnctimes 'commutes-with-conjugate) t) ;; generally I think users will want this
48 (setf (get '%signum
'commutes-with-conjugate
) t
) ;; x=/=0, conjugate(signum(x)) = conjugate(x/abs(x)) = signum(conjugate(x))
49 ;; Trig-like functions and other such functions
51 (setf (get '%cosh
'commutes-with-conjugate
) t
)
52 (setf (get '%sinh
'commutes-with-conjugate
) t
)
53 (setf (get '%tanh
'commutes-with-conjugate
) t
)
54 (setf (get '%sech
'commutes-with-conjugate
) t
)
55 (setf (get '%csch
'commutes-with-conjugate
) t
)
56 (setf (get '%coth
'commutes-with-conjugate
) t
)
57 (setf (get '%cos
'commutes-with-conjugate
) t
)
58 (setf (get '%sin
'commutes-with-conjugate
) t
)
59 (setf (get '%tan
'commutes-with-conjugate
) t
)
60 (setf (get '%sec
'commutes-with-conjugate
) t
)
61 (setf (get '%csc
'commutes-with-conjugate
) t
)
62 (setf (get '%cot
'commutes-with-conjugate
) t
)
63 (setf (get '$atan2
'commutes-with-conjugate
) t
)
65 (setf (get '%jacobi_cn
'commutes-with-conjugate
) t
)
66 (setf (get '%jacobi_sn
'commutes-with-conjugate
) t
)
67 (setf (get '%jacobi_dn
'commutes-with-conjugate
) t
)
69 (setf (get '%gamma
'commutes-with-conjugate
) t
)
70 (setf (get '$pochhammer
'commutes-with-conjugate
) t
)
74 (setf (get '$matrix
'commutes-with-conjugate
) t
)
75 (setf (get 'mlist
'commutes-with-conjugate
) t
)
76 (setf (get '$set
'commutes-with-conjugate
) t
)
80 (setf (get 'mequal
'commutes-with-conjugate
) t
)
81 (setf (get 'mnotequal
'commutes-with-conjugate
) t
)
82 (setf (get '%transpose
'commutes-with-conjugate
) t
)
86 (setf (get '$max
'commutes-with-conjugate
) t
)
87 (setf (get '$min
'commutes-with-conjugate
) t
)
89 ;; When a function has the conjugate-function property, use a non-generic function to conjugate it.
90 ;; The argument to a conjugate function for an operator op is the CL list of arguments to op. For
91 ;; example, the conjugate function for log gets the argument for log, not the expression log(x).
92 ;; It would be a bit more efficient if a conjugate function received the full expression--that
93 ;; way for a pure nounform return (for example, return conjugate(log(x))), a conjugate function
94 ;; would not not need to apply the operator to the argument to the conjugate function, instead it
95 ;; could simply paste ($conjugate simp) onto the expression.
97 ;; Not done: conjugate-functions for all the inverse trigonometric functions.
99 ;; Trig like and hypergeometric like functions
101 (setf (get '%log
'conjugate-function
) 'conjugate-log
)
102 (setf (get '%plog
'conjugate-function
) 'conjugate-plog
)
103 (setf (get 'mexpt
'conjugate-function
) 'conjugate-mexpt
)
104 (setf (get '%asin
'conjugate-function
) 'conjugate-asin
)
105 (setf (get '%acos
'conjugate-function
) 'conjugate-acos
)
106 (setf (get '%atan
'conjugate-function
) 'conjugate-atan
)
107 (setf (get '%atanh
'conjugate-function
) 'conjugate-atanh
)
108 (setf (get '%asec
'conjugate-function
) 'conjugate-asec
)
109 (setf (get '%acsc
'conjugate-function
) 'conjugate-acsc
)
111 (setf (get '%bessel_j
'conjugate-function
) 'conjugate-bessel-j
)
112 (setf (get '%bessel_y
'conjugate-function
) 'conjugate-bessel-y
)
113 (setf (get '%bessel_i
'conjugate-function
) 'conjugate-bessel-i
)
114 (setf (get '%bessel_k
'conjugate-function
) 'conjugate-bessel-k
)
116 (setf (get '%hankel_1
'conjugate-function
) 'conjugate-hankel-1
)
117 (setf (get '%hankel_2
'conjugate-function
) 'conjugate-hankel-2
)
118 (setf (get '%log_gamma
'conjugate-function
) 'conjugate-log-gamma
)
120 ;; conjugate of polylogarithm li & psi
121 (setf (get '$li
'conjugate-function
) 'conjugate-li
)
122 (setf (get '$psi
'conjugate-function
) 'conjugate-psi
)
125 (setf (get '%sum
'conjugate-function
) 'conjugate-sum
)
126 (setf (get '%product
'conjugate-function
) 'conjugate-product
)
128 ;; Return true iff Maxima can prove that z is not on the
129 ;; negative real axis.
131 (defun off-negative-real-axisp (z)
132 (setq z
(trisplit z
)) ; split into real and imaginary
133 (or (eql t
(mnqp (cdr z
) 0)) ; y # 0
134 (eql t
(mgqp (car z
) 0)))) ; x >= 0
136 (defun on-negative-real-axisp (z)
137 (setq z
(trisplit z
))
138 (and (eql t
(meqp (cdr z
) 0))
139 (eql t
(mgrp 0 (car z
)))))
141 (defun off-negative-one-to-onep (z)
142 (setq z
(trisplit z
)) ; split z into real and imaginary parts
144 (eq t
(mnqp (cdr z
) 0)) ; y # 0
145 (eq t
(mgrp (car z
) 1)) ; x > 1
146 (eq t
(mgrp -
1 (car z
))))) ; -1 > x
148 (defun in-domain-of-asin (z)
149 (setq z
(trisplit z
)) ; split z into real and imaginary parts
150 (let ((x (car z
)) (y (cdr z
))) ;z = x+%i*y
152 (eq t
(mnqp y
0)) ; y # 0
154 (eq t
(mgrp x -
1)) ; x > -1
155 (eq t
(mgrp 1 x
)))))) ; x < 1
157 ;; Return conjugate(log(x)). Actually, x is a lisp list (x).
159 (defun conjugate-log (x)
161 (cond ((off-negative-real-axisp x
)
162 (take '(%log
) (take '($conjugate
) x
)))
163 ((on-negative-real-axisp x
)
164 (add (take '(%log
) (neg x
)) (mul -
1 '$%i
'$%pi
)))
165 (t (list '($conjugate simp
) (take '(%log
) x
)))))
168 ;; Return conjugate(plog(x)); again, x is the CL list (x).
169 (defun conjugate-plog (x)
171 (cond ((off-negative-real-axisp x
)
172 (take '(%plog
) (take '($conjugate
) x
)))
173 ((on-negative-real-axisp x
)
174 (add (take '(%plog
) (neg x
)) (mul -
1 '$%i
'$%pi
)))
175 (t (list '($conjugate simp
) (take '(%plog
) x
)))))
177 ;; Return conjugate(x^p), where e = (x, p). Suppose x isn't on the negative real axis.
178 ;; Then conjugate(x^p) == conjugate(exp(p * log(x))) == exp(conjugate(p) * conjugate(log(x)))
179 ;; == exp(conjugate(p) * log(conjugate(x)) = conjugate(x)^conjugate(p). Thus, when
180 ;; x is off the negative real axis, commute the conjugate with ^. Also if p is an integer
181 ;; ^ commutes with the conjugate.
183 ;; We don't need to call $ratdisrep before checking if p is a declared integer--the
184 ;; simpcheck at the top level of simp-conjugate does that for us. So we can call
185 ;; maxima-integerp on p instead of using $featurep.
187 ;; The rule that is commented out is, I think, correct, but I'm not sure how useful it is and
188 ;; the testsuite plus the share testsuite never use this rule. For now, let's keep
191 ;; Running the testsuite plus the share testsuite calls conjugate-mexpt 63,441
192 ;; times. This is far more times than all the other conjugate functions. Of these
193 ;; calls, the exponent is an integer 63,374 times. So for efficiency, we check
194 ;; (maxima-integerp p) first.
196 ;; The case of a nounform return only happens 9 times. For the nounform return, the power has
197 ;; been simplified at the higher level. So at least for running the testsuite, we shouldn't
198 ;; worry all that much about re-simplifying the power for the nounform return.
200 (defun conjugate-mexpt (e)
201 (let ((x (first e
)) (p (second e
)))
202 (cond ((or (maxima-integerp p
) (off-negative-real-axisp x
))
203 (power (take '($conjugate
) x
) (take '($conjugate
) p
)))
204 ;((on-negative-real-axisp x) ;conjugate(x^p) = exp(-%i %pi conjugate(p)) (-x)^p
205 ; (setq p (take '($conjugate) p))
206 ; (mul (power '$%e (mul -1 '$%i '$%pi p)) (power (mul -1 x) p)))
208 (list '($conjugate simp
) (power x p
))))))
210 (defun conjugate-sum (e)
211 (if (and ($featurep
(third e
) '$real
) ($featurep
(fourth e
) '$real
))
212 (take '(%sum
) (take '($conjugate
) (first e
)) (second e
) (third e
) (fourth e
))
213 (list '($conjugate simp
) (simplifya (cons '(%sum
) e
) t
))))
215 (defun conjugate-product (e)
216 (if (and ($featurep
(third e
) '$real
) ($featurep
(fourth e
) '$real
))
217 (take '(%product
) (take '($conjugate
) (first e
)) (second e
) (third e
) (fourth e
))
218 (list '($conjugate simp
) (simplifya (cons '(%product
) e
) t
))))
220 (defun conjugate-asin (x)
222 (if (in-domain-of-asin x
) (take '(%asin
) (take '($conjugate
) x
))
223 (list '($conjugate simp
) (take '(%asin
) x
))))
225 (defun conjugate-acos (x)
227 (if (in-domain-of-asin x
) (take '(%acos
) (take '($conjugate
) x
))
228 (list '($conjugate simp
) (take '(%acos
) x
))))
230 (defun conjugate-acsc (x)
232 (if (off-negative-one-to-onep x
) (take '(%acsc
) (take '($conjugate
) x
))
233 (list '($conjugate simp
) (take '(%acsc
) x
))))
235 (defun conjugate-asec (x)
237 (if (off-negative-one-to-onep x
) (take '(%asec
) (take '($conjugate
) x
))
238 (list '($conjugate simp
) (take '(%asec
) x
))))
240 (defun conjugate-atan (x)
243 (setq xx
(mul '$%i x
))
244 (if (in-domain-of-asin xx
)
245 (take '(%atan
) (take '($conjugate
) x
))
246 (list '($conjugate simp
) (take '(%atan
) x
)))))
248 ;; atanh and asin are entire on the same set; DLMF http://dlmf.nist.gov/4.37.F1 and
249 ;; http://dlmf.nist.gov/4.23.F1
251 (defun conjugate-atanh (x)
253 (if (in-domain-of-asin x
) (take '(%atanh
) (take '($conjugate
) x
))
254 (list '($conjugate simp
) (take '(%atanh
) x
))))
256 ;; Integer order Bessel functions are entire; thus they commute with the
257 ;; conjugate (Schwartz refection principle). But non-integer order Bessel
258 ;; functions are not analytic along the negative real axis. Notice that DLMF
259 ;; http://dlmf.nist.gov/10.11.E9 isn't correct; we have, for example
260 ;; conjugate(bessel_j(1/2,-1)) =/= bessel_j(1/2,conjugate(-1))
262 (defun conjugate-bessel-j (z)
263 (let ((n (first z
)) (x (second z
)))
264 (if (or ($featurep n
'$integer
) (off-negative-real-axisp x
))
265 (take '(%bessel_j
) (take '($conjugate
) n
) (take '($conjugate
) x
))
266 (list '($conjugate simp
) (simplifya (cons '(%bessel_j
) z
) t
)))))
268 (defun conjugate-bessel-y (z)
269 (let ((n (first z
)) (x (second z
)))
270 (if (off-negative-real-axisp x
)
271 (take '(%bessel_y
) (take '($conjugate
) n
) (take '($conjugate
) x
))
272 (list '($conjugate simp
) (simplifya (cons '(%bessel_y
) z
) t
)))))
274 (defun conjugate-bessel-i (z)
275 (let ((n (first z
)) (x (second z
)))
276 (if (or ($featurep n
'$integer
) (off-negative-real-axisp x
))
277 (take '(%bessel_i
) (take '($conjugate
) n
) (take '($conjugate
) x
))
278 (list '($conjugate simp
) (simplifya (cons '(%bessel_i
) z
) t
)))))
280 (defun conjugate-bessel-k (z)
281 (let ((n (first z
)) (x (second z
)))
282 (if (off-negative-real-axisp x
)
283 (take '(%bessel_k
) (take '($conjugate
) n
) (take '($conjugate
) x
))
284 (list '($conjugate simp
) (simplifya (cons '(%bessel_k
) z
) t
)))))
286 (defun conjugate-hankel-1 (z)
287 (let ((n (first z
)) (x (second z
)))
288 (if (off-negative-real-axisp x
)
289 (take '(%hankel_2
) (take '($conjugate
) n
) (take '($conjugate
) x
))
290 (list '($conjugate simp
) (simplifya (cons '(%hankel_1
) z
) t
)))))
292 (defun conjugate-hankel-2 (z)
293 (let ((n (first z
)) (x (second z
)))
294 (if (off-negative-real-axisp x
)
295 (take '(%hankel_1
) (take '($conjugate
) n
) (take '($conjugate
) x
))
296 (list '($conjugate simp
) (simplifya (cons '(%hankel_2
) z
) t
)))))
298 (defun conjugate-log-gamma (z)
300 (if (off-negative-real-axisp z
)
301 (take '(%log_gamma
) (take '($conjugate
) z
))
302 (list '($conjugate simp
) (take '(%log_gamma
) z
))))
304 ;; conjugate of polylogarithm li[s](x), where z = (s,x). We have li[s](x) = x+x^2/2^s+x^3/3^s+...
305 ;; Since for all integers k, we have conjugate(x^k/k^s) = conjugate(x)^k/k^conjugate(s), we
306 ;; commute conjugate with li.
307 (defun conjugate-li (z)
308 (let ((s (take '($conjugate
) (first z
))) (x (take '($conjugate
) (second z
))))
309 (take '(mqapply) `(($li array
) ,s
) x
)))
311 (defun conjugate-psi (z)
312 (let ((s (take '($conjugate
) (first z
))) (x (take '($conjugate
) (second z
))))
313 (take '(mqapply) `(($psi array
) ,s
) x
)))
315 ;; When all derivative variables & orders are real, commute the derivative with
317 (defun conjugate-derivative (z)
318 (cond ((every #'manifestly-real-p
(cdr z
))
319 (setq z
(cons (take '($conjugate
) (car z
)) (cdr z
)))
320 (simplifya (cons (list '%derivative
) z
) t
))
322 (list '($conjugate simp
) (simplifya (cons (list '%derivative
) z
) t
)))))
324 (setf (get '%derivative
'conjugate-function
) #'conjugate-derivative
)
326 ;; When a function maps "everything" into the reals, put real-valued on the
327 ;; property list of the function name. This duplicates some knowledge that
328 ;; $rectform has. So it goes.
330 (setf (get '%imagpart
'real-valued
) t
)
331 (setf (get 'mabs
'real-valued
) t
)
332 (setf (get '%realpart
'real-valued
) t
)
333 (setf (get '%carg
'real-valued
) t
)
334 (setf (get '$ceiling
'real-valued
) t
)
335 (setf (get '$floor
'real-valued
) t
)
336 (setf (get '$mod
'real-valued
) t
)
337 (setf (get '$unit_step
'real-valued
) t
)
338 (setf (get '$charfun
'real-valued
) t
)
341 ;; The function manifestly-real-p makes some effort to determine if its input is
344 ;; manifestly-real-p isn't a great name, but it's OK. Since (manifestly-real-p '$inf) --> true
345 ;; it might be called manifestly-extended-real-p. A nonscalar isn't real.
347 ;; There might be some advantage to requiring that the subscripts to a $subvarp
348 ;; all be real. Why? Well li[n] maps reals to reals when n is real, but li[n] does
349 ;; not map the reals to reals when n is nonreal.
351 (defun manifestly-real-p (e)
356 (not (manifestly-pure-imaginary-p e
))
357 (not (manifestly-complex-p e
))
358 (not (manifestly-nonreal-p e
)))
359 (and (consp e
) (consp (car e
)) (get (caar e
) 'real-valued
)) ;F(xxx), where F is declared real-valued
360 (and ($subvarp e
) (manifestly-real-p ($op e
)))))) ;F[n], where F is declared real-valued
362 ;; The function manifestly-pure-imaginary-p makes some effort to determine if its input is
365 (defun manifestly-pure-imaginary-p (e)
371 (and (symbolp e
) (kindp e
'$imaginary
) (not ($nonscalarp e
)))
372 (and ($subvarp e
) (manifestly-pure-imaginary-p ($op e
)))))
373 ;; For now, let's use $csign on constant expressions only; once $csign improves,
374 ;; the ban on nonconstant expressions can be removed.
375 (and ($constantp e
) (not (eq '$und e
)) (not (eq '$ind e
)) (eq '$imaginary
($csign e
))))))
377 ;; Don't use (kindp e '$complex)!
379 (defun manifestly-complex-p (e)
381 (or (and (symbolp e
) (decl-complexp e
) (not ($nonscalarp e
)))
383 (and ($subvarp e
) (manifestly-complex-p ($op e
)) (not ($nonscalarp e
))))))
385 (defun manifestly-nonreal-p (e)
386 (and (symbolp e
) (or (member e
`($und $ind t nil
)) ($nonscalarp e
))))
388 ;; We could make commutes_with_conjugate and maps_to_reals features. But I
389 ;; doubt it would get much use.
391 (defun simp-conjugate (e f z
)
393 (setq e
(simpcheck (cadr e
) z
)) ; simp and disrep if necessary
395 (cond ((complexp e
) (conjugate e
)) ; never happens, but might someday.
396 ((manifestly-real-p e
) e
)
397 ((manifestly-pure-imaginary-p e
) (mul -
1 e
))
398 ((or (manifestly-nonreal-p e
) ($mapatom e
))
399 (list '($conjugate simp
) e
))
401 ((op-equalp e
'$conjugate
) (car (margs e
)))
403 ((and (symbolp (mop e
)) (get (mop e
) 'real-valued
)) e
)
405 ((and (symbolp (mop e
)) (get (mop e
) 'commutes-with-conjugate
))
406 (simplify (cons (list (mop e
)) (mapcar #'(lambda (s) (take '($conjugate
) s
)) (margs e
)))))
408 ((setq f
(and (symbolp (mop e
)) (get (mop e
) 'conjugate-function
)))
409 (funcall f
(margs e
)))
411 ;;subscripted functions
412 ((setq f
(and ($subvarp
(mop e
)) (get (caar (mop e
)) 'conjugate-function
)))
413 (funcall f
(append (margs (mop e
)) (margs e
))))
416 (list '($conjugate simp
) e
))))
