@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.

Here is the procedure for splitting the hex string into the three fields: sign, exponent, and significand (aka mantissa). The info is returned in two forms: In the nested record Raw, each bit string is simply converted to a nonnegative integer. In the enveloping record, the info is presented in a natural form of signed integers S and E such that the represented number is S*2^E. The conversion between the two forms is that from the raw form, the represented number is (-1)^b*(1+s/2^52)*2^(e-1023) where b is the sign bit, s the 52-bit significand, and e the 11-bit exponent. The cases e=0 and e=2^11-1 are special cases. Note that the raw exponent is never negative.

The procedure can also be passed an hfloat, which'll first be converted to a hex string by the previous procedure.

SplitHF:= proc(H::Or(hfloat, And(string, 16 &under length)))
description 
   "Split a double-precision hardware float (IEEE-754 standard), or corresponding"
   " hex string, into sign, exponent, and significand fields."
;
option 
   `Author: Carl J Love <carl.j.love@gmail.com> 8-Feb-2019`,
   `Reference: https://en.wikipedia.org/wiki/Double-precision_floating-point_format`
   #Reference documents all bit-position magic numbers.
;
local raw_sg, raw_sf, raw_ex, sg, sf, ex;
    (raw_ex, raw_sf):= sscanf(`if`(H::hfloat, `hf<->hex`(H), H), "%3x%13x")[];
    raw_ex:= irem(raw_ex, 2^11, raw_sg);
    sg:= (-1)^raw_sg;
    (ex,sf):= (raw_ex-2^10-52+1, sg*(2^52+raw_sf));
    Record(
       "mantissa"= sf, 
       "exponent"= ex, 
       "rational"= 
          sg*
          `if`(
             #The exponent being all 0 or all 1 bits are special cases:
             #   +/- 0, +/- infinity, undefined (NaN), underflow, overflow. This
             #   code doesn't distinguish those latter three, and +/- appear alike.
             #   These distinctions can be made by examining the Raw form.
             raw_ex in {0, 2^11-1},
             `if`(raw_sf = 0, `if`(raw_ex = 0, 0, infinity), undefined),
             Scale2(sf, ex)
          ),
       "Raw"= Record("signbit"= raw_sg, "significand"= raw_sf, "exponent"= raw_ex)
  )
end proc
:     

@radaar No, there's no limit on "memory used". It can be several times the amount of memory on your system without causing any issue. Here's an analogy: You can take a long trip in your car, using 100 liters of fuel, without any issue, even though your fuel tank only holds 40 liters. The number is only useful for comparing two programs. It does not represent the amount of memory in use at any one time.

Anyway, the memory numbers reported for this particular program are small. It's nothing to worry about.

@emendes There was a bug in my procedure, now fixed. Please download the code and try again.

Regarding simulation: I think that it'd be fairly easy to simulate basic four-function arithmetic using objects to represent the floats. What is your purpose? I assume that it's pedagogical. Are you trying to teach how round-off errors can lead to anomalies? The object floats would be represented as pairs of integers (mantissa and exponent). 

I'll start writing the splitting code.

@mmcdara Vote up (largely for your Reply immediately above).

I believe that the decline in efficiency is precisely due to the reason that you said: It's "stupid" (efficiency-wise) to use range= -1..1 when the standard deviation is small; that range is too wide. Also, it's stupid (accuracy-wise) to use range= -1..1 when the standard deviation is large; that range is too narrow. In particular, if range= a..b and either CDF(X, evalf(a)) = 0.0 or CDF(X, evalf(b)) = 1.0, then the range is too wide. Ideally, a should be slightly >= Quantile(X, eps) and b slightly <= Quantile(X, 1-eps) where eps is akin to a "machine epsilon": a small positive number such that evalf(1 - eps) < 1 (strict inequality!) and such that neither Quantile is +/- Float(infinity).

Regarding "adaptive": It's adaptive to an extent; GIGO applies. If you want it to be adaptive in this respect, then don't specify the range; it will then choose a suitable one. Not specifying the range is equivalent to specifying range= deduced. That having been said, for this example there will likely be cases where it fails to calculate an appropriate range and requests that the user does it.

@FrancescoLR It seems to me that to achieve any real benefit from parallelism, you'd also need multiple disk drives operating in parallel. Otherwise, I'd think that the disk drive rather than the CPU would be the bottleneck. I'm not sure about this, because perhaps buffering the file output to RAM ameliorates the bottleneck effect.