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
、私のコードにバグがあるか、それはアルゴリズムのですか?
どこから来る 'i'ん与えますか? – AdamO
@AdamOそれは 'iteration'ではなく、' i'であると仮定されています –
あなたのコードを約20分間演奏した後、正確な答えはなく、説明以外おそらく浮動小数点演算の制限に起因します。私は負の数を見始める前に、ほぼ50回反復作業をするようにしました。何回も繰り返すと、 'z'の' 2^iteration'項が非常に大きくなり、 'x^2 - y^2'項が小さくなり、丸めなどが始まります。あなたが見る負の数は単なる人工物です。 –