complex conjugate and modulus

[John Lapeyre]

Hi,

I am trying to write a function to rewrite
z*conjugate(z) as abs(z)^2 , and do this when z is
other expressions, as well.

Based on my attempts and reading the archives, I think this is
impossible.

What I have succeeds for the following,

declare([y,z], complex);
conj_to_abs( 2 * y * conjugate(y) * z * conjugate(z) * w * conjugate(w) );

which gives :

2*w^2*abs(y)^2*abs(z)^2

-------
But fails for

conj_to_abs(cos(z) * conjugate(cos(z)));

which gives
cos(z)*cos(conjugate(z))

This last result is mathematically correct, but I dont know how to
stop the replacement conjugate(cos(z)) --> cos(conjugate(z)), which
happens before it is passed to my routine.

-------
It also fails for

conj_to_abs(sqrt(z)*conjugate(sqrt(z)));

which gives

conjugate(sqrt(z))*sqrt(z)

In this case, my pattern matching routine refuses to recognize
sqrt(z) = sqrt(z) as true.
I put in several print statements right before the comparison, and
it still fails. At the command line, I can do

is ( part( conjugate(sqrt(z)), 1) = sqrt(z));

with the result 'true'. But it fails in my routine. It must have to
do with some internal representation.

If I set simp:false globally, then this example *does* give the correct
result. But if I set simp:true, then
abs(sqrt(z))^2 is simplified to z , which is incorrect.

Still I don't understand why the matching routine fails above as I
wrote. Maybe I could define another function for abs that has the properties
that I want.

--------

This one does give the desired result

conj_to_abs( foo(z) * conjugate(foo(z)));

with

abs(foo(z))^2,

as long as I globally set simp:false,

otherwise, after the value is returned, it is simplified (incorrectly)
to
foo(z)^2

I know there has been quite a bit of talk about how foo(z) and cos(z)
are treated, with no easy solution at hand.

---------------

I have a couple of questions.

1) What kind of facilities are there for dumping various
representations of an expression ? There must be something
like this, because it would be so useful, but I don't know
of such a thing.

