정밀 매우 높은 여기

Question

나는 표기법을 사용합니다 2^(1/3)의 지속적인 부분을 찾기

는 다음을 계산하여 숫자의 지속적인 부분을 찾을 수있다 정의를 적용하지만, 적어도 O (n) 비트의 메모리가 필요하며, , ... a _n 실제로는 훨씬 더 나쁩니다. 이중 부동 소수점 정밀도를 사용하면 , ... 을 찾을 수만 있습니다.

a, b, c가 유리수 인 경우 1/(a ​​+ b * 2 ^1/3 + c * 2)와 같은 고유 한 합계 p, q, r이 존재한다는 사실을 사용할 수도 있습니다. 2/3) = X + Y * 2 1/3 + Z * 2 즉 2/3 그래서

I 절대 정밀도로 X, Y 및 Z를 나타내면 부스트 합리적인 lib를 사용하여 바닥을 얻을 수 있습니다 (x + y * 2 ^1/3 + z * 2 ^2/3) 정확히 2의 정확도를 사용하는 경우에만 ^1/3 및 2 ^2/3입니다. 왜냐하면 저는이 값이 실제 값의 1/2 이내 여야하기 때문입니다. 불행하게도 x, y, z의 분자와 분모는 상당히 빠르게 증가합니다. 보통의 수레를 사용하면 오류가 빨리 쌓입니다.

나는 을 계산 할 수 있었다 이런 식으로하는 ... 시간 아래에서 하지만, 어떻게 든 티카는 것을 할 수 이초한다. 여기에 당신은 높은 정확도로^(1/3) (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

내 코드입니다.

여기는 파이썬에서 Halley 반복을 사용하여 고정 점에서 2^(1/3)을 계산합니다. 데드 코드는 주사위가없는 Newton 반복을 통해 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) 

Answer 2

라그랑 알고리즘 알고리즘 누스 책 컴퓨터 프로그래밍, 4.5.3 분석 유클리드 알고리즘 섹션 2 집 (예 13 (P)의 기술. 375에서, 예를 들어 설명한다

제 3 판).

f은 실수 계수 만의 실수가 x0 > 1 인 정수 계수의 다항식이라고합시다. 그런 다음 Lagrange 알고리즘은 연속 분수 인 x0의 연속적인 지수를 계산합니다.

나는 cf([-2, 0, 0, 1], 10000)를 실행하는 컴퓨터에 파이썬

그것은 1 초 정도 걸립니다

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 

에서 그것을 구현했습니다. (계수는 다항식 x^3 - 2에 해당하며 실제 루트는 2^(1/3)입니다. 출력은 Wolfram Alpha의 결과와 일치합니다.

주의 할

빠르게 매우 큰 정수가 함수 내에서 평가 된 다항식의 계수. 그래서이 접근법은 다른 언어에서 bigint 구현을 필요로합니다. (순수 파이썬 3이 그것을 다루지 만, 예를 들어 numpy는 그렇지 않습니다.)