实现算法来计算R中的R
我试图在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,有没有在我的代码中的错误或者是它的算法?
你的代码只是因为它不符合wiki页面中写的算法而变得混乱起来。正确的版本是这样的:
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
好的......我不会想到检查代码本身的准确性+1 –
谢谢你的回答。正如我在我的问题的第一句中写的,我试图实现算法的(给定)变体,但由于结果我不确定它的准确性。这就是为什么我在问题结束时询问R代码是否存在错误,或者问题是算法本身!所以如果变化超过25次迭代是错误的,但变化的实施是可以的,这将回答我的问题,并对我非常有帮助。那么你是否在实现“错误”算法的过程中发现了任何问题? – amadeusamadeus
@amadeusamadeus如果你想调试算法,我认为math.stackexchange.com的人应该能够帮助你。我不是100%确定你是如何派生出来的......它*看起来像Householder的一些反正切值的算法,它给出了$ \ pi $的一小部分。无论如何,你的算法匹配精美格式的公式,所以没有问题。最后要考虑的一件事:大多数算法不是基于像您所做的那样运行固定迭代而终止,而是当该值达到目标精度时终止。 50次迭代可能会达到浮点精度的限制... – AdamO
哪里'i'从何而来? – AdamO
@AdamO它应该是'iteration',而不是'i' –
在玩了大约20分钟的代码之后,我没有确切的答案,而且我不确定是否有解释可能是由于浮点运算的局限性。在我开始看到负数之前,我已经完成了将近50次的迭代。我认为经过这么多迭代之后,'z'中的'2 ^迭代'项变得如此之大,并且'x^2 - y^2'项变得如此之小,以至于四舍五入等开始起作用。你看到的负数只是一个神器。 –