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.

• user2305: This is not tail recursive. May give stackoverflow on very large number
• Will Ness: @user2305 no, Haskell is lazy. ... | n mod` d == 0 = d : factorize d (n div d) ...` is guarded recursion, as in "guarded by a lazy constructor (:)". it is perfectly fine and idiomatic, and won't cause any stack overflows.

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.

• Sᴀᴍ Onᴇᴌᴀ: Are you saying that the code by yasar (the OP) does not "get primes and get prime factors of a number."? or is that about the code you posted? Is that code supposed to be an improvement? If so, please explain why that is an improvement (unless your statement is intended to do that).

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