Edited at

[memo] 『二平方の定理』の素数調査

More than 5 years have passed since last update.

とあるブログで『二平方の定理』というものを見て、面白いなぁと思った

のですが、手計算だと面倒だったのでスクリプトで該当素数を一覧表示。

ruby two_square_theorem.rb [調査限界値(省略時:257)]


two_square_theorem.rb

# coding: UTF-8

require "prime"

def theorem?
lambda { |n|
Prime.prime?(n, Prime::EratosthenesGenerator.new) && n % 4 == 1
}
end

def quadrable
lambda { |n|
v = Math.sqrt(n)
[ v == v.floor.to_f, v.to_i ]
}
end

def prover
squares = {}
lambda { |n|
1.upto(Math.sqrt(n/2)) { |x|
squares[x] ||= (x ** 2)
res, y = quadrable[ n - squares[x] ]
return [x, y] if res
}
return nil
}
end

def output(args = nil, &blk)
blk && blk[args]
end

if __FILE__ == $0

n = ( ARGV.shift || 257 ).to_i

case n
when theorem?
puts output( prover[n] || output { puts "#{n} is not resolved."; nil } || exit(0) ) { |x, y|
print "[ #{n} = #{x}^2 + #{y}^2 = #{x**2} + #{y**2} ] : "
n == x**2+y**2
} ? "OK" : "NG"
else
output { puts "#{n} is not (4n+1)'s Prime." }
end

puts "\n=== [ List ] ==="
l = Math.log10(n).ceil

Prime.each(n, Prime::EratosthenesGenerator.new).grep(theorem?).grep(prover).to_enum.with_index(1) { |n, i|
puts output(prover[n]) { |x, y|
xx, yy = x**2, y**2
fmt = "No. %#{l}d [ %#{l}d = %#{l-1}d^2 + %#{l-1}d^2 = %#{l}d + %#{l}d ] : "
print fmt % [ i, n, x, y, xx, yy ]
n == xx + yy
} ? "OK" : "NG"
}

end



実行例

ruby theorem.rb

[ 257 = 1^2 + 16^2 = 1 + 256 ] : OK

=== [ List ] ===
No. 1 [ 5 = 1^2 + 2^2 = 1 + 4 ] : OK
No. 2 [ 13 = 2^2 + 3^2 = 4 + 9 ] : OK
No. 3 [ 17 = 1^2 + 4^2 = 1 + 16 ] : OK
No. 4 [ 29 = 2^2 + 5^2 = 4 + 25 ] : OK
No. 5 [ 37 = 1^2 + 6^2 = 1 + 36 ] : OK
No. 6 [ 41 = 4^2 + 5^2 = 16 + 25 ] : OK
No. 7 [ 53 = 2^2 + 7^2 = 4 + 49 ] : OK
No. 8 [ 61 = 5^2 + 6^2 = 25 + 36 ] : OK
No. 9 [ 73 = 3^2 + 8^2 = 9 + 64 ] : OK
No. 10 [ 89 = 5^2 + 8^2 = 25 + 64 ] : OK
No. 11 [ 97 = 4^2 + 9^2 = 16 + 81 ] : OK
No. 12 [ 101 = 1^2 + 10^2 = 1 + 100 ] : OK
No. 13 [ 109 = 3^2 + 10^2 = 9 + 100 ] : OK
No. 14 [ 113 = 7^2 + 8^2 = 49 + 64 ] : OK
No. 15 [ 137 = 4^2 + 11^2 = 16 + 121 ] : OK
No. 16 [ 149 = 7^2 + 10^2 = 49 + 100 ] : OK
No. 17 [ 157 = 6^2 + 11^2 = 36 + 121 ] : OK
No. 18 [ 173 = 2^2 + 13^2 = 4 + 169 ] : OK
No. 19 [ 181 = 9^2 + 10^2 = 81 + 100 ] : OK
No. 20 [ 193 = 7^2 + 12^2 = 49 + 144 ] : OK
No. 21 [ 197 = 1^2 + 14^2 = 1 + 196 ] : OK
No. 22 [ 229 = 2^2 + 15^2 = 4 + 225 ] : OK
No. 23 [ 233 = 8^2 + 13^2 = 64 + 169 ] : OK
No. 24 [ 241 = 4^2 + 15^2 = 16 + 225 ] : OK
No. 25 [ 257 = 1^2 + 16^2 = 1 + 256 ] : OK