gaussian_prime.rb

# encoding: utf-8
require 'prime'



def fermat p, x, y

	k = (x**2 + y**2)/p

	b = (x % k)
	b = ((b - k).abs < b) ? (b - k) : b

	d = (y % k)
	d = ((d - k).abs < d) ? (d - k) : d

	u = (x*b + y*d)/k
	v = (x*d - y*b)/k 

	return [u, v]

end


def q_residue p
	p.times do |x|
		r = (x*x + 1) % p 

		if r == 0 then
			return x
		end
	end
end

=begin
p = 29

x = q_residue(p)
y = 1
k = (x*x + y*y) /p 
puts "#{k}*#{p} = #{x}^2 + #{y}^2"

while k > 1 

	res = fermat(p, x, y)
	x = res[0]
	y = res[1]
	k = (x*x + y*y) /p 
	puts "#{k}*#{p} = #{x}^2 + #{y}^2"

end
=end


def get_xy p

	x = q_residue(p)
	y = 1
	k = (x*x + y*y) /p 
	#puts "#{k}*#{p} = #{x}^2 + #{y}^2"

	while k > 1 

		res = fermat(p, x, y)
		x = res[0]
		y = res[1]
		k = (x*x + y*y) /p 
		#puts "#{k}*#{p} = #{x}^2 + #{y}^2"

	end

	return [x, y]
end


ARRAY_SIZE = 100

arr = Array.new(ARRAY_SIZE+1){ Array.new(ARRAY_SIZE+1, 0) }


#=begin
Prime.each(100000) do |p|
	if p == 2 then
=begin
		p Complex(1, 1)
		p Complex(1, -1)
=end
		arr[1][1] = 1

	elsif (p % 4) == 1 then
		t = get_xy(p)
		x = t[0].abs
		y = t[1].abs
=begin
		p Complex(x, y)
		p Complex(-x, y)
		p Complex(x, -y)
		p Complex(-x, -y)
		p Complex(y, x)
		p Complex(y, x)
		p Complex(y, -x)
		p Complex(y, -x)
=end
		if x <= ARRAY_SIZE && y <= ARRAY_SIZE then
			arr[x][y] = 1
			arr[y][x] = 1
		end

	else
=begin
		p Complex(p, 0)
		p Complex(0, p)
		p Complex(-p, 0)
		p Complex(0, -p)
=end
		if p <= ARRAY_SIZE then
			arr[p][0] = 1
			arr[0][p] = 1
		end

	end
end
#=end


(ARRAY_SIZE*2-1).times do |y|
	yy = y - (ARRAY_SIZE-1)
	(ARRAY_SIZE*2-1).times do |x|
		xx = x - (ARRAY_SIZE-1)
		if arr[xx.abs][yy.abs] == 1 then
			print "■"
		else
			print "　"
		end
	end
	puts ""
end