Project Euler Problem 023

import Data.List
import qualified Data.IntSet as S

--IntegerではなくIntを使ってといてみた
limit :: Int
limit = 28123
isabundant :: Int -> Bool
isabundant n = n < (subtract n $ product $ map (fromIntegral . f) $ group $ factor (fromIntegral n)) where
  f [x] = x+1
  f l@(x:_) = ((x^(length l+1))-1) `quot` (x-1) --等比数列の和の公式利用
abundants :: [Int]  
abundants = filter isabundant [1..limit]  
sets :: [Int] -> [Int] -> S.IntSet -> S.IntSet
sets [x] ys s = foldl (\a c -> S.insert (x+c) a) s ys
sets (x:xs) ys s = sets xs ys (foldl (\a c -> S.insert (x+c) a) s ys)
ansets :: S.IntSet
ansets = sets abundants abundants S.empty
-- 4179871
-- コンパイルして8sec
main :: IO ()
main = print $ sum $ [1..limit] \\ S.elems ansets



--以下003からのコピペ。
primes :: [Integer]
primes = 2 : ([3,5..] `minus` join [[p*p, p*p+2*p..] | p <- primes']) 
  where
    primes' = 3 : ([5,7..] `minus` join [[p*p, p*p+2*p..] | p <- primes'])   
    join  ((x:xs):t)    = x : union xs (join (pairs t))
    pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t                      
    union (x:xs) (y:ys) = case (compare x y) of 
      LT -> x : union  xs  (y:ys)
      EQ -> x : union  xs     ys 
      GT -> y : union (x:xs)  ys
    union  xs     []    = xs
    union  []     ys    = ys
    minus (x:xs) (y:ys) = case (compare x y) of 
      LT -> x : minus  xs  (y:ys)
      EQ ->     minus  xs     ys
      GT ->     minus (x:xs)  ys
    minus  xs     _     = xs
--必要に応じてgroupする
factor :: Integer -> [Integer]
factor n = factorimpl n primes where
  factorimpl n pri@(p:xs) =
    if p*p>n then [n]
    else if n`rem` p == 0 then
      p:factorimpl (n `quot` p) pri
         else
           factorimpl n xs
divsum :: Integer -> Integer
divsum n =  (product. map tohi .group . factor) n - n where
    tohi [x] = x+1
    tohi l@(x:xs) = (x^(length l+1)-1) `quot` (x-1)
isPrime :: Integer -> Bool
isPrime 1 = False
isPrime n = imp n primes where
  imp n (x:xs) = if x*x>n then True
                 else if rem n x==0 then False
                      else imp n xs      
{-
別解。コンパイルして37sec、あまり早い解ではない
import Control.Monad hiding (join)
import Data.List hiding (union)
--003より
primes :: [Integer]
primes = 2 : ([3,5..] `minus` join [[p*p, p*p+2*p..] | p <- primes']) 
  where
    primes' = 3 : ([5,7..] `minus` join [[p*p, p*p+2*p..] | p <- primes'])   
    join  ((x:xs):t)    = x : union xs (join (pairs t))
    pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t                      
    union (x:xs) (y:ys) = case (compare x y) of 
      LT -> x : union  xs  (y:ys)
      EQ -> x : union  xs     ys 
      GT -> y : union (x:xs)  ys
    union  xs     []    = xs
    union  []     ys    = ys
    minus (x:xs) (y:ys) = case (compare x y) of 
      LT -> x : minus  xs  (y:ys)
      EQ ->     minus  xs     ys
      GT ->     minus (x:xs)  ys
    minus  xs     _     = xs
--必要に応じてgroupする
factor :: Integer -> [Integer]
factor n = factorimpl n primes where
  factorimpl n pri@(p:xs) =
    if p*p>n then [n]
    else if n`rem` p == 0 then
      p:factorimpl (n `quot` p) pri
         else
           factorimpl n xs
divsum :: Integer -> Integer
divsum n =  (product. map tohi .group . factor) n - n where
    tohi [x] = x+1
    tohi l@(x:xs) = (x^(length l+1)-1) `quot` (x-1)
abundant :: [Integer]
abundant = filter (\x->x-(divsum x)<0) [1..28123]
nonex :: Integral a => [a] -> a -> [a]
nonex lis nn = f lis [1..nn] where
  f [] acc = acc
  f (x:xs) acc = f xs $ filter (\y -> y`rem`x/=0) acc

abundants :: (Num a, Ord a) => [a] -> a -> [a]
abundants [] limit= []  
abundants (x:xs) limit= let rest = abundants xs limit in
  merge ((map (+x) (x:xs))) rest --ここだけなんだけど、これじゃダメな気がするな・・→ダメじゃなかった！
  where
    merge xs [] = xs
    merge [] ys = ys
    merge (x:xs) (y:ys) =
      if xlimit then []
        else x:merge xs (y:ys)
      else if x==y then x:merge xs ys
           else if y>limit then []
                else y:merge (x:xs) ys
abundantlist :: [Integer]
abundantlist = abundants abundant 28123


--4179871
main :: IO ()
main = print $ sum $ nonex abundantlist 28123  
-- -}