import System.Random

factors :: Int -> [Int]
factors n = factor [x | x <- [2..n], odd x, millerRabin x 1] n

factor :: [Int] -> Int -> [Int]
factor [] _ = []
factor (x:xs) n
    | mod n x == 0  = x:factor' (div n x) x
    | mod n x /= 0  = factor xs n

factor' :: Int -> Int -> [Int]
factor' n p
    | mod n p == 0          = p:factor' (div n p) p
    | mod n p /= 0 && n > 1 = [n]
    | otherwise             = []

millerRabin :: Int -> Int -> Bool
millerRabin n k
    | n < 1          = error "Invalid integer."
    | n > 3 && k > 0 = let a  = (\mi ma -> head $! take ma
                                (randomRs (mi, ma) (mkStdGen 100))) 2 (n - 2)
                           ds = findDS $! (n-1)
                           ad = (fromIntegral a) ^ fst ds
                           x :: Integer
                           x  = mod ad (fromIntegral n)
                        in if x == 1 || x == (fromIntegral $! (n - 1)) then
                                millerRabin n $! (k - 1)
                           else mrLoop (fromIntegral x) n k (snd ds) 1
    | otherwise = True

mrLoop :: Int -> Int -> Int -> Int -> Int -> Bool
mrLoop x n k s r
    | r < s     = let xx = mod (x * x) n in if xx == 1 then False
                     else if xx == (n - 1) then millerRabin n $! (k - 1)
                     else mrLoop xx n k s $! (r + 1)
    | otherwise = False

findDS :: Int -> (Int, Int)
findDS d
    | mod d 2 == 0 = let ds = findDS $! (div d 2) in (fst ds, (snd ds) + 1)
    | otherwise    = (d, 0)

concatPrimes :: [Int] -> Int
concatPrimes []        = 0
concatPrimes (x:xs)
    | millerRabin x 1  = let ys = concatPrimes xs in if ys == 0 then x
                            else read ((show x) ++ (show ys))
    | otherwise        = concatPrimes xs

homePrime :: Int -> Int
homePrime n = let xs = factors n in if length xs > 1 then
                    homePrime $ concatPrimes xs else n

euclidMullin :: Int -> [Int]
euclidMullin 1 = [2]
euclidMullin n = let xs = euclidMullin (n-1)
                 in xs ++ [((\p -> head $ factors p) (product xs + 1))]

main=do
    let x = 9
    let y = 7
    putStrLn $ "Home Prime of " ++ (show x) ++ ": " ++ (show $ homePrime x)
    putStrLn $ "Euclid-Mullin sequence of " ++ (show y) ++ " " ++
        (show $ euclidMullin y)