@Carl Love 

BTW, I always had the impression that op(L) is faster than L[] for lists, but it seems that the latter is a bit faster.

P.S. Do you have an example where a local for the variable in seq is mandatory? (probably only for multiple kernels [or threads?] in some circumstances).

P.P.S. (related to this thread).  I wanted to see the code for combinat:-choose. When called with 1 argument L, it returns Choose1(L).

I was not able to identify Choose1 (via opaquemodules etc) . Do you have any idea about it?
The life was easier without modules and objects :-)

@Carl Love 

Yes, thank you, I have corrected.

@Carl Love 

But for this problem you don't need arbitrary-length arithmetic. A simple routine to add two lists containing the digits is enough and it will be faster; the reversing is automatic in this approach.

@Carl Love 

For such problems needing computer speed and memory rather than the symbolic power of Maple,
a simple C program seems to be more suitable.

The database of exact solutions for Einstein's equations seems to be (almost?) 100% done.
It would be interesting an estimate (%) for the material contained in Abramowitz and Stegun/NIST
(or maybe Gradshteyn & Ryzhik). 

.

What kind of decomposition do you hope for this example? A factorization of the denominator is needed and this implies the complicated roots of a cubic (supposing that the indeterminates are x,y).

@one man 
My point was (is) that a "true" spiral on a cylinder does not have self-intersections.
You can convince yourself:
with(plots):with(plottools):
display(cylinder([0,0,0], 1, 6, style=surface, color=yellow),
spacecurve([cos(s/20*cos(s)),sin(s/20*cos(s)),3+s/20*sin(s),s=0..8*Pi],
color=red,thickness=5),scaling=constrained,axes=none);

Animated version:

@one man 

In Mathematics the definitions are clear. If there are non-equivalent definitions for a notion, we must say which one is used.
For example, here is something that is not spiral in the sense of "monotone curvature of constant sign".
(You should say what definition is satisfied by your spiral if your animation refers to a spiral). 

r:= t + 4*sin(t/2):
plot( [r*cos(t), r*sin(t), t=0..8*Pi]);

@one man 

Are you sure that you have a "true" spiral here? According to the definition it should have a monotone curvature of constant sign.

@acer 

It seems that codegen[optimize] is not always efficient. For

E3:=(x^2-x+2-y)^25+x+(1+y)^25:

it does not work.
[  (1+y)^25 annihilates the collect "trick" ].

@acer 

It is interesting to see how/when simplify/size is able to shrink an expression (polynomial here):

restart;
n:=75:
E1 := add(i*x[i],i=1..n)*add(i^2*x[i],i=1..n)+add(x[i]^2,i=1..n):
E2:=(x^2-x+2)^30+x+1:
E3:=(x^2-x+2-y)^25+x+1+y:
E:=[E1,E2,E3]:
                          
for e0 in E do
e1:=expand(e0):
e2:=simplify(e1,size):
print(length~([e0,e1,e2]));
od:

                      [3009, 84700, 44505]
                        [34, 1287, 1287]
                       [42, 15924, 6938]

@TomM 

Such a phenomenon is normal when working with floats.
If you convert the foats to exact rationals (i.e. 0.333 becomes 333/1000), then the coefficient of x^9 happens to be exactly 0
(btw, why didn't you rename your vars for us?, a search & replace is so simple ...).

Using floats, the coeff of x^9 may of course be a (usually) small nonzero (due to cancelations, catastrophic or not).
But in practice when a float such as -53417.4875669262 appears as input, it usually means that the actual quantity could have been -53417.48756692623... or -53417.48756692618... and in such a case the coefficient of x^9 could be of course <> 0. It is very simple!

@Kitonum 

An unpleasant fact is that the plot3d routine is not able to draw this cylinder properly
when z is given explicitely.

expr:=(x-y)^2+(y-z)^2+(z-x)^2 = 3:
z1,z2:=solve(expr,z):
plot3d([z1,z2],x=-3..3,y=-3..3,grid=[100,100], scaling=constrained);

plot3d([z1,z2],x=-3..3,y=x-sqrt(2)..x+sqrt(2),grid=[100,100], scaling=constrained);

@Bachatero 

Unfortunately this is theoretically impossible due to Richardson's_theorem.

Note that even now Maple is not useless here:

E:=-n*mu+ln(product(mu^x[i]/factorial(x[i]), i = 1 .. n)):
E4:=eval(E,n=4):
extrema( E4,{},mu,'sol'):sol;


It remains to extrapolate ...