ピタゴラス数は超高速で求まる

ピタゴラス数の高速な求め方について、あれこれ考えた。 ピタゴラス数とは、x*x+y*y==z*z を満たす整数の３つ組 x,y,z のすべての組み合わせを作ってからやると、すごい量になる。 せめて、xとyとかだけにすれば良いか？

などということは全然考えずに、整数 m, n (m>n>0) に対して、 (m*m-n*n, 2*m*n, m*m+n*n) の３つ組がピタゴラス数になるというのを使った。

もちろん、互いに素を除かないとダメだが、それはプログラム中でちょこっと対応。

これですべてのピタゴラス数を見つけることができる証明は省略！ ネットで調べるなり、整数論の本を見るなりして探しましょう。

ということで、新しいプログラム。

triads         :: Int -> [(Int,Int,Int)]
triads n       =  [(x,y,z) | (x,y,z) <- map toabc (doubles (findlimit n)), z<=n, gcd x y == 1]

findlimit    :: Int -> Int
findlimit n  = fromIntegral (truncate (sqrt (fromIntegral n)))

doubles      :: Int -> [(Int,Int)]
doubles n    =  [(x,y) | x <- [1..n-1], y <- [x+1..n]]

toabc        :: (Int,Int) -> (Int,Int,Int)
toabc (n,m)  = (m*m-n*n,2*m*n,m*m+n*n)

これを実行すると、どうなるだろうか。

*Main> triads 20
[(3,4,5),(15,8,17),(5,12,13)]

ということで、20以下のピタゴラス数が全部求まっている。 出てくる順番は、m,n に依存した順番なので、分かり辛いかも。 ピタゴラス数としての順番にしたい人は、自分で変更してください。

さて、どこまでできるだろうか。 全部表示されても困るので、ピタゴラス数の個数だけを求めてみよう。

*Main> (length.triads) 100
16
*Main> (length.triads) 1000
158
*Main> (length.triads) 10000
1593
*Main> (length.triads) 100000
15919
*Main> (length.triads) 1000000
159139
*Main>

こんな感じで、非力な32ビットノートパソコンでも、簡単に求まるのだ。

1000000(100万)までの最後のピタゴラス数は

*Main> (last.triads) 1000000
(1413,998284,998285)

となる。y+1==z の場合になっている。まあ、これは偶然だが。

ということで、Haskellは非常に高速！ というより、計算の仕方が良かっただけかな。

数学的にうまいやり方が見つかっている場合には、使いましょう。 知らないと馬鹿なログラムを書いてCPUを無駄遣いします。 １万倍、あるいはそれ以上の差が出ることはよくあります。