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

とあるブログで『二平方の定理』というものを見て、面白いなぁと思った
のですが、手計算だと面倒だったのでスクリプトで該当素数を一覧表示。

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