Prime factorization in Haskell

Question

I am a Haskell beginner. Here is my function to find prime factors of a number

primes = 2:takePrimes [3, 5 ..]
    where takePrimes (x:xs) = let smallPrimes = untilRoot x primes
                              in  if 0 `notElem` (map (mod x) smallPrimes)
                                  then x:takePrimes xs
                                  else takePrimes xs

untilRoot n = takeWhile (\x -> x*x < n)

firstPrimeDivisor :: Integer -> Integer
firstPrimeDivisor x = head [p | p <- primes, x `mod` p == 0]

primeFactor 1 = []
primeFactor x = let p = firstPrimeDivisor x
                in p:primeFactor (x `div` p)

How does it look?

Answer 1

It looks fine, but there are a few issues with it:

1. When we check if a number is a prime, we need to test all possible divisors up its root inclusively. If you print the primes list, you'll see that 9 is there, which is obviously a bug. untilRoot n = takeWhile (\x -> x*x < n) should be untilRoot n = takeWhile (\x -> x*x <= n). However, this bug does not affect the correctness of your factorization algorithm.

2. It does not scale well. In the firstPrimeDivisor function, all primess up to x are checked in the worst case. That's why it works slowly even for moderately large prime numbers(for instance, 10^9 + 7).

Here a more efficient solution which checks only O(sqrt(n)) divisors in the worst case:

factorize :: Integer -> Integer -> [Integer]
factorize _ 1 = [] 
factorize d n 
    | d * d > n = [n]
    | n `mod` d == 0 = d : factorize d (n `div` d)
    | otherwise = factorize (d + 1) n

primeFactors :: Integer -> [Integer]
primeFactors = factorize 2

It is also much more concise. The idea behind it is to prove two statements first and then write a very simple code based on them:

1. The smallest divisor of any natural number greater than two is a prime. Let's assume that it is not the case and d = p * q (p, q > 1), where d is the smallest divisor of n. But then p < d and p is a divisor of n. Thus, d is not the smallest one. This statement also shows why your original solution is correct, even though the primes list is contains some composite numbers.

2. Any composite number has a divisor that does not exceed its square root. The proof by contradiction is very straightforward so I will omit it.

Answer 2

Kraskevich's solution has one problem which is get the prime factors in set of 2,3,4,... In fact , we can get the primes in the set of 2,3,5,7,9...(2 and odd numbers), the bellow code avoids unnecessary operations on even numbers(greater than 2).

factors' :: Integral t => t -> [t]
factors' n
  | n < 0 = factors' (-n)
  | n > 0 = if 1 == n
               then []
               else let fac = mfac n 2 in fac : factors' (n `div` fac)
  where mfac m x
          | rem m x == 0 = x
          | x * x > m    = m -- if this line code can be matched, this is to say m can not be divided by 2,3,5,7,...n, n is the largest odd number less than sqrt(m). In other words, it is consistent with the definition of prime numbers.
          | otherwise    = mfac m (if odd x then x + 2 else x + 1) -- get factor in (2,3,5,7,9...)

In this problem, we should pay attention to

1. give a number n, x bellow to the set of 2,3,4,...,n-1(in above code, we use its subset(2 and even number)). if n is the smallest number that can be divisible by x，then n is always a prime number.
2. give a number n, we can define n = p1 * p2 * p3 *...* pi * pj * ... pi is prime, pi <= pj so we can get n/p1 = p2 * p3 * ...pi * pj..., what is p1? p1 is the number which gets in 1.

Answer 3

Here's a Sieve of Eratosthenes implementation from a literate Haskell page:

primes :: [Integer]
primes    = [2, 3, 5] ++ (diff [7, 9 ..] nonprimes)
  where
    nonprimes :: [Integer]
    nonprimes = foldr1 f $ map g $ tail primes
        where
            f (x:xt) ys = x : (merge xt ys)
            g p         = [ n * p | n <- [p, p + 2 ..]]

    merge :: (Ord a) => [a] -> [a] -> [a]
    merge xs@(x:xt) ys@(y:yt) =
      case compare x y of
        LT -> x : (merge xt ys)
        EQ -> x : (merge xt yt)
        GT -> y : (merge xs yt)

    diff :: (Ord a) => [a] -> [a] -> [a]
    diff xs@(x:xt) ys@(y:yt) =
      case compare x y of
        LT -> x : (diff xt ys)
        EQ -> diff xt yt
        GT -> diff xs yt