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