我试图在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'项变得如此之小,以至于四舍五入等开始起作用。你看到的负数只是一个神器。 –