@Markiyan Hirnyk 

Here you did not provide an order for the variables, so a the program has chosen [ю, щ] for you
and has decomposed the solution set into 7 disjoint cells.

@Markiyan Hirnyk 

RegularChains is used internally by SemiAlgebraic . You need RegularChains only for special manipulations.
As I said, the order of the variables determines the type of the cells. I think that it is clear.

@Markiyan Hirnyk 

As solid people I can say that the first CPV is 3.80259236620486251646441207770...
and the second one is - infinity.

@Markiyan Hirnyk 

Interesting examples. It is possible to compute them using series; direct computations as in your answer would need a very large value for Digits and probably a prohibitive amount of time.

@Markiyan Hirnyk 

This way for a small variation like

f:=exp(-sec(t))*cos(t)/(-1/4+sin(t)^2) + 1/sin(t-Pi/6)^5;

is complicated.

@Carl Love 

The simplification of floor is very weak. Even
simplify(floor(x)) assuming x<1,x>0;
remains unevaluated.

This is not a big problem because the isolation of the singularity usually works.

Unfortunately there is a bug in series which I did not mention in the answer.

series(1/(sin(x)),x=0, 1);
Produces the same error.

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.