# The following work is based on https://en.wikipedia.org/wiki/Bailey-Borwein-Plouffe_formula/
# Omar Aceval Garcia, 03/27/2025

# nthPiUnreduced Computes the n-th hexadecimal digit of pi. The usual digits are denoted 0,1,2,3,4,5,6,7,8,9,A,B,C,D,E,F
# An output of 15 for example is interpreted as F

nthPiUnreduced := proc(n) local k,portion:
portion := 4*add(16^(n-k)/(8*k+1), k=0..n) -
2*add(16^(n-k)/(8*k+4), k=0..n) - 
add(16^(n-k)/(8*k+5), k=0..n) - 
add(16^(n-k)/(8*k+6), k=0..n):
return floor(16*(frac(portion))):
end:

# The reason it is called "unreduced" is because originally I was reducing the 16^(n-k) terms modulo the denominator, 
# since only the fractional part of the big sum matters. But actually it turned out to be way faster to just let maple
# work with the big numbers. If you want, here is that version:

# If you want the original code I was playtesting, it's below.

# Quickly computing x^r modulo n
modexp :=proc(x,r,n) local p,k:
if n = 1 then return 0: fi:
p:= 1:
for k from 1 to r do
p := p*x:
if p > n then 
p := modp(p, n):
fi:
od:
return p:
end:

# It turns out frac(-0.1) spits out -0.1, AKA NOT the fractional part of -0.1, which is 0.9 (grr!)
realfrac := proc(x):
if x>=0 then return frac(x):
else return 1+frac(x):
fi:
end:

nthPi := proc(n) local k,portion:
portion := 4*add(modexp(16, n-k, 8*k+1)/(8*k+1), k=0..n) -
2*add(modexp(16, n-k, 8*k+4)/(8*k+4), k=0..n) - 
add(modexp(16, n-k, 8*k+5)/(8*k+5), k=0..n) - 
add(modexp(16, n-k, 8*k+6)/(8*k+6), k=0..n):
return floor(16*(realfrac(portion))):
end: