Project Euler Problem 094

import Data.List
import Data.Maybe
import Data.Ratio
-- 全既約ピタゴラス数生成
pita1 :: (Num t, Ord t) => [t] -> [(t, t, t)]
pita1 brange= [(m*m-n*n,2*m*n,m*m+n*n)|m<-brange,n<-brange,n t -> [(t, t, t)]
pita2 dst   = [(m*m-n*n,2*m*n,m*m+n*n)| n<-[1..dst], m<-[(n+1),(n+3)..dst], gcd m n == 1]
-- 生成行列を使った方法。これも既約なものだけを生成し定数倍は生成しない
pita3 :: [(Integer, Integer, Integer)]
pita3 = imp (3,4,5)
  where
    imp (a,b,c) = (a,b,c):
                  concat [imp (a',b',c')|
                          (a', b', c') <- [(a-2*b+2*c, 2*a-b+2*c, 2*a-2*b+3*c),
                                           (a+2*b+2*c, 2*a+b+2*c, 2*a+2*b+3*c),
                                           (-a+2*b+2*c, -2*a+b+2*c, -2*a+2*b+3*c)]]

-- この既約バージョンを使って、上限が分かってるならこんな感じで整数倍を試せる
pitatimes :: (Enum t, Num t) => [t] -> t -> [[t]]
pitatimes l ulimit= [[s,2*s..ulimit]| s<-l]
-- 3辺の和が問題の場合はm*m-n*n+2*m*n+m*m+n*n=2*m*(m+n)を直接利用しても良い

n :: Integer
n=10^9
un :: Integer
un=(n+2)`quot`6
isSquare :: Integral a => a -> Bool
isSquare n = (round . sqrt $ fromIntegral n) ^ 2 == n
sum1 :: Integer
sum1 = sum $ map (\x -> 3*x+1) $ filter (\a -> isSquare((3*a+1)*(a-1))) [2..((n-2)`quot`3)]
sum2 :: Integer
sum2 = sum $ map (\x -> 3*x-1) $ filter (\a -> isSquare((3*a-1)*(a+1))) [2..((n+2)`quot`3)]
ans1 :: Integer
ans1 = sum1 + sum2

-- 518408346
-- コンパイルして6minぐらいかかる。boardにはもっと高速な解が乗っている
main :: IO ()
main = print ans1