2) I guess there is not much opportunity for fine grained simplification.
In order to get my routine to work, I have to turn of simplification in
blocks. Then I have to simplifiy the loop variable "i" explicitly to
avoid getting part(e,1+1). (This took a long while to figure out, in part
because the error messages don't give line numbers in the input file.) Or
is there something better than simp:false ?

Following is the code I am using, which I am sure is not very idiomatic:

Thanks,

John

-------------------
/* test if an expression and its conjuate appear in a list of args
to "*" */
conj_multp(e) := block( [p1, p2, clis, noclist, cmatch : false, simp:false],
if (not atom(e)) and (op(e) = "*") then (
clis : sublist(args(e),lambda([e], (not atom(e)) and op(e) = conjugate)),
noclis : sublist(args(e),lambda([e], atom(e) or not (op(e) =
conjugate))),
for p1 in clis do (
for p2 in noclis do (
print ( " multp full ", p1," -- ", p2),
print ( " multp parts >> ", part(p1,1), " -- ", p2),
/*
print ( " multp ", op(part(p1,1)), args(part(p1,1)), op(p2),
args(p2)),
print ( " multp inpart ", op(inpart(p1,1)), args(inpart(p1,1)),
op(p2), args(p2)),
*/
if part(p1,1) = p2 then (
/* print( " got match"), */
cmatch : true,
return(true)
)
)
),
return( cmatch)
)
);



block ([simp:false],
local(ee),
matchdeclare(ee, conj_multp),
defrule(conj_mult_rule, ee,
block ( [i,p1, clis, noclis, newclis : [], mflag ],
print("entering rule >> ", ee),
clis : sublist(args(ee),lambda([e], (not atom(e)) and op(e) =
conjugate)),
noclis : sublist(args(ee),lambda([e], atom(e) or not (op(e) =
conjugate))),
print(clis,noclis),
for p1 in clis do (
mflag : false,
for i thru length(noclis) do (
i : ratsimp(i),
print(" in rule ", p1),
print(" in rule part ", part(p1,1)),
print(" in rule noclis full i=",i, noclis),
print(" in rule noclis ", part(noclis,i)),
if part(p1,1) = part(noclis,i) then (
noclis : substpart(abs(part(noclis,i))^2,noclis,i),
mflag : true
)
),
if mflag = false then newclis : cons(p1,newclis)
),
print( newclis, " -- " , noclis),
newclis : append(newclis,noclis),
if ( length(newclis) > 1) then return( apply("*", newclis)),
return (part(newclis,1))
)
),
conj_to_abs(e) := block([simp:false] ,apply1( e, conj_mult_rule))
);

[Barton Willis]

Two suggestions:

(1) Instead of setting simp to false, try subst('fake_conjugate,
'conjugate, ...). After processing, revert the substitution. It might
help to declare fake_conjugate to be additive and multiplicative (but
it might hurt too).

(2) Instead of apply1 + a rule, try using ratsubst.

Here is some CL code and some tests--expect bugs!

;; start of file conjugate_to_abs.lisp ------------

(mfuncall '$declare 'non-simp-conjugate '$additive)
(mfuncall '$declare 'non-simp-conjugate '$multiplicative)

(defmspec $conj_to_abs (e)
(setq e ($substitute 'non-simp-conjugate '$conjugate (cadr e)))
($substitute '$conjugate 'non-simp-conjugate (conj-to-abs e e)))

(defun conj-to-abs (e p)
(cond ((or ($mapatom e) ($mapatom p)) e)

((and (not ($mapatom p)) (eq (mop p) 'non-simp-conjugate))
(setq p (cadr p))
($ratsubst (power (take '($cabs) p) 2) (mul p (take '(non-simp-conjugate)
p)) e))

(t
(mapcar #'(lambda (s) (setq e (conj-to-abs e s))) (cdr p))
e)))

;; end of file conjugate_to_abs.lisp ------------

Examples:

(%i12) load("conjugate_to_abs.lisp")$
(%i13) declare([x,y,z],complex)$

(%i14) conj_to_abs(x * conjugate(5*x));
(%o14) 5*cabs(x)^2

(%i15) conj_to_abs(x^3 * conjugate(5*x));
(%o15) 5*x^2*cabs(x)^2

(%i16) conj_to_abs( 2 * y * conjugate(y) * z * conjugate(z) * w * conjugate
(w));
(%o16) 2*cabs(w)^2*cabs(y)^2*cabs(z)^2

(%i17) conj_to_abs(sqrt(z)*conjugate(sqrt(z)));
(%o17) cabs(sqrt(z))^2

(%i18) conj_to_abs( foo(z) * conjugate(foo(z)));
(%o18) cabs(foo(z))^2

(%i19) conj_to_abs(z * conjugate(42 + z));
(%o19) cabs(z)^2+42*z

(%i20) conj_to_abs((z + 42) * conjugate(z));
(%o20) 42*conjugate(z)+cabs(z)^2

(%i21) conj_to_abs(z^5 * conjugate(z + 12));
(%o21) z^4*cabs(z)^2+12*z^5

(%i22) conj_to_abs(z^5 * conjugate(z^2+12));
(%o22) z^3*cabs(z^2)^2+12*z^5

Here, it would be better if non-simp-conjugate was not additive :(

(%i23) conj_to_abs((x + y) * conjugate(x + y));
(%o23) x*conjugate(y)+cabs(y)^2+conjugate(x)*y+cabs(x)^2

Again, my code is almost surely buggy, but maybe it is a start for
something
that will work for you. One problem is

(%i11) conj_to_abs(6.7b0 + conjugate(x));
`rat' replaced 6.7B0 by 67/10 = 6.7B0
(%o11) (10*conjugate(x)+67)/10

I don't know of a way around that---for doubles, keepfloat takes care of
this,
but not for big floats.

Barton

-----maxima-***@math.utexas.edu wrote: -----

>I am trying to write a function to rewrite
> z*conjugate(z) as abs(z)^2 , and do this when z is
>other expressions, as well.

[John Lapeyre]

For some reason, this lispcode does not work on my system. I
was using maxima 5.13 with sbcl, the lisp below failed to
load.

I just upgraded to the latest debian package.
maxima 5.15.0, GCL 2.6.7
Now the code loads silently, but calls to conj_to_abs fail

(%i1) load("conj_to_abs.lisp");
(\%o1) \verb|conj_to_abs.lisp|
(%i2) declare([x,y,z],complex)$

(%i3) conj_to_abs(x * conjugate(5*x));

Maxima encountered a Lisp error:

Error in LISP:LAMBDA-CLOSURE [or a callee]: The variable   is unbound.

---------------
> ;; start of file conjugate_to_abs.lisp ------------
>
> (mfuncall '$declare 'non-simp-conjugate '$additive)
> (mfuncall '$declare 'non-simp-conjugate '$multiplicative)
>
> (defmspec $conj_to_abs (e)
> (setq e ($substitute 'non-simp-conjugate '$conjugate (cadr e)))
> ($substitute '$conjugate 'non-simp-conjugate (conj-to-abs e e)))
>
> (defun conj-to-abs (e p)
> (cond ((or ($mapatom e) ($mapatom p)) e)
>
> ((and (not ($mapatom p)) (eq (mop p) 'non-simp-conjugate))
> (setq p (cadr p))
> ($ratsubst (power (take '($cabs) p) 2) (mul p (take '(non-simp-conjugate)
> p)) e))
>
> (t
> (mapcar #'(lambda (s) (setq e (conj-to-abs e s))) (cdr p))
> e)))
>
> ;; end of file conjugate_to_abs.lisp ------------

[Barton Willis]

-----maxima-***@math.utexas.edu wrote: -----

>For some reason, this lispcode does not work on my system. I
>was using maxima 5.13 with sbcl, the lisp below failed to
>load.

With CVS maxima + SBCL 1.0.13 my code loads and my examples OK. Maybe
fixing
it isn't worthwhile; the Maxima code from Stavros is likely better; with my
code things like

(%i27) e : sin(x) * conjugate(sin(x));
(%o27) sin(x)*sin(conjugate(x))

Doesn't work:

(%i28) conj_to_abs(''e);
(%o28) sin(x)*sin(conjugate(x))

Works:

(%i29) conj_to_abs(sin(x) * conjugate(sin(x)));
(%o29) cabs(sin(x))^2

Barton

[Barton Willis]

P.S. Yet another version:

(%i33) declare(z,complex);
(%o33) done

(%i34) conj_to_abs (e) := block([args],
if mapatom(e) then e
elseif ?mtimesp(e) then (
args : args(e),
for p in args do e : ratsubst('cabs(p)^2, p * conjugate(p),e),
e)
else apply(op(e), map('conj_to_abs, args(e))));

(%i36) sin(z) * conjugate(sin(z));
(%o36) sin(z)*sin(conjugate(z))

(%i37) conj_to_abs(%);
(%o37) cabs(sin(z))^2

(%i38) ev(%,cabs);
(%o38) cosh(imagpart(z))^2*sin(realpart(z))^2+sinh(imagpart(z))^2*cos
(realpart(z))^2

Barton

[Stavros Macrakis]

On Sat, Aug 9, 2008 at 12:10 AM, John Lapeyre <***@johnlapeyre.com> wrote:
> I am trying to write a function to rewrite z*conjugate(z) as abs(z)^2 , and do this when z is other expressions, as well.
>
> Based on my attempts and reading the archives, I think this is impossible.

That's a bold claim!!

What I would suggest is to *not* try to evade or override built-in
Maxima simplification, but to work with it. If Maxima simplifies
sin(x)*conjugate(sin(x)) to sin(x)*sin(conjugate(x)), then the latter
form is what you need to simplify.

Here is one approach. Note that conjsimp_pair makes an assumption
which you may want to revisit. You may or may not want to replace
standalone conjugate(w) with abs(2)^2/w, but that's easy enough to
add.

/* conjsimp replaces z^a*conjugate(z^a) by abs(z)^(2*a) in an expression */
/* Copyright 2008 Stavros Macrakis; licensed under LGPL */

conjugate_to_abs(ex):=
block([inflag:true,
args,
prod: 1], /* return value */
if mapatom(ex) then ex
elseif part(ex,0)="*"
then (args: args(ex),
while (n:find_pair(conjsimp_pair,args))#false do
(args: delete(n[1],delete(n[2],args)),
prod: prod*n[3]),
prod*xreduce("*",args))
else map(conjugate_to_abs,ex))$

/* find first i<j s.t. f(l[i],l[j]) is not false
and return [ l[i], l[j], f(l[i],l[j]) ]
if no such i<j, then return false */
find_pair(f,l):=
block([lenl:length(l),val],
catch(
(for i thru lenl-1 do
for j:i+1 thru lenl do
if false # val:f(l[i],l[j])
then throw([l[i],l[j],val]),
false)))$

/* assumes that z2 is the one that contains the conjugate, which it
usually is */
conjsimp_pair(be1,be2):=
(be1: base_exp(be1),
be2: base_exp(be2),
if be1[1]=conjugate(be2[1])
then abs(be1[1])^(2*be2[2]) * be1[1]^(be1[2]-be2[2])
else false)$

base_exp(ex):=
block([inflag:true],
if mapatom(ex) or part(ex,0) # "^"
then [ex,1]
else args(ex))$

>>>>>>>>>>>>examples<<<<<<<<<<<<<<<
z*conjugate(z) => abs(z)^2
w*z*conjugate(z) => w*abs(z)^2
conjugate(w)*z*conjugate(z) => conjugate(w)*abs(z)^2
sin(z)*sin(conjugate(z)) => sin(z)^2
z^3/conjugate(z) => z^4/abs(z)^2
conjugate(w)*conjugate(z)/z^2 => conjugate(w)*abs(z)^2/z^3

[Stavros Macrakis]

In my previous message, I posted some code that I hope will help with
your immediate problem. Here, some discussion on your specific
questions.

On Sat, Aug 9, 2008 at 12:10 AM, John Lapeyre <***@johnlapeyre.com> wrote:

> This last result is mathematically correct, but I dont know how to
> stop the replacement conjugate(cos(z)) --> cos(conjugate(z)), which
> happens before it is passed to my routine.

As a general rule, it is a bad idea to block Maxima's simplifications,
which standardize the form of expressions so that further
simplifications can be applied. For example, if conjugate(cos(x))
weren't simplified to cos(conjugate(x)), Maxima wouldn't be able to
recognize that conjugate(cos(x)) - cos(conjugate(x)) = 0. Or that
sqrt(x)=x^(1/2).

That said, I understand that it is frustrating when these automatic
simplifications put things in a form that you don't want or that is
harder to work with.

> In this case, my pattern matching routine refuses to recognize
> sqrt(z) = sqrt(z) as true.

Without trying to debug in detail, I would guess that this is because
you've turned simplification off. The standard internal
representation of sqrt(x) is x^(1/2); the transformation is done by
simplification.

> But if I set simp:true, then abs(sqrt(z))^2 is simplified to z , which is incorrect.

This is a bug. Sorry. I believe Dieter Kaiser is looking into this
and related bugs.

> 1) What kind of facilities are there for dumping various representations of an expression?

The visible form of the expression is what "part" works off of. So
part(a/b) => "/". The internal form is what "inpart" works off of, so
inpart(a/b,0) => "*" and inargs(a/b) => [a, 1/b] where
inargs(ex):=block([inflag:true],args(ex)). You can also look at the
internal (Lisp-based) representation with ?print(expr). As a general
rule, if you are manipulating expressions programmatically, I would
strongly recommend using inpart and inargs; part is extremely
inefficient because it reformats the *whole* expression every time it
looks at it....

> 2) I guess there is not much opportunity for fine grained simplification.
> In order to get my routine to work, I have to turn of simplification in
> blocks. ...Or is there something better than simp:false ?

simp:false is almost always a bad thing to do, because everything in
Maxima depends on simplification, even arithmetic as you can see.
Many parts of Maxima do not work correctly with simp:false.

The problem with abs(sqrt(z))^2 is not something you will be able to
work around by turning simplification off.

Hope this is helpful.

-s