読者です 読者をやめる 読者になる 読者になる

猫尾製作所

あまりアテにしないでね

sqrt(d) の連分数展開

数学 プログラミング 連分数 代数学

コード

# sqrt(d) の連分数展開

# 最大公約数(2つの整数)
def gcd2(a, b)
	if a == 0 || b == 0 then
		return 0
	end
	
	while a % b != 0 do
		c = a % b
		a = b
		b = c
	end
	
	return b
end

# 最大公約数(3つの整数)
def gcd3(a, b, c)
	return gcd2(gcd2(a, b), c)
end

# 引数が与えられていなければ異常終了
if ARGV.length != 1 then
	puts "Please give n as integer."
	exit -1
else
	n  = ARGV[0].to_i
end

# 配列を用意
a = Array.new()
b = Array.new()
c = Array.new()
d = Array.new()

# 漸化式の初期値
a[0] = (Math.sqrt(n)).floor
b[0] = 1
c[0] = -a[0]
d[0] = 1

i = 0

# とりあえず10000回以内で打ち切ることを前提とする
while i < 10000 do
	puts "[#{i}] : #{a[i]}, #{b[i]}, #{c[i]}, #{d[i]}\n"
	
	# p, q, r, s を求める
	s = b[i]*b[i]*n - c[i]*c[i]
	q = b[i]*d[i]
	p = ((b[i]*d[i]*Math.sqrt(n) - c[i]*d[i]) / s).floor
	r = - c[i]*d[i] - p*s
	g = gcd3(q, r, s) # 最大公約数
	g = g == 0 ? 1 : g
	
	i += 1
	d[i] = s / g
	b[i] = q / g
	c[i] = r / g
	a[i] = p
	
	for j in 0..(i-1) do
		if a[i] == a[j] && b[i] == b[j] && c[i] == c[j] && d[i] == d[j] then
			puts "Circulation : [#{j}]..[#{i-1}]"
			exit 0
		end
	end
end

実行例

% ruby sqrtd.rb 7
[0] : 2, 1, -2, 1
[1] : 1, 1, -1, 3
[2] : 1, 1, -1, 2
[3] : 1, 1, -2, 3
[4] : 4, 1, -2, 1
Circulation : [1]..[4]

% ruby sqrtd.rb 77
[0] : 8, 1, -8, 1
[1] : 1, 1, -5, 13
[2] : 3, 1, -7, 4
[3] : 2, 1, -7, 7
[4] : 3, 1, -5, 4
[5] : 1, 1, -8, 13
[6] : 16, 1, -8, 1
Circulation : [1]..[6]

% ruby sqrtd.rb 177
[0] : 13, 1, -13, 1
[1] : 3, 1, -11, 8
[2] : 3, 1, -10, 7
[3] : 2, 1, -12, 11
[4] : 8, 1, -12, 3
[5] : 2, 1, -10, 11
[6] : 3, 1, -11, 7
[7] : 3, 1, -13, 8
[8] : 26, 1, -13, 1
Circulation : [1]..[8]

番号のすぐ右に出ている整数が連分数展開の結果です(あとは計算用に用いたダミー変数)。Circulation : [m]..[n] とは m 番目の項から n 番目の項までが循環するという意味合いです。