BBP算法所需的工作精度是多少？

Question

BBP算法所需的工作精度是多少？

7

我想在低内存环境下计算圆周率的第n位。由于我没有小数点，所以这个使用Python的仅整数BBP算法是一个很好的起点。我只需要一次计算一个数字。如何确定我可以设置D的最小值，即“工作精度的位数”？D=4给了我很多正确的数字，但是有几个数字会偏差1。例如，使用精度为4计算第393位数字，得到0xafda，从中提取出数字0xa。然而，正确的数字是0xb。无论我将D设置得多高，似乎测试足够数量的数字都会找到一个公式返回不正确值的数字。

我尝试在数字“接近”另一个数字时提高精度，例如0x3fff或0x1000，但是找不到“接近”的好定义；例如，在第9798位计算会得到0xcde6，它与0xd000并不非常接近，但正确的数字是0xd。 有人能帮我确定使用此算法计算给定数字需要多少工作精度吗？ 谢谢，编辑
参考：

精度（D）   第一个错误数字
-------------   ------------------
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]) )

- tba

2个回答

0

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)

- sosei

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- mdma · Accepted Answer

无论我将D设置多高，似乎测试足够数量的数字总能找到一个返回错误值的数字。如果您测试足够数量的数字，则始终会出现错误-该算法不使用任意精度，因此最终会出现舍入误差。当数字不变时，无限迭代将很难确定所需的最小精度。最好的方法是通过经验确定它，理想情况下是通过与已知正确源进行比较，并增加数字精度直到匹配，或者如果没有正确的源，则从您的最大精度开始（我猜测是14，因为第15位数字几乎总是包含舍入误差）。编辑：更准确地说，该算法包括一个循环-从0..n，其中n是要计算的数字。循环的每次迭代都会引入一定量的误差。在循环了足够多次之后，误差将侵入您正在计算的最重要的数字，因此结果将是错误的。维基百科文章使用了14位数字的精度，这足以正确计算10 ** 8位数。正如您所示，较少的精度会导致错误更早地发生，因为精度较低，误差会在较少的迭代中变得可见。净结果是，我们可以正确计算数字的n值随着精度越少而变得更低。如果您具有D个十六进制数字的精度，则为D * 4位。每次迭代都会引入0.5位的误差到最不重要的位中，因此进行2次迭代时，LSB可能是错误的。在求和期间，这些误差被添加，因此被累积。如果加总的错误数量达到了最重要的数字中的LSB，则提取的单个数字将是错误的。粗略地说，那就是N> 2 **（D-0.75）时。（根据某些对数基数的正确性。）通过经验推断您的数据，似乎近似拟合为N〜（2 **（2.05 * D）），尽管数据点很少，因此可能无法准确预测。您选择的BBP算法是迭代的，因此计算序列中的数字将需要越来越长的时间。要计算0..n位数字，需要O（n ^ 2）步骤。维基百科文章给出了一个公式，用于计算第n个数字，它不需要迭代，只需要指数和有理数。这不会遭受与迭代算法相同的精度损失，您可以按需要以恒定时间（或者最坏情况下对数类型，具体取决于指数与模数的实现）计算pi的任何数字，因此计算n位数字将需要O（n）时间，可能是O（n log n）。