Problem 152

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

大きな素因数を持つ数から考えていき、残りが２，３になったら全探索を開始。

import Number (primes, isPrime, merge)
import Data.List (groupBy, sortBy)
import Data.Ratio (Ratio, denominator, (%))
import Data.Ord (comparing)
import Data.Function (on)
lim :: Integer
lim = 80
invSq :: (Real a) => a -> Rational
invSq x = 1/((^2).toRational$x)
sumInvSq :: [Integer] -> Rational
sumInvSq = sum.map invSq
multi :: Integer -> [Integer]
multi p = filter indivisible$[p,2*p..lim]
where qs = takeWhile(<=lim).dropWhile(<=p)$primes
indivisible x = all ((/=).mod x) qs
indivDenom :: (Real a) => Integer -> [a] -> Rational -> [[a]]
indivDenom p s x = [t | t <- subsequences s, let d = denominator.(+x).sum.map invSq$ t, mod d p /= ]
subsequences            :: [a] -> [[a]]
subsequences []         =  [[]]
subsequences (x:xs)     =  subsequences xs ++ map (x:) (subsequences xs)
main :: IO ()
main = let p = last.takeWhile (<=div lim 2)$ primes
in print.length.search$ [(,([[]],p))]
search :: [(Rational, ([[Integer]], Integer))] -> [[Integer]]
search [] = []
search ((x,(s,3)):ts) = [u++v | u<-s, v<-search' (1%2-x) pow23] ++ search ts
search ((x,(s,p)):ts) = search$ gather next ++ ts
where q = head.filter isPrime$ [p-1,p-2..2]
next = [(x+sumInvSq u,(map (u++) s,q)) | u <- indivDenom p (multi p) x]
gather = map tie.groupBy ((==) `on` fst).sortBy (comparing fst)
tie (y:ys) = (fst y, (concatMap (fst.snd) $ y:ys, snd.snd$ y))
search' :: Rational -> [Integer] -> [[Integer]]
search'  _ = [[]]
search' _ [] = []
search' x (y:ys) | x < invSq y = search' x ys
| x > sumInvSq (y:ys) = []
| otherwise = map (y:) (search' (x-invSq y) ys) ++ search' x ys
pow23 :: [Integer]
pow23 = takeWhile (<=lim)$2:3:map (*2) pow23 `merge` map (*3) pow23