AFAIK Maple does not approximate CPV integrals.

@Mariusz Iwaniuk 

@Carl Love 

1. I mentioned d+20 just for safety (e.g. any modern processor also has an extended precision (80 bits) ).  My point was that using hardware floats a hex digit of Pi could be impossible to determine using BBP even for rather small d; for example if this digit is an F and it is surrounded by many F's.  I agree that such situations should be rare, but a warning is OK.

2. In your SSj the hfloats are not used in the first sum which is actually the bottleneck, so increasing Digits by 3 has practically no impact.

@Axel Vogt 

Such situations are actually inherent! I am aware that BBP may miss some trailing hex digits for  some d's, but it was strange that this happened for small d's.

Assume that starting at some position d, the hex digits of Pi are:
3FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF1234...
(such a d should exist because Pi is probably a normal number).
Now, using hfoat computations, BBP(d+20)  could be any of  FFFFFFFFF  or 000000000
depending on how the round-off errors are combined.
So, e.g. the hex digits  d+20, d+21, ...   simply cannot be computed/guaranteed with BBP and hfloats!

All these are special (particular) cases.
It is not possible to express the general solution in terms of real functions.
Just think at the much more simple case of the Laplace equation  diff(f(x,y),x,x) + diff(f(x,y),y,y) = 0.
Or, even better, at the Cauchy-Riemann equations: 
diff(f(x,y),x) - diff(g(x,y),y) = 0, diff(f(x,y),y) + diff(g(x,y),x) = 0.

@Carl Love 

The 9th hex digit is not always correct (I have checked this). I'd recommend to return only 8 hex digits in the procedure BBP (i.e. ModuleApply).

Edit: unfortunately the 8th digit is also sometimes incorrect. Please revise the code.

Edit2: I think I found the problem. ModuleApply must have

Digits:= trunc(evalhf(Digits)) + 3;

Note that +2  is not enough.
Now the results are correct at least for d <= 5000.
The 9th digit is still wrong.
I am not sure about d>5000  :-(

@Carl Love 

Nice! Vote up!

P.S. Why do you say that  16 &^ 0  mod  4;   is a bug?  It is indeed  0 &^ 0  mod 4; 
(so, a convention, even if 0^0 = 1).

The fact that the hypercube [0,1]^d is densifiable is obvious: take the (m+1)^d points having all coordinates in the set {0,1/m,2/m,...,m/m} and take any polygonal curve passing thru all of them. With this approach, the "t" can be found at once.

Of course, the "cosine" curve g is more elegant, but the problem is that the "proof" you mentioned is only a sketch. It is not very hard to obtain a real proof such that given any of the m^d small hypercubes H, find (effectively) t in [0,1] for which g(t) \in H.
This is the way to go, and the Maple code is not needed then.

So, you have a kind of "almost" space filling curve  (see wiki).
- What is the source for this curve?
- Do you know to prove mathematically the claimed property? Or you want some kind of "incomplete" proof using Maple?
Or maybe you want to use this curve for something else?

@tomleslie 

It is not difficult to explain.
When SolveEquations solves f(x)=0, it actually minimizes an associated function F(x)  (such as F(x) = f(x)^2).

A minimization with constraints g(x)<=0  is much more difficult (mainly for several variables). It is not possible in general just to minimize F and then select the solutions satisfying g(x)<=0  because in most cases the set of such solutions could be empty. It is faster for the user to select himself the desired solutions if he knows that such solutions exist and they are not too many.

@Earl 

You have mentioned the tangency. In this case of a not differentiable curve it should be interpreted as having at least a common point and disjoint interiors. This can be done but it's of course more complicated.

So, without shades of grey and without thickness? Then make a custom legend using textplot, line, pointplot.

@weidade37211 

@Earl 

If you really want exact tangency with a circle  you should solve the problem separately.
Since the curve is not differentiable for t=k*Pi (k in Z), these points need special atention, fsolve/DS only is not enough.

@tomleslie 

Combining yours with mine:

r1:=t->piecewise(t<Pi, sqrt(exp((t+(1/4)*Pi)*(1/2))),
                t<2*Pi, sqrt(exp((9*Pi/4-t)*(1/2)))
               ):
r:=t->r1(t - 2*Pi*floor(t/(2*Pi))):
plots:-animate(plot, 
[[[r(t-u)*cos(t),r(t-u)*sin(t),t=0..2*Pi] ,
[[0,r(Pi/2-u)*sin(Pi/2)],[0,2+r(Pi/2-u)*sin(Pi/2)]]   ],thickness=[1,6]],
u=0..2*Pi, scaling=constrained);

@Adam Ledger 

For this specific s()  it does not matter too much; it's only a shift of the sequence with one position. 
Of course for a recurrence  s(n) = f(n, s(n-1)) rather than  s(n) = f(s(n-1)),  the starting index could be important.