http://projecteuler.net/index.php?section=problems&id=86
l(m2 - n2) か 2lmn のどちらかがM。
primes = 2:3:5:7:[ p + 6 | p <- tail (tail primes) ] calc_exp n p | mod n p /= 0 = (n, 0) | otherwise = (m, e + 1) where (m, e) = calc_exp (div n p) p factorize :: Int -> [Int] -> [(Int,Int)] factorize 1 ps = [] factorize n (p:ps) | e > 0 = (p,e):(factorize m ps) | otherwise = factorize m ps where (m, e) = calc_exp n p val f = foldr (\(p,e) y -> p^e * y) 1 f divisors [] = [([], [])] divisors ((p,e):fs) = [ (neg ((p,e1):a),neg ((p,e-e1):b)) | e1 <- [0..e], (a,b) <- divisors fs ] neg [] = [] neg ((p,e):a) = if e == 0 then neg a else (p,e):(neg a) divs n = [ (val f1, val f2, val f3) | (f1, f2, f3) <- g f ] where f = factorize n primes g f = [ (f1,f3,f4) | (f1,f2) <- divisors f, (f3,f4) <- divisors f2 ] ds p | odd p = sum (map num1 a) | otherwise = (sum (map num1 a)) + (sum (map num2 b)) where a = divs p b = divs (div p 2) num1 (l,s,t) | s > t && odd s && odd t && gcd s t == 1 = f a b | otherwise = 0 where a = l * s * t b = l * ((s * s - t * t) `div` 2) num2 (l,m,n) | m > n && odd (m + n) && gcd m n == 1 = f a b | otherwise = 0 where a = 2 * l * m * n b = l * (m * m - n * n) f a b | a < b && b < 2 * a = a - (div (b - 1) 2) | b < a = div b 2 | otherwise = 0 limit = 10^6 a = map ds [1..] b = 0:[ m + n | (m,n) <- zip a b ] main = print (fst (head (filter (\(_,m) -> m > limit) (zip [0..] b))))
