@bsoudmand I get manifestly real expressions (with Re and Im resolved). My guess is that you have an earlier version of Maple such that evalc won't automatically "map over" my fancily labelled expressions. Try this:

simplify(evalc~([Re,Im](Ec))) assuming positive;

Kitonum's Answer shows that evalc is not necessary in this case (a rational function is an easy case). However, in general I think it will be necessary.

@acer I was wrong. I was misled by some spurious timings that I can't reproduce.

Instead of cot~(...), you could use map[evalhf, inplace](cot, ...). It would be much faster, but your point is still valid.

@emendes I just posted, as a separate Answer, a cycle detector that works over the exact algebraic numbers.

@emendes Okay, here's an efficient equivalent to Mathematica's NestList:

NestList:= proc(f, x, n::nonnegint)
local R:= rtable(0..n, [x]), k;
   for k to n do R[k]:= f(R[k-1]) od:
   [seq(R)]
end proc:

I also just posted (two Replies above) some alternatives that don't use for but require Maple 2018+ for the embedded assignment.

@emendes It is difficult (although not impossible) to effectively use seq for a recursive evaluation, i.e., for creatiing a sequence whose terms are computed from previous terms. Why don't you want to use for?

Here's a very easy (and very inefficient) code that does it using seq:

f:= y-> 4*y*(1-y):
x0:= (5-sqrt(5))/8:
seq((f@@p)(x0), p= 1..9);

In Maple 2018, that can be improved with embedded assignment:

f:= y-> 4*y*(1-y):
x:= (5-sqrt(5))/8:
seq((x:= f(x)), p= 1..9);

Or, noting that p is now irrelevant, seq can be replaced with $:

f:= y-> 4*y*(1-y):
x:= (5-sqrt(5))/8:
'(x:= f(x))' $ 9;

Syntax note: An embedded assignment always requires its own parentheses, even if they seem redundant or there's no possibility of ambiguity.

@vv That works in this case, which is a simple polynomial. If you step it up even one notch, to a rational function, and use a high value of p, I suspect that you'd quickly either consume all your memory or cause a kernel failure.

Or, f could be purely numeric, for example, based on a dsolve(..., numeric) result.

Regarding unary :-: I thought that Acer wasn't as clear as possible about this point, so I'm rephrasing it. When any symbol is prefixed with :-  with no left operand for the :-, it refers to the global instance of that symbol, not to any module or package member or procedure local or top-level local that may otherwise have the same name.

Regarding with: I strongly recommend losing the habit of using it. For me, there are only two acceptable cases for using it:

1. When using Maple as a "desktop calculator" for work that will not be saved.
2. When using an older package (such as Tolerances) that rebinds operator symbols such as `+`. Newer packages can (and should) provide the same functionality by using objects and overloading the operators.

A desire to use "the short-form name" is not an acceptable reason to me. It sickens me whenever I see that mentioned in a help page. Using short-form names leads to hard-to-read code because when using multiple packages, the reader can't easily tell which package (if any) the name comes from.

@Preben Alsholm It's probably worth noting that the average of the two is the well-known Trapezoid Rule, in this case applied to potentially unevenly spaced ordinates.

@Carl Love Acer has noted a bug with *~ when used with Tolerances: It uses the global * rather than the overloaded one. While my solutions above are still correct for this example, in light of this bug, the following will be safer:

subsindets(Z, {:-`*`, :-`+`}, e-> `if`(e:::-`*`, `*`, `+`)(e));

@acer You wrote the following code:

restart;
c:=piecewise(u+v<Pi, 1/6, u+v<2*Pi, 2/3, 1):
opts:=coords=spherical,color=c,grid=[49,49],style=surface:
plots:-display(plot3d([1,u,v],u=0..2*Pi-2*v,v=0..Pi,opts),
               plot3d([1,u,v],u=2*Pi-2*v..2*Pi,v=0..Pi,opts));

Your point being that using MESHes that are nonrectangular wrt their coordinate system can ameliorate the "sawtooth" effect, requiring a much smaller grid than that required to do it on rectangular MESHes. I point this out because it wasn't obvious to me from just reading your code; I needed to execute it.

I think that that point is better illustrated by seeing the mesh: Just remove style= surface. (Although the plot per se looks better without seeing the mesh.)

@acer Acer, you're correct that I overestimated the capability of VV's solution. It does only work for zgradient. So, the only function that can be displayed with continuous color variation over anchor colors is (x,y,z)-> z, which is a bit silly, since the values of that function are already evident on the plot. Also, I was wrong that the large grid is needed; it works just as well with a small grid.

Vote up for the way that you asked the Question. I very often vote up Questions without comment, but I want to give extra-special praise to your style.

@acer Vote up. The idea of replacing one COLOR structure with another, the latter created with colorscheme, is ingenious. Ever since colorscheme was introduced, I'd been trying to come up with an easy way to show continuous color variation with a small set of "anchor" colors (the [black, red, yellow, white] in this case). Your method is the way. VV also shows a way, but it comes at the expense of a 300x300 grid. While a grid of that size poses no problem to the kernel, it's often a problem for the GUI, which becomes sluggish, and that sluggishness remains as long as the plot is in the worksheet, even when off-screen. 

It's always been remarkable to me that the plot renderer is capable of showing continuous color variation with the variation seeming to be at the pixel level, regardless of the grid.

@denbkh Sorry, I misinterpretted your usage of "..." and your usage of "if".^[*1] You meant that n is required to be a power of 2 (in order for the expression to be defined at all); whereas, your "if" led me to think that n being a power of 2 was being considered as a special case. In that case, we're dealing with a finite indefinite sum whose last term is 1. Then the sum is indeed 2*n-1, as you said.

The exact result 2*n-1 can be obtained from this rsolve:

RS:= rsolve({S(0)=0, T(1)= n, T(x)= T(x-1)/2, S(x)= S(x-1)+T(x)}, {T(x), S(x)});
simplify(eval(S, eliminate(eval(RS, {S(x)= S, T(x)=1}), {S,x})[1]));

The only difference from my previous  rsolve is that the initial partial sum is 0, not 1.

^[*1]I think that the misinterpretation is largely my fault.