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))
);