Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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 ""
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.