Consider the task of sorting a list of complex floating-point numbers by magnitude.
First Attempt
The usual method to do this in Maple is with the sort procedure. By passing a boolean-valued function that computes then compares the magnitudes of two complex numbers, we can sort the list. The following procedure shows how this is accomplished.
sort1 := proc(L)
return sort(L, proc(z1,z2) abs(z1) <= abs(z2) end proc);
end proc:
A disadvantage of this approach is that the absolute-value procedure is called twice every time a pair of numbers is compared. For a long list, the time spent in the absolute value routine dominates the computation time.
Second Attempt
To avoid the overhead of repeatedly computing absolute values, we can first compute the magnitude of each number. Because we want the sorted list of complex numbers, not their magnitudes, we need a way to associate a magnitude with its corresponding complex value. The usual technique is to pair the two in a sublist. The list of pairs are sorted by passing a comparison procedure that extracts and compares the magnitude from each pair. After sorting, the complex values are unpacked from the list of pairs. The following procedure implements this technique in a minimal fashion.
sort2 := proc(L)
local l;
return [seq](l[2]
,l = sort([seq]([abs(l),l], l=L)
,(p1,p2) -> p1[1] <= p2[1]));
end proc:
A drawback to this procedure is that the comparison procedure is user-defined, hence slower than the built-in sort comparison.
Third Attempt
We can sort a list of floats (the magnitudes) with the builtin sort comparison procedure 'numeric' (also known as `<`). However, we need a way to associate each magnitude (float) with the specific complex value that generated it. The setattribute and attributes procedures provide the desired capability. Use setattribute to set the original complex number as an attribute to the computed magnitude. Sort the list, using the builtin `<` comparison procedure, then use the attributes procedure to convert the sorted list of magnitudes to the corresponding complex values. The following procedure implements the technique.
Note that this does not work with a list of Gaussian integers (complex numbers with integer components); setattribute only works with some expression types. It works with floats, but not integers, complex integers, or complex floats. To ensure that there are no problems, we should consider declaring the parameter L to be a list of complex floats (L::list(complex(float))).
sort3 := proc(L)
local l;
return map(attributes,sort([seq](setattribute(abs(l),l), l=L),`<`));
end proc:
Verification
Verify that each procedure sorts properly.
n := 10:
L := [seq](RandomTools:-Generate(complex(float(range=0..1
,method=uniform)))
,i=1..n):
sort1(L);
sort2(L);
sort3(L);
Timing Comparison
Time the three alternatives with a fairly long list, 10,000 complex floats.
n := 10^4:
L := [seq](RandomTools:-Generate(complex(float(range=0..1,method=uniform)))
,i=1..n):
time(sort1(L));
time(sort2(L));
time(sort3(L));
20.650
2.510
1.030
Both sort2 and sort3 are significantly faster than sort1, the naive implementation. Using setattribute more than doubles the speed again. For shorter lists the speed increases are less dramatic, the time depending on the number of calls to the comparison routine. For a merge sort, which the Maple sort routine uses, the number of comparisons for a list of n elements is
add(ceil(evalf(log[2](3/4*k))), k=1..n);
See Knuth, Donald, "The Art of Computer Programming," Vol 3, 2nd Edition, p. 186.
Suggestion
While the setattribute trick works with this example, it is unfortunate that it cannot be more generally applied. It might be useful if the syntax of sort were extended so that we can use the built-in comparison routine to compare particular elements of sublists. For example, it would nice to write
sort4 := proc(L)
local l;
return [seq](l[2]
,l = sort([seq]([abs(l),l], l=L),`<`[1]));
end proc:
where the `<`[1] tells the builtin function to numerically compare the first item in a list of lists.
Comments
Magnitude
Joe,
It is a great idea to use attributes.
Replacing
abs(x)withop(1,x)^2+op(2,x)^2seems to give some time improvement in all examples. It can't be applied to real x though.faster magnitude comparison
Another possibility, which works with pure reals and imaginaries, is (Re^2+Im^2)(z), or a variant. Here is a comparison of several methods for generating the magnitude (or square of magnitude).
restart; use RandomTools in zs := [seq](Generate(complex(float(range=0..1,method=uniform))),i=1..10^3); end use: n := 2; time(proc() to 10^n do map(z -> z*conjugate(z),zs) od end()); time(proc() to 10^n do map(Re^2+Im^2, zs) od end()); time(proc() to 10^n do map(z -> Re(z)^2+Im(z)^2, zs) od end()); time(proc() to 10^n do map(abs,zs) od end()); time(proc() to 10^n do map(z->op(1,z)^2+op(2,z)^2,zs) od end()); time(proc() to 10^n do [seq(Re(z)^2+Im(z)^2, z=zs)] od end()); time(proc() to 10^n do [seq(op(1,z)^2+op(2,z)^2,z=zs)] od end()); 7.570 6.259 6.879 12.050 7.079 5.410 5.610Using seq is somewhat faster than map. The fastest here, Re(z)^2+Im(z)^2, is also robust, it works for pure reals and imaginaries.
seq vs. map
One seq seems to be faster than map. For 2 seqs, the situation may be different - I am not sure. For example, multiple timings of the following 2 sorting procedures didn't show a significant preference of one of them,
Re: Magnitude
It can't be applied to real x though.
Nor to purely imaginary x either.
This trick is awesome.
It took me a while to realize just how useful this trick really is. If you want to sort integers you can write an integer n as SFloat(n,0). I think that has pretty low overhead. Needless to say, I am now using it to sort monomials by degree. It is substantially faster than accessing an element of a list. Thanks for the great technique.
Nice extension
Using SFloat to convert an integer into an "equivalent" software float, so that it can take an attribute, is a nice extension. Thanks for sharing it. For those wondering what Roman is doing, here is the idea (I assume),
sortpolys := proc(polys,vars) local p; return map(attributes ,sort([seq](setattribute(SFloat(degree(p,vars),0),p) ,p = polys))); end proc: sortpolys([x^3+x, 1, x^2+x+1],x); 2 3 [1, x + x + 1, x + x]However, I wouldn't think that there would be much speed advantage of this compared to my method 2 (where the keys are precomputed but pairs are used for the comparison) except for fairly long lists. How many monomials do you have to sort?
Is there a similar method for making any numeric quantity attributable while permitting it to be directly compared with other such numbers? SFloat can handle floats but not fractions:
SFloat(2.3, 0); 2.3 SFloat(2/3, 0); Error, invalid arguments for Float constructorI'd be reluctant to use evalf because of the possibility of round-off.
Thanks again
> How many monomials do you have to sort?
1-2 thousand, in a loop :) This is a program I have been optimizing for quite some time. Your tip gave a notable improvement, and now I think I should release the program. I am indebted to you for this great technique, and you are credited in the source code.
> Is there a similar method for making any numeric quantity attributable?
I have no idea. I think you would have to use evalf.
Thanks again,
Roman
A few thousand monomials is
A few thousand monomials is more than I would have guessed, and enough to have an effect, as you mention. It's neat that the technique has a useful application---my write up on sorting complex values was merely a demonstration vehicle though a related problem prompted the "discovery". Over the years I've looked for a good use for attributes, but hadn't recognized any until now. They are used a fair amount in the Slode package, but I haven't looked into how.