A Perfect Night For Bananafish

だがそれを語るには人生は短すぎる

円周率を求める

Mathクラスで定義された定数PIとついでにアークタンジェントで計算したπ。
module Math

include Math
p PI #=>3.14159265358979
p atan2(1, 1)*4 #=>3.14159265358979

円周率を求める

とりあえず小数点以下8桁まで

require 'open-uri'
open("http://www.google.com/search?q=pi") { |f|
	while line=f.gets
		if /pi =(.*?)</=~line
			p $1.to_f
			f.close
			exit
		end
	end
} 

反省はしていない。

円周率を求める

もう少しまじめに計算してみる。
wikipedia:円周率にあるarcsin(x)のテイラー展開にx=1/2を代入した
 \frac{\pi}{6}=\sum_{n=0}^\infty\frac{(2n)!}{(2^{4n+1})(n!)^2{(2n+1)}

階乗を求める - A Perfect Night For Bananafishで作った階乗を求めるfact(n)を使って、
pi.rb

#階乗
def fact(n)
	f=1
	if (n==0)
		return f
	end
	(1..n).each{ |i|
		f=f*i
	}
	return f
end

def pi(n)
	t=0
	0.upto(n) { |i|
		t +=fact(2*i).to_f/((2**(4*i+1))*(fact(i)**2)*(2*i + 1))
	}
	return t*6
end

(1..30).each { |i|
	puts "#{i}:	#{pi(i)}"
}

整数同士の除算にならないように分子に.to_fをつけて浮動小数に変えている。
n=30ぐらいまで計算してみると

>ruby pi.rb
1:	3.125
2:	3.1390625
3:	3.14115513392857
4:	3.14151117234003
5:	3.14157671577487
6:	3.14158942531912
7:	3.14159198235838
8:	3.14159251115786
9:	3.14159262287062
10:	3.14159264687556
11:	3.14159265210589
12:	3.14159265325874
13:	3.14159265351534
14:	3.14159265357293
15:	3.14159265358595
16:	3.14159265358891
17:	3.14159265358959
18:	3.14159265358975
19:	3.14159265358978
20:	3.14159265358979
21:	3.14159265358979
22:	3.14159265358979
23:	3.14159265358979
24:	3.14159265358979
25:	3.14159265358979
26:	3.14159265358979
27:	3.14159265358979
28:	3.14159265358979
29:	3.14159265358979
30:	3.14159265358979

n=20ぐらいで落ち着くようだ。
nが75以上に増やすとwarning: Bignum out of Float rangeに、n=86以上でNaN(Not a number)になった。
少数以下の桁数を稼ぐためのものもそのうちやってみよう。

モンテカルロ法で円周率を求める

長さ1の正方形に囲まれた扇形(四分の一円)を書き、ランダムに点を打った点が扇形の中にある場合の数の割合を4倍すればπが求まるというもの。
詳しい説明はgoogle:モンテカルロ法 円周率で。
montecarlo_pi.rb

include Math
def pi(n)
	r=0
	n.times{
		if sqrt(rand()**2+rand()**2)<1
			r+=1
		end
	}
	return 4.0*r/n
end

N=100
Trial=10000
t0=Time.now
p=0
(1..N).each {|i|
	p+=pi(Trial)
}
puts p/N
puts Time.now-t0

10000回点を打って求めたπを100回実行して平均をとったものおよびかかった時間を出力する。
というのを5回実行してみた

>ruby montecarlo_pi.rb
3.141872
3.109
>ruby montecarlo_pi.rb
3.138012
3.105
>ruby montecarlo_pi.rb
3.144284
3.14
>ruby montecarlo_pi.rb
3.140536
3.117
>ruby montecarlo_pi.rb
3.140592
3.136

出力の上が円周率で下がかかった時間。かかった時間が円周率っぽいのは奇跡。

ガウス=ルジャンドルのアルゴリズムで円周率を求める

wikipedia:ガウス=ルジャンドルのアルゴリズムをそのまま。
gl_pi.rb

include Math
def gl_pi(n)
	a=1
	b=sqrt(1/2.0)
	t=1/4.0
	p=1
	(1..n).each{
		x=(a+b)/2
		y=sqrt(a*b)
		t=t-p*(a-x)**2
		a=x
		b=y
		p=2*p
	}
	return ((a+b)**2)/(4*t)
end
(1..5).each { |i|
	puts gl_pi(i)
}

とりあえず5回まで実行。収束がはやい。

>ruby gl_pi.rb
3.14057925052217
3.14159264621354
3.14159265358979
3.14159265358979
3.14159265358979