Many thanks to all who commented on my code, and especially to VV for his extensive testing for accuracy. Here is some improved code. The numeric computation runs in pure evalhf mode, so it's much faster, and it's probably more accurate also. I got 6 correct digits starting at the millionth position, the accuracy verified from the table at the end of the paper. I also simplified the conversion to a hex string to a single sprintf.

It should be pretty straightforward to convert this to Matlab.

SSj:= proc(d::nonnegint, j::posint)
local k;
   Frac(
      Frac(add(ModExp(16, d-k, 8*k+j)/(8*k+j), k= 0..d)) 
      + add(1/(8*k+j)/16^(k-d), k= d+1..d+14)
   )
end proc:

ModExp:= proc(b, N, m)
#Compute b^N mod m
local t:= 1, n:= N, r:= 1;
   while t < n do t:= 2*t end do;
   do
      if n >= t then r:= irem(b*r, m); n:= n-t end if;
      if t=1 then return r end if;
      t:= iquo(t, 2);
      r:= irem(r^2, m)
   end do
end proc:
      
Frac:= proc(x::numeric)  
local f:= x - trunc(x);
   `if`(f >= 0, f, 1+f)
end proc:

Pi_BaileyBorweinPlouffe:= proc(d::posint)
#Prints hexadecimal digits of Pi starting at an arbitrary position, d.
#9 digits will be printed because
#that's the limit of hardware-float arithmetic according to the paper.
#The accuracy of all 9 digits is not guaranteed.
local x, dm1:= d-1;
   Digits:= trunc(evalhf(Digits));
   x:= evalhf(16^9*Frac(4*SSj(dm1,1) - 2*SSj(dm1,4) - SSj(dm1,5) - SSj(dm1,6)));
   sprintf("%09X", round(x))
end proc:

What's the approximate size, n? What's the type of the coefficients (floats, integers, etc.)?

As far as I know, the only thing available for scheduling is Threads:-Sleep.

I'm not familiar with DifferentialGeometry, so I can't answer your question. I just wanted to point out that your posted file shows Maple giving a response to each of your commands. Perhaps those responses are wrong, but that's different than Maple not responding to your commands.

@briceM Sorry about that. The M -~ -1 needs to be changed to M -~ (-1). I already corrected the code in the Answer, so you can just copy-and-paste that.

Considering that you have numerous terms with the fractional exponents 2.41, 1.41, or 0.41, it seems likely that solutions are complex. If you remove the range restrictions of the variables and include the complex option, you'll get solutions.

@Lali_miani The set E is ordered, but the order is chosen by Maple for its own convenience. I'm guessing that your version of Maple is before Maple 12, and thus the order is based on addresses and may vary between sessions. Later versions of Maple use a set ordering that doesn't vary between sessions, but may vary between versions. You should never rely on the set order except to know that sets will be ordered, but in a way that you can't control.

@Kitonum Imagine looking up a word in a dictionary or a name in a phone book. I mean actual books, not something on a computer. Surely you'd make use of the fact that the words or names are sorted rather than looking at every entry. Checking that the entries are sorted would also be a huge waste of time. Rather, we must assume that they are sorted. If our search happens to prove that assumption false, maybe we report that or maybe we fail silently; that's not very important.

Being someone who answers many Questions here, I think that sometimes we should consider why the Question is being asked, what "bigger picture" is the asker trying to learn, what some unseen third party (such as an instructor) hopes the asker will learn, etc. That's what I mean by pedagogic spirit. After all, don't we teach the same courses as those unseen instructors?

@Kitonum I don't know if you are any sort of instructor or professor, as are many of the regular posters on MaplePrimes. If so, then surely you realize that while your Answer is technically correct, it violates the pedagogic "spirit" of the Question. To begin with, you've made no use of the fact that the list is sorted.

@audi001 You've failed to define dS__c and you've defined dS__a twice. Surely this was unintentional. Once you've set the correct dS__c, run the solve again. This time, there's a good chance that it'll run out of memory, that there'll be no solution, or that it simply takes too much time. In any of those cases, repost and we'll see what we can do. We'll start by setting some of the 10 parameters to 0 or 1 and trying to solve again.

I deleted your Question that was nearly the same as the one above.

"Warning, solutions may have been lost" is a common response from solve. It generally indicates that solve suspects that there are solutions, or that there are more solutions than it has found, but it can't find them due to limitations of its algorithms. It has nothing to do with running out of memory.

Please follow the instructions in my Answer below.

Hamid: While Mehdi has shown how to use if-then in the loop only because you explicitly asked about that, I must add that there's no good reason to do that.

@tomleslie Sorry, I don't know that. Doing the appropriate solve commands explicitly, I get the result whether I use decimals or fractions.

@Adam Ledger This is much simpler: floor(5^n/4) = (5^n-1)/4 for nonnegative integer n. Therefore, sum(5^n*floor(5^n/4), n= 1..N) = 5/96*(5*25^N - 6*5^N + 1). It's a shame that Maple can't automatically get that initial simplification for floor. The summation is trivial after that, being the difference of two geometric series.

The error is unrelated to what you show in the screenshot. We'd need to see the whole worksheet. Use the green up-arrow on the toolbar in the MaplePrimes editor to upload a worksheet. My first guess is that x1 and or x2 have been previously defined.

@sand15 No, Maple is not implicitly selecting a specific real interval for x or y when solve(sin(x)=y, x) returns arcsin(y). Both the domain and co-domain of arcsin are the complex numbers, as is true for nearly all of Maple's mathematical functions.