發現的2 ^（1/3）連分數非常高的精度

Question

在這裏，我將使用符號

有可能通過計算它，然後找到一個數的連分數應用的定義，但至少需要O（n）位的內存找到一個 _{...一個 n，實際上它是一個更糟糕的。使用雙浮點精度，只能找到一個，一個 ... 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，即

所以，如果我表示的x，y和z以絕對精度使用升壓理性LIB我可以得到地板（X + Y * 2 ^1/3 + Z * 2 ^三分之二）準確地只使用雙精度2 ^1/3和2 ^2/3，因爲我只需要它在真值的1/2以內。不幸的是，分子和X，Y和Z的分母大幅增長快，如果你使用普通的浮動，而不是錯誤很快堆積起來。

這樣我就能夠計算，一個 ...在一個小時內一個，但不知何故，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

代碼你可能有更多的運氣計算2 ^（1/3），以高精確度，然後試圖從獲得的部分繼續使用區間運算以確定精度就足夠了。

這裏是用Python我在這個刺，使用哈雷迭代計算2 ^（1/3）的固定點。死代碼是通過牛頓迭代比Python更有效地計算定點倒數的嘗試 - 沒有骰子。從我的機器

時間是約三十秒花大多試圖從固定點表示連分數。

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) 

Answer 2

拉格朗日算法

該算法在Knuth的書計算機程序設計，第2卷（出埃及記13 4.5.3分析歐幾里得算法的部分，p的藝術。375中描述第3版）。

讓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沒有。）