我正在尝试实现 Brent-Salamin algorithm 的变体在 R 中。它在前 25 次迭代中运行良好,但随后出现意外,返回负结果。
我要实现的算法是:
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 的新手,我的代码中有错误还是算法?
最佳答案
您的代码完全是一团糟,因为它与维基页面中所写的算法不一致。正确的版本如下所示:
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
关于r - 实现一种算法来计算 R 中的 pi,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/47994742/