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 ((/=0).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 /= 0]

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$ [(0,([[]],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' 0 _ = [[]]
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