jacobsthal.rb

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 のときは
						# 平方剰余の第２補充則より計算する
						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"