Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
# 無限降下法 (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 ""
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.