我試圖在R中實現Brent-Salamin algorithm的變體。它在前25次迭代中運行良好,但後來出乎意料地返回負結果。實現算法來計算R中的R
算法我想要實現的是:
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
其中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)
}
輸出:
> mypi(10)
[1] 3.141593
> mypi(20)
[1] 3.141593
> mypi(50)
[1] -33.34323
所以,我是新來的R,有沒有在我的代碼中的錯誤或者是它的算法?
哪裏'i'從何而來? – AdamO
@AdamO它應該是'iteration',而不是'i' –
在玩了大約20分鐘的代碼之後,我沒有確切的答案,而且我不確定是否有解釋可能是由於浮點運算的侷限性。在我開始看到負數之前,我已經完成了將近50次的迭代。我認爲經過這麼多迭代之後,'z'中的'2 ^迭代'項變得如此之大,並且'x^2 - y^2'項變得如此之小,以至於四捨五入等開始起作用。你看到的負數只是一個神器。 –