Earlier today, I was re-reading an old post when I was inspired. The post asked how one could use typesetting rules to make the call of a procedure typeset as a binary operation. The example used in the post was something along the lines of:
a <symbol> b;
Inspired by a response to that post, I started looking at the help page for the
define command and doing some testing in Maple. I was ready to make a reply to that post, when I stumbled upon a mind-boggling error. In the time it's taken me to track down that error (and run various errands), the original post has vanished from my view. Perhaps it still exists somewhere, but for the life of me, I cannot find it; therefore, I post to the next most appropriate place to describe and discuss this shortcoming: This blog.
Perhaps my next blog post will be a discussion of keeping blogs on this site...
I had wanted to include an example of using
define() to define an operator
&d to calculate the (Euclidean) distance between two points (given as lists). This simple procedure turned out to be extremely complicated, since define() does not seem to have been designed to work with anything more than basic objects/structures. Let me walk you through what I was trying to do and what's really going on.
My first attempt to define this operator looked something like this:
define( '`&d`' , 'orderless', 'conditional'(&d( a::list, b::list ) = simplify( sqrt( add( (a[i]-b[i])^2, i=1..nops(a) ) ) ) , nops(a)=nops(b) ) );
When I would test my operator
&d, the results were consistently incorrect.
[1,2,3,4,5] &d [4,3,2,1,0]; would produce
A lot of head-scratching and a little digging later, infolevel shed some light onto the situation: By setting infolevel[all]:=3; , I could see that the pattern being matched/placed was simply
sqrt( (a-b)^2 )
My apologies for not knowing how to get nice Maple results pretty-printed in here. I guess I should eventually learn... Also, my apologies for not bothering to upload a file related to this post.
Why was this? Only one term in the sum? After a little more beating Maple with a stick, it occured to me: Maple is doing a simplification to the rule for &d when the define() command is run. Looking back at my
define() command and noticing that a and b are simply names in that definition, it's obvious why there's only one term under the radical: a is just a name, and
nops(a) returns 1. Maple is evaluating far, far too early.
So then, of course I wanted to get around this problem. I tried to slow down the evaluation by creating multiple "helper operators." Looking at the process of finding the (Euclidean) distance between two points, the first thing to do is to subtract the points pairwise. Then, one must square each difference. Why not just have &d do the subtraction and pass that list (as Maple can natively perform element-wise subtraction of one list from another, or whatever you may call it) to another operator &e which does element-wise squairing? This did not work too well:
define( '`&d`', 'orderless', &d( a::list, b::list ) = &e(a-b) ) ;
define( '`&e`', &e( a::list ) = map(`^`, a, 2) ); #Square each element of a list
[1,2,3,4,5] &d [4,3,2,1,0]; produces
Again, the premature evaluation had thwarted me; a was just a name, and the map command produced a^2. Try as I might to build auxiliary operators, it was just too much to try to work with an abstract list!
You will notice that from this point on, there is no more checking for equal-length lists. I discarded that part of the operator because Maple will throw an error if one tries to subtract one list of a particular length from another list of a different length.
Yesterday, while hacking The Mathematics Survival Kit ebook, I learned about the lovely command
eval(x,n::posint);, which evaluates the expression x to only n levels. Perhaps this would help. I tested it out on a simple operator:
define( '`&f`', &f( a::list ) = eval( map(`^`, a, 2), 1) ); #Square each element of a list
Success! In a most unintuitive move, one must force Maple to delay the evaluation in order for an operand defined via
define() to process objects which are not simple constants.
Therefore, in order to use
define() to make an operator which takes two lists and returns the (Euclidean) distance between the two points (points in Euclidean n-space) represented by those lists, one must do something similar to the following:
define( '`&d`', 'orderless', &d( a::list, b::list ) = eval( simplify( sqrt( `+`( op( map( `^`, a-b, 2 ) ) ) ) ), 1 ) );
Just to note at this moment: I have been specifying &d to be an "orderless" operator. I am fully aware that this is unnecessary, as &d(a,b)=&d(b,a) by the definition of the (Euclidean) distance function. I am specifying "orderless" so that if one were to change the definition of &d to allow for symbolic input and not just lists, then the commutativity would be preserved.
Finally, I should note one interesting tidbit I picked up along the road to solving this problem. Instead of a helper operator &e, I define a procedure Q by:
Q := proc( a::list )
If I were to define my &d in this way:
define( '`&d`', 'orderless', &d( a::list, b::list ) = Q(a-b) );
define() would cause an error, because it is trying to pass the object
a-b to the procedure Q, which is expecting a list. If, instead, I were to define my &d in this way:
define( '`&d`', 'orderless', &d( a::list, b::list ) = `Q`(a-b) );
then the evaluation of Q is held until &d is evaluated. Using backticks to delay the evaluation of a "helper operator" works exactly the same way; the problem is that the backticks do not delay the evaluation of the
map() command (along with, I imagine, many others).
Let me make one final note: DO NOT REPLY TO THIS POST TELLING ME HOW I COULD HAVE DONE ALL OF THIS USING PROCEDURES!! The entire point of this post is to discuss the unintuitive shortcomings of the define() command. If you would like to do the same thing by defining a procedure, then all you have to do is:
`&d` := proc( a::list, b::list )
return simplify( sqrt( `+`(op(map(`^`,a-b,2))) ) );