fermat_descent.rb

# 無限降下法 (descent)
# 
# input:  (x, y, k, p)  s.t. kp = x^2 + x^2
# output: (x, y, k, p)  s.t.  p = x^2 + x^2, k = 1
def descent(x, y, k, p)

	if k == 1
		return [x, y, 1, p]

	else
		# 計算の途中過程もコンソールに出力する
		result = [x, y, k, p]
		puts pretty_format(result)

		# x (resp. y) の mod k における絶対最小剰余 x_mod_k (resp. y_mod_k) をとる
		x_mod_k = x % k
		if x_mod_k > (x_mod_k - k).abs
			x_mod_k = x_mod_k - k
		end
		y_mod_k = y % k
		if y_mod_k > (y_mod_k - k).abs
			y_mod_k = y_mod_k - k
		end

		# (x, y, k) -> (next_x, next_y, next_k)
		next_k = (x_mod_k**2 + y_mod_k**2) / k
		next_x = (x * x_mod_k + y * y_mod_k) / k
		next_y = (x * y_mod_k - y * x_mod_k) / k

		# (next_x, next_y, next_k) に対して descent を再帰的に実行する
		return descent(next_x, next_y, next_k, p)
	end

end



# 平方剰余を探索する
# 
# input: 
#     a: integer, 
#     p: prime 
# output: 
#     x: integer  s.t. x^2 == a (mod p)
def qr a, p
	p.times do |x|
		if (x*x - a) % p == 0
			return x
		end
	end
end



# 出力結果の見た目を整形
def pretty_format result
	x = result[0]
	y = result[1]
	k = result[2]
	p = result[3]

	return "   #{k} * #{p} = #{x.abs}^2 + #{y.abs}^2"

end



p = 2017
x = qr(-1, p)
y = 1
k = (x**2 + y**2) / p

puts "p = #{p}"
result = descent(x, y, k, p)
puts pretty_format(result)
puts ""


p = 20171201
x = qr(-1, p)
y = 1
k = (x**2 + y**2) / p

puts "p = #{p}"
result = descent(x, y, k, p)
puts pretty_format(result)
puts ""