R에서 파이를 계산하는 알고리즘 구현

Question

Brent-Salamin algorithm의 변형을 R에 구현하려고합니다. 처음 25 번의 반복에서는 잘 작동하지만 예기치 않게 작동하여 부정적인 결과를 반환합니다.

제가 구현하려는 알고리즘은 다음과 같습니다

n은 현재 반복이다

initial values: 
x_0 = 1; y_0 = 1/sqrt(2); z_0 = 1/2 

x_n = (x_n-1 + y_n-1)/2 
y_n = sqrt(x_n-1 * y_n-1) 
z_n = z_n-1 - 2^n * (x_n^2-y_n^2) 

p_n = (2 * x_n^2)/z_n 

.

더 아름답게 서식이 지정된 수식은 here입니다.

내가 알아 낸 코드는 다음과 같습니다

mypi <- function(n){ 

    x = 1 
    y = 1/sqrt(2) 
    z = 1/2 
    iteration = 0 

    while(iteration < n){ 
    iteration = iteration + 1 

    newx = (x + y)/2 
    y = sqrt(x * y) 
    x = newx 
    z = z-(2^iteration * (x^2 - y^2)) 
    p = (2 * x^2)/z 
    } 

    return(p) 
} 

출력 : 내가 R에 새로운 오전 그래서으로

> mypi(10) 
[1] 3.141593 
> mypi(20) 
[1] 3.141593 
> mypi(50) 
[1] -33.34323 

, 내 코드에 버그가 있거나 알고리즘입니까?

Answer 1

위키 페이지에 작성된 알고리즘과 일치하지 않기 때문에 코드가 단순히 망가집니다. 원래 올바른 버전은 다음과 같습니다

mypi <- function(n){ 

    x = 1 
    y = 1/sqrt(2) 
    z = 1/4 
    p <- 1 

    iteration = 0 

    while(iteration < n){ 
    iteration = iteration + 1 

    newx = (x + y)/2 
    y = sqrt(x * y) 
    # x = newx 
    # z = z-(2^iteration * (x^2 - y^2)) 
    z = z- p* (x-newx)^2 
    p = 2*p 
    x = newx 
    } 

    (newx + y)^2/(4*z) 
} 

는

> mypi(10) 
[1] 3.141593 
> mypi(20) 
[1] 3.141593 
> mypi(50) 
[1] 3.141593