BBP算法所需的工作精度是多少？
我想在内存有限的环境下计算圆周率的第 n 位数字。因为我不能使用小数，所以这个只用整数的BBP算法在Python中对我来说是个很好的起点。我只需要一次计算一个圆周率的数字。我该如何确定可以将D设置为最低值，也就是“工作精度的位数”？
当D=4时，我能得到很多正确的数字，但有几个数字会错一位。例如，计算第393位时，精度为4给我的是0xafda，从中我提取出的数字是0xa。然而，正确的数字是0xb。
无论我把D设置得多高，似乎总能找到一个测试的数字，公式返回的值都是错误的。
我尝试在数字“接近”另一个数字时提高精度，比如0x3fff或0x1000，但找不到“接近”的好定义；例如，计算第9798位时，我得到的是0xcde6，这和0xd000并不是很接近，但正确的数字是0xd。
有没有人能帮我弄清楚使用这个算法计算某个数字需要多少工作精度？
谢谢，
编辑
供参考：
precision (D) first wrong digit ------------- ------------------ 3 27 4 161 5 733 6 4329 7 21139 8+ ???
请注意，我是一次计算一个数字，例如：
for i in range(1,n):
D = 3 # or whatever precision I'm testing
digit = pi(i) # extracts most significant digit from integer-only BBP result
if( digit != HARDCODED_PI[i] ):
print("non matching digit #%d, got %x instead of %x" % (i,digit,HARDCODED_PI[i]) )
2 个回答
在编程中，有时候我们会遇到一些问题，特别是在使用某些工具或库的时候。这些问题可能会让我们感到困惑，尤其是当我们不太了解这些工具的工作原理时。
比如说，有人可能在使用某个库的时候，发现它的某个功能没有按照预期工作。这时候，我们就需要去查找相关的资料，看看是否有其他人遇到过类似的问题，或者是否有解决方案。
在这个过程中，我们可能会在网上找到一些讨论，比如在StackOverflow上，其他开发者分享他们的经验和解决办法。这些讨论通常会提供一些代码示例，帮助我们更好地理解问题所在。
总之，遇到问题时，不要害怕去寻求帮助，很多时候，其他人已经遇到过类似的情况，并且愿意分享他们的解决方案。
from typing import TypeVar
from gmpy2 import mpz, mpq, powmod as gmpy2_powmod, is_signed as gmpy2_is_signed
__all__ = ['PiSlice']
Integer = TypeVar('Integer', int, mpz)
class PiSlice:
'''
References
----------
"BBP digit-extraction algorithm for π"
https://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula
'''
version = '1.0.0'
def __spigot(self, p: Integer, a: Integer, accuracy: mpq) -> mpq:
def search_junction(p: Integer, a: Integer) -> Integer:
n = mpz(0)
divisor = 8 * p + a
while 16 ** n < divisor:
n += 1
divisor -= 8
return p - (n - 1)
p = mpz(p)
junction = search_junction(p, a)
s = 0
divisor = a
for k in range(junction):
s += mpq(gmpy2_powmod(16, p - k, divisor), divisor)
divisor += 8
for n in range(mpz(p - junction), -1, -1):
if (intermediate := mpq(16 ** n, divisor)) >= accuracy:
s += intermediate
divisor += 8
else:
return s
n = mpz(1)
while (intermediate := mpq(mpq(1, 16 ** n), divisor)) >= accuracy:
s += intermediate
n += 1
divisor += 8
return s
def __init__(self, p: Integer):
'''
'''
self.p = p
def raw(self, m: Integer) -> Integer:
'''
Parameters
----------
m: Integer
Sets the number of slices to return.
Return
------
random_raw: Integer
Returns a hexadecimal slice of Pi.
'''
p = self.p
spigot = self.__spigot
accuracy = mpq(1, 2 ** (mpz(m + 64) * 4)) #64 is the margin of accuracy.
sum_spigot = 4 * spigot(p, 1, accuracy) - 2 * spigot(p, 4, accuracy) - spigot(p, 5, accuracy) - spigot(p, 6, accuracy)
proper_fraction_of_sum_spigot = mpq(sum_spigot.numerator % sum_spigot.denominator, sum_spigot.denominator)
if gmpy2_is_signed(proper_fraction_of_sum_spigot):
proper_fraction_of_sum_spigot += 1
return mpz(mpz(16) ** m * proper_fraction_of_sum_spigot)
无论我把D设置得多高，似乎只要测试足够多的数字，就会找到一个公式返回错误值的情况。
如果你测试的数字足够多，你总是会遇到错误，因为这个算法不支持任意精度，所以最终会出现四舍五入的错误。
无限循环直到数字不再变化，这样很难确定计算特定数字所需的最小精度。
最好的办法是通过实验来确定，理想情况下是和一个已知正确的来源进行比较，逐渐增加数字的精度，直到结果匹配。如果没有正确的来源，可以从你能达到的最大精度开始（我猜是14位，因为第15位几乎总是会有四舍五入的错误）。
编辑：更准确地说，这个算法包含一个循环 - 从0到n，其中n是要计算的数字。每次循环都会引入一定量的误差。经过足够多次循环后，误差会影响到你正在计算的最重要的数字，因此结果会出错。
维基百科的文章使用14位精度，这足以正确计算第10**8位数字。正如你所展示的，精度越低，错误出现得越早，因为精度不足，经过较少的循环就能看到误差。最终结果是，能够正确计算某个数字的n值在精度降低时会变得更低。
如果你有D个十六进制数字的精度，那就是D*4位。每次迭代，最不重要的位会引入0.5位的误差，因此经过2次迭代，最不重要的位可能会出错。在求和过程中，这些误差会被累加。如果累加的误差达到了最重要数字的最不重要位，那么你提取的单个数字就会出错。大致来说，当N > 2**(D-0.75)时，就会出现这种情况。（具体要根据某个对数基数来修正。）
根据你的数据进行经验推断，似乎有一个大致的关系是N=~(2**(2.05*D))，不过数据点不多，所以这可能不是一个准确的预测。
你选择的BBP算法是迭代的，因此计算序列中的数字会越来越耗时。要计算数字0到n，将需要O(n^2)的步骤。
维基百科的文章提供了一个计算第n位数字的公式，这个公式不需要迭代，只需要指数运算和有理数。这不会像迭代算法那样失去精度，你可以根据需要以常数时间（或者在最坏情况下是对数时间，具体取决于模运算的实现）计算出任意位的π，因此计算n位数字将需要O(n)时间，可能是O(n log n)。