Suche nach der Fortsetzung Anteil von 2^(1/3) bis sehr hohe Präzision

Question

Hier werde ich die Schreibweise

Es ist möglich, die Verwendung fortgesetzt Anteil einer Zahl zu finden, indem sie dann die Berechnung Anwenden der Definition, aber das erfordert mindestens O (n) Bits des Gedächtnisses, um eine , eine n zu finden, in der Praxis ist es ein viel schlechteres. Unter Verwendung der doppelten Gleitkommapräzision kann nur ein , ein ... ein gefunden werden. Eine Alternative besteht darin, die Tatsache zu verwenden, dass, wenn a, b, c rationale Zahlen sind, eindeutige Rationalitäten p, q, r existieren, so dass 1/(a ​​+ b * 2 ^1/3 + c * 2 2/3) = x + y * 2 1/3 + z * 2 ^2/3, nämlich

wenn ich also repräsentieren x, y und z auf absolute Präzision mit dem Boost rationale Lib kann ich Stock bekommen (x + y * 2 ^1/3 + z * 2 ^2/3) genau nur mit doppelter Genauigkeit für 2 ^1/3 und 2 ^2/3 weil ich brauche es nur innerhalb 1/2 des wahren Wertes zu sein. Leider wachsen die Zähler und Nenner von x, y und z erheblich schneller, und wenn Sie stattdessen normale Gleitkommazahlen verwenden, häufen sich die Fehler schnell an.

Auf diese Weise konnte ich eine , ein berechnen ... ein in weniger als einer Stunde, aber irgendwie kann Mathematica tun, dass in 2 Sekunden. Hier ist mein Code als Referenz

#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

Sie werden vielleicht mehr Glück haben Rechen 2^(1/3) mit hohen Genauigkeit und dann zu versuchen, die weiteren Fraktion aus, dass abzuleiten, Intervall-Arithmetik, um festzustellen, ob die Genauigkeit ausreichend ist.

Hier ist mein Stich in Python, mit Halley-Iteration 2^(1/3) in Fixpunkt zu berechnen. Der tote Code ist ein Versuch, Festkomma-Kehrwerte effizienter als Python mittels Newton-Iteration zu berechnen - kein Würfel.

Timing von meiner Maschine ist etwa dreißig Sekunden, verbrachte meist versucht, die Fortsetzung Fraktion aus der Festkommadarstellung zu extrahieren.

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

Der Lagrange Algorithmus

Der Algorithmus wird zum Beispiel in Knuth Buch The Art of Computer Programming, Band 2 (Ex 13 in Abschnitt 4.5.3 Analyse von Euklids Algorithmus, p. 375 in 3. Auflage).

Lassen Sie f ein Polynom ganzzahliger Koeffizienten sein, deren einzige reelle Wurzel eine irrationale Zahl x0 > 1 ist. Dann berechnet der Lagrange-Algorithmus die fortlaufenden Quotienten des fortlaufenden Bruchteils von x0.

Ich setzte es in 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 

Es ca. 1 s auf meinem Computer cf([-2, 0, 0, 1], 10000) zu laufen dauert. (Die Koeffizienten entsprechen dem Polynom x^3 - 2, dessen einzige reelle Wurzel 2^(1/3) ist.) Die Ausgabe stimmt mit der von Wolfram Alpha überein.

Caveat

Die Koeffizienten der Polynome innerhalb der Funktion schnell sehr groß werden ganze Zahlen ausgewertet. Dieser Ansatz benötigt also eine Bigint-Implementierung in anderen Sprachen (Pure python3 behandelt das, aber zum Beispiel numpy nicht.)