fermat_descent2.rb

require 'prime'


# 無限降下法 (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



# 出力結果の見た目を整形
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



# 繰り返し二乗法により
# a^e (mod p) を計算する
def power_mod a, e, p
	sum = 1
	x = a

	# e を二進展開しつつべき乗していく
	while e > 0 
		if e % 2 == 1
			sum = (sum * x) % p
		end

		e = e / 2
		x = (x*x) % p
	end

	return sum
end



# 法 p における (-1) の平方根を探索する
# 
# input: 
#     p: prime 
# output: 
#     x: integer  s.t. x^2 == -1 (mod p)
def sqrt_negative_one_modulo p
	# 法 p で平方非剰余になる奇素数 b を探す
	b = 3
	while b < p
		if Prime.prime?(b) && qr(b, p) == -1
			break
		end
		b += 2	
	end

	# オイラーの基準により, b^{(p-1)/4} が求める平方根
	return power_mod(b, (p-1)/4, p)

end



# 平方剰余の相互法則を用いて「p が法 q の平方剰余 (quadratic residue) かどうか」を判定する
#
# input:  
#    p, q: odd primes
# output: 
#    p が q の 平方剰余 => +1, 平方非剰余 => -1
def qr p, q
	sgn = 1
	if p % 4 == 3 && q % 4 == 3
		sgn = -1
	end

	# 平方剰余の相互法則を使って p と q を入れ替えて
	# 「q が法 p で平方剰余かどうか」を判定する
	a = q % p  # a == q (mod p)

	p.times do |x|
		if (x*x - a) % p == 0
			return sgn
		end
	end 

	return -sgn
end




p = 2017
x = sqrt_negative_one_modulo(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 = sqrt_negative_one_modulo(p)
y = 1
k = (x**2 + y**2) / p

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