3

这里我将使用符号

在此处输入图像描述

可以通过计算它然后应用定义来找到一个数字的连分数，但这需要至少 O(n) 位的内存才能找到一个0，一个1 ...一个n，实际上它是一个多更差。使用双浮点精度只能找到0、1 ... a 19。

另一种方法是使用这样一个事实：如果 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 有理库将 x、y 和 z 表示为绝对精度，我只能使用 2 1/3 的双精度准确地获得 floor(x + y*2 1/3 + z *2 2/3 )和2 2/3因为我只需要它在真实值的 1/2 以内。不幸的是，x、y 和 z 的分子和分母增长得相当快，如果你使用常规浮点数，错误会很快堆积起来。

通过这种方式，我能够在一小时内计算出0、1 ... 10000，但不知何故，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;
}
4

2 回答 2

2

拉格朗日算法

例如，该算法在 Knuth 的著作 The Art of Computer Programming, vol 2 中进行了描述（Ex 13 in section 4.5.3 Analysis of Euclid's Algorithm, p. 375 in 3rd edition）。

让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

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

警告

在函数内部计算的多项式的系数很快就变成了相当大的整数。所以这种方法需要一些其他语言的 bigint 实现（纯 python3 处理它，但例如 numpy 没有。）

于 2016-12-17T22:39:51.903 回答
1

您可能会更幸运地计算 2^(1/3) 以达到高精度，然后尝试从中推导出连分数，使用区间算术来确定精度是否足够。

这是我在 Python 中的尝试，使用 Halley 迭代计算 2^(1/3) 的定点。死代码是通过牛顿迭代比 Python 更有效地计算定点倒数的尝试——没有骰子。

我的机器上的时间大约是 30 秒，大部分时间都在尝试从定点表示中提取连分数。

prec = 40000
a = 1 << (3 * prec + 1)
two_a = a << 1
x = 5 << (prec - 2)
while True:
    x_cubed = x * x * x
    two_x_cubed = x_cubed << 1
    x_prime = x * (x_cubed + two_a) // (two_x_cubed + a)
    if -1 <= x_prime - x <= 1: break
    x = x_prime

cf = []
four_to_the_prec = 1 << (2 * prec)
for i in range(10000):
    q = x >> prec
    r = x - (q << prec)
    cf.append(q)
    if True:
        x = four_to_the_prec // r
    else:
        x = 1 << (2 * prec - r.bit_length())
        while True:
            delta_x = (x * ((four_to_the_prec - r * x) >> prec)) >> prec
            if not delta_x: break
            x += delta_x
print(cf)
于 2016-12-17T19:18:27.900 回答