http://projecteuler.net/index.php?section=problems&id=143

やっと解けた。

import Control.Monad
import Data.List
import qualified Data.Map as M
import qualified Data.IntSet as S
eisenstein :: Int -> [(Int,Int)]
eisenstein upp = do n <- [1..sqrt' upp -1]
m <- [n+1..min (div (upp-n*n) (2*n)).sqrt'$ upp + n^2]
guard$ mod (m-n) 3 /= 
guard$ gcd n m == 1
let a = (m+n)*(m-n)
b = n*(2*m+n)
k <- [1..div upp $ a+b]
return (k*a,k*b)
where sqrt' = floor.sqrt.fromIntegral
torricelli :: M.Map Int S.IntSet -> Int -> [Int]
torricelli eisen p = case M.lookup p eisen of
Nothing -> []
Just qs -> do q <- S.toList qs
r <- findR (S.toList qs) q
return$ p+q+r
where findR qs q = case M.lookup q eisen of
Nothing -> []
Just rs -> S.toList. S.intersection rs . S.fromList $ dropWhile(<=q) qs
p143 n = nub.takeWhile(<n).sort.concatMap (torricelli e)$ [1..div n 3]
where e = foldl insert' M.empty .eisenstein$  n
insert' mp (a,b) = let [p,q] = sort [a,b]
in M.insertWith f p (S.singleton q) mp
f new old = let [n] = S.toList new
in S.insert n old
main = print.sum.p143$11*10^4

本質的な部分は

p^2+q^2+pqが平方数になるような(p,q)を探す

であると思う。

これ、ピタゴラス数に似ているんですよ。

実は、ピタゴラス数のような生成式が存在すると。

それを利用して解いた。

brute forceとは比べ物にならないくらい速くなった。

数学ってすごいね、って思った。