Поиск непрерывной дроби 2^(1/3) до очень высокой точности

Question

Здесь я буду использовать обозначение

можно найти цепную дробь числа, вычисляя его тогда применяя определение, но для этого требуется, по меньшей мере, O (n) бит памяти, чтобы найти , a ... a _n, на практике это намного хуже. Используя двойную точность с плавающей запятой, можно найти только , a ... a .

Альтернатива заключается в том, что если a, b, c - рациональные числа, то существуют единственные рациональные p, q, r такие, что 1/(a ​​+ b * 2 ^1/3 + c * 2 ^2/3) = х + у * 2 ^1/3 + Z * 2 ^2/3, а именно

Так что, если я представляю х, у и г до абсолютной точностью используя boost boost lib, я могу получить пол (x + y * 2 ^1/3 + z * 2 ^2/3) точно только с использованием двойной точности для 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 с использованием итерации Halley для вычисления 2^(1/3) в неподвижной точке. Мертвый код - попытка вычислить встречные обратные точки более эффективно, чем Python через итерацию Newton - без кубиков.

Сроки с моей машины составляют около тридцати секунд, в основном тратятся на извлечение оставшейся доли из представления с фиксированной точкой.

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

Алгоритм Лагранжа

алгоритм описан, например, в книге Кнута Искусство программирования, том 2 (Ex 13 в разделе 4.5.3 Анализ алгоритма Евклида, стр. 375 в 3-е издание).

Пусть f является многочленом целочисленных коэффициентов, единственным вещественным корнем которого является иррациональное число x0 > 1. Затем алгоритм Лагранжа вычисляет последовательные отношения цепной доли x0.

Я реализовал его в питон

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.

Caveat

Коэффициенты полиномов оцениваются внутри функции быстро становятся достаточно большими целыми числами. Таким образом, этот подход нуждается в некоторой реализации bigint на других языках (с ним работает Pure python3, но, например, numpy не делает.)