Project Euler Problem 095

import Data.List
factors :: Integer -> [[Integer]]
factors  = factors' [2..]
  where
    factors' (x:xs) n | n < x*x = [[n]]
                      | mod n x == 0 = (map (x:).factors' (x:xs))(div n x) ++ factors' xs n
                      | mod n x /= 0 = factors' xs n
-- from 023                                       
divisors :: Integral a => a -> [a]
divisors x = 1:divs [] 2 x
  where
    divs sofar d x | d*d>x = sofar
                   | d*d == x = d:sofar
                   | x `mod` d == 0 = divs (d:(x `div` d):sofar) (1+d) x
                   | True = divs sofar (1+d) x
dsum :: Integral a => a -> a
dsum x = sum $ divisors x
dsum2 :: Integer -> Integer
dsum2 x = product lis - x
  where
    fact = factorize x
    lis = map (\(a,n) -> (a^(n+1)-1)`quot`(a-1)) fact

-- from 021
factorize :: Integer -> [(Integer, Int)]
factorize 1 = [(1, 0)]
factorize n = format $ factorize' n primes
    where
      factorize' n ps@(p : ps')
          | p * p > n    = [n]
          | rem n p == 0 = p : factorize' (div n p) ps
          | otherwise    = factorize' n ps'
      format ps = [(x, length xs) | xs@(x : _) <- group ps]
-- from 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

chainraw :: Integer -> [Integer]
chainraw n = imp n []
  where
    imp n _ | n>=(10^6) = [0]
    imp 1 _ = [1]
    imp 6 _ = [6]
    imp 28 _ = [28]
    imp 496 _ = [496]
    imp 8128 _ = [8128]
    imp 33550336 _ = [33550336]
    imp 8589869056 _ = [8589869056]
    imp n a@(x:xs) | n==x = a
                   | n`elem`a = [-1]
                   | otherwise = imp (dsum2 n) (a++[n])
    imp n a = imp (dsum n) (a++[n])
chain :: Integer -> (Int, Integer)
chain n = (length a,minimum a) where a = chainraw n

ans1imp :: (Int, Integer)
ans1imp = maximum [chain n | n <- [1..10^6]]

ans1 :: Integer
ans1 = let (a,b) = ans1imp in b

-- 14316
main :: IO ()
main = print ans1