3

ここでは、表記法を使用します

ここに画像の説明を入力

数値を計算してから定義を適用することで、数値の連分数を見つけることができますが、 0 、 a 1 ... a n を見つけるには少なくとも O(n) ビットのメモリが必要です。悪い。倍精度浮動小数点を使用すると、0、1 ... 19しか検出できません。

別の方法として、a,b,c が有理数である場合、1/(a+b*2 1/3 +c*2 2/3 ) = x +y*2 1/3 +z*2 2/3、つまり

ここに画像の説明を入力

したがって、ブースト有理ライブラリを使用して 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を 1 時間以内に計算できましたが、どういうわけか 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 (セクション 4.5.3 Analysis of Euclid's Algorithm、p. 375、第 3 版の例 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

私のコンピュータでは を実行するのに約 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) を高い精度で計算してから、区間演算を使用して精度が十分かどうかを判断し、そこから連分数を導出しようとする運が良いかもしれません。

ハレー反復を使用して、固定小数点で 2^(1/3) を計算します。デッド コードは、固定小数点の逆数を Python よりも効率的に Newton 反復を使用して計算しようとする試みです。ダイスはありません。

私のマシンからのタイミングは約 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 に答える