Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
require 'prime'
# 平方剰余の相互法則を用いて「p が法 q の平方剰余 (quadratic residue) かどうか」を判定する
#
# input:
# p, q: odd prime
# 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
# mod p のヤコブスタール和 S(a) を計算する
def jacobsthal a, p
sum = 0
p.times do |x|
xx = (x*(x*x + a)) % p
symbol = 0
if xx != 0
# xx は p と互いに素
symbol = 1
# ルジャンドル記号を計算するため xx を素因数分解
factors = Prime.prime_division(xx)
factors.each do |factor|
q = factor[0]
e = factor[1]
# 指数が偶数のときにルジャンドル記号を計算する
if q >= 2 && e % 2 == 1
if q != 2
# q が奇素数のときは
# 平方剰余の相互法則より計算する
symbol *= qr(q, p)
else
# q = 2 のときは
# 平方剰余の第2補充則より計算する
if p % 8 == 1 || p % 8 == 7
symbol *= +1
else
symbol *= -1
end
end
end
end
end
# ルジャンドル記号を足し合わせる
sum += symbol
end
return sum
end
x = jacobsthal( 1, 2017 )
y = jacobsthal( 5, 2017 )
puts "S(1) = #{x}"
puts "S(5) = #{y}"
puts "2017 = (#{x}/2)^2 + (#{y}/2)^2"
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.