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

@Magma Okay, granted, my comment was impolite. Since my primary purpose is that readers learn from my posts, how can I help you to learn efficient Maple, and reasonable coding practices (regardless of language)?

If you want to use GF, much better code would be

select(
   p-> andmap((g,F)-> F:-order(F:-input(g)) = p-1, [2,3], GF(p,1)),
   select(isprime, [seq(2001..2999, 2)])
);

And even if you want to go old-school with a while loop, either of these is far better than what you had:

p:= 2000: 
L:= table():
while (p:= nextprime(p)) < 3000 do
   G:= GF(p,1);
   if G:-order(G:-input(2)) = p-1 and G:-order(G:-input(3)) = p-1 then
      L[p]:= p
   end if
end do:
[entries(L, indexorder, nolist)];
#
# Or
#
p:= 2000: u:= 0:
L:= Vector():
while (p:= nextprime(p)) < 3000 do
   G:= GF(p,1);
   if (G:-order@G:-input)~({2,3}) = {p-1} then
      L((u:= u+1)):= p
   end if
end do:
[seq(L)];

I thought that your followup entitled "not working..." was appropriate as a separate Question, so I made it so, with the title "[Finite-field arithmetic] not working..." Please continue the discussion of that issue there. And in this thread, please continue discussing how to use Maple to work with the quotient ring GF(p,k) / f where k > 1 and f is a reducible univariate polynomial over GF(p,k).

@Magma That's ridiculous. It seems that you've learned nothing from my Answer.

@bellaqueame I've rewritten your ElemRing:

ElemRing:= proc(
   f::depends(polynom(polynom(integer, tauh), x)), 
   x::name,
   tauh::specfunc(polynom(integer, _Z), RootOf),
   p::prime
)
local i, n:= p^(degree(op(tauh))*degree(f,x))-1, ring:= rtable(0..n);
   for i to n do 
      ring[i]:= Nextpoly(ring[i-1], x, tauh) mod p
   od;
   [seq(ring)]
end proc: 

Please execute this in Math Input (aka 1D Input) or in a Code Edit Region because I suspect that the depends parameter modifier doesn't work in 2D Input. Your versions have numerous levels of unnecessary, nested list-forming brackets. I think that this is the reason that they don't work. I don't feel like disentangling all your nested lists, because they're simply not necessary.

Note that I changed the parameters (and their order) from your ElemRing. Like before, the x (a global unassigned name) should be passed! It's not necessary to pass k because it can be deduced from tauh. The order is just my style preference and follows the argument order of Nextpoly. It is traditional in Maple that the field characteristic p be the last parameter.

@bellaqueame From your "invalid power" error, I assume that you're using 2D Input. The error can be corrected by changing

<seq(a^(` q `^k), k= 0..n)>^+

to

<seq(a^(` q `^k), k= 0..n)>^%T

and also make the same change on the following line.

I strongly discourage the use of 2D Input. I use Maple Input (aka 1D Input).

I haven't looked at your code yet.