7

我正在寻找在低内存环境中计算 Pi 的第n个数字。由于我没有可用的小数，Python 中的这种仅整数的 BBP 算法是一个很好的起点。我一次只需要计算一位 Pi。如何确定我可以设置的最低 D，“工作精度的位数”？

D=4 给了我很多正确的数字，但有几个数字会相差一位。例如，计算精度为 4 的数字 393 给我 0xafda，我从中提取数字 0xa。但是，正确的数字是 0xb。

无论我将 D 设置得多高，似乎测试足够多的数字都会发现公式返回的值不正确。

当数字与另一个数字“接近”时，我尝试提高精度，例如 0x3fff 或 0x1000，但找不到任何好的“接近”定义；例如，在数字 9798 处计算得到 0x c de6 ，它与 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]) )
4

2 回答 2

3

无论我将 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)。

于 2010-05-29T13:33:55.173 回答
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)
于 2021-08-04T12:12:19.803 回答