algorithm - 以非常高的精度找到 2^(1/3) 的连分数

Question

这里我将使用符号

可以通过计算然后应用定义来找到一个数的连分数，但这至少需要 O(n) 位内存才能找到 a₀，a ₁ ... a_n，实际上它更糟。使用 double float 只能找到 a₀、a₁ ... a₁₉。

另一种方法是利用以下事实:如果 a、b、c 是有理数，则存在唯一的有理数 p、q、r 使得 1/(a+b*2^1/3 +c*2^2/3) = x+y*2^1/3+z*2^2/3,即

因此，如果我使用 boost rational 库以绝对精度表示 x、y 和 z，我可以获得 floor(x + y*2^1/3+z*2^2/3) 准确地只对 2^1/3 和 2^2/3 使用 double 因为我只需要它在真值的 1/2 以内.不幸的是，x、y 和 z 的分子和分母增长得相当快，如果您使用常规 float ，错误会迅速堆积。

通过这种方式，我能够在不到一个小时的时间内计算出 a₀、a₁ ... a₁₀₀₀₀，但不知何故 mathematica 可以在 2 秒内完成。这是我的代码供引用

#include <iostream>

#include <boost/multiprecision/cpp_int.hpp>
namespace mp = boost::multiprecision;

int main()
{
    const double t_1 = 1.259921049894873164767210607278228350570251;
    const double t_2 = 1.587401051968199474751705639272308260391493;
    mp::cpp_rational p = 0;
    mp::cpp_rational q = 1;
    mp::cpp_rational r = 0;
    for(unsigned int i = 1; i != 10001; ++i) {
        double p_f = static_cast<double>(p);
        double q_f = static_cast<double>(q);
        double r_f = static_cast<double>(r);
        uint64_t floor = p_f + t_1 * q_f + t_2 * r_f;
        std::cout << floor << ", ";
        p -= floor;
        //std::cout << floor << " " << p << " " << q << " " << r << std::endl;
        mp::cpp_rational den = (p * p * p + 2 * q * q * q +
                                4 * r * r * r - 6 * p * q * r);
        mp::cpp_rational a = (p * p - 2 * q * r) / den;
        mp::cpp_rational b = (2 * r * r - p * q) / den;
        mp::cpp_rational c = (q * q - p * r)     / den;
        p = a;
        q = b;
        r = c;
    }
    return 0;
}

Answer 1

拉格朗日算法

该算法在 Knuth 的书 The Art of Computer Programming，第 2 卷(第 4.5.3 节欧几里得算法分析第 375 页第 3 版中的 Ex 13)中有描述。

设f 是整数系数的多项式，其唯一实根是无理数x0 > 1。然后拉格朗日算法计算x0的连分式的连续商。

我是用python实现的

def cf(a, N=10):
    """
    a : list - coefficients of the polynomial,
        i.e. f(x) = a[0] + a[1]*x + ... + a[n]*x^n
    N : number of quotients to output
    """
    # Degree of the polynomial
    n = len(a) - 1

    # List of consecutive quotients
    ans = []

    def shift_poly():
        """
        Replaces plynomial f(x) with f(x+1) (shifts its graph to the left).
        """
        for k in range(n):
            for j in range(n - 1, k - 1, -1):
                a[j] += a[j+1]

    for _ in range(N):
        quotient = 1
        shift_poly()

        # While the root is >1 shift it left
        while sum(a) < 0:
            quotient += 1
            shift_poly()
        # Otherwise, we have the next quotient
        ans.append(quotient)

        # Replace polynomial f(x) with -x^n * f(1/x)
        a.reverse()
        a = [-x for x in a]

    return ans

在我的电脑上运行 cf([-2, 0, 0, 1], 10000) 大约需要 1 秒。 (系数对应于多项式 x^3 - 2，其唯一实根为 2^(1/3)。)输出与 Wolfram Alpha 的输出一致。

警告

函数内部评估的多项式系数很快变成非常大的整数。所以这种方法需要其他语言的一些 bigint 实现(纯 python3 处理它，但例如 numpy 没有。)