Prime factors in Haskell

Question

I'm new to Haskell.

How to generate a list of lists which contains prime factors of next integers?

Currently, I only know how to generate prime numbers:

primes = map head $ iterate (\(x:xs) -> [y | y<-xs, y `mod` x /= 0 ]) [2..]

Hangup

How far have you got? For example you could try writing a function to factorize a single number. – Harald

Answer 1

A simple approach to determine the prime factors of n is to

• search for the first divisor d in [2..n-1]
• if D exists: return d : primeFactors(div n d)
• otherwise return n (since n is prime)

Code:

prime_factors :: Int -> [Int]

prime_factors 1 = []
prime_factors n
  | factors == []  = [n]
  | otherwise = factors ++ prime_factors (n `div` (head factors))
  where factors = take 1 $ filter (\x -> (n `mod` x) == 0) [2 .. n-1]

This obviously could use a lot of optimization (search only from 2 to sqrt(N), cache the prime numbers found so far and compute the division only for these etc.)

UPDATE

A slightly modified version using case (as suggested by @user5402):

prime_factors n =
  case factors of
    [] -> [n]
    _  -> factors ++ prime_factors (n `div` (head factors))
  where factors = take 1 $ filter (\x -> (n `mod` x) == 0) [2 .. n-1]

Decoteau

instead of factors == [], it is more idiomatic (and perhaps more efficient) to use case factors of [] -> ... – Sandwich

note that foo == [] introduces an Eq constraint in the type of foo (which is actually not needed here) while case foo of [] -> ... does not. – Xyster

An alternative to factors == [] but still using guards would be to just use null factors instead. – Proud

there is no need too check all [2..n-1] numbers. You can only check all [2..sqrt(n)], because if p | n, then p <= sqrt(n). – Kramlich

@RottenBrain I know (I mentioned it in my answer as a possible optimization) - I didn't want to spoil the fun for the OP :-) – Decoteau

@Frank Schmitt, oh sorry. You really have mentioned this optimisation :) – Kramlich

without the sqrt it's just terrible, and even with it it's suboptimal. – Rampageous

@WillNess Then please go ahead and post your own optimized solution :-) – Decoteau

please see https://mcmap.net/q/743357/-my-haskell-solution-to-euler-3-is-inefficient. :) Also, haskell.org/haskellwiki/…. (and much more verbose https://mcmap.net/q/743358/-racket-programming-where-am-i-going-wrong) – Rampageous

@WillNess Impressive. Although I think that for repeated isPrime() calls for large numbers, a solution with memoization of the prime numbers found so far should be faster. – Decoteau

(also, check this out: #24003819 - a nice one-liner, though also suboptimal). – Rampageous

btw your code can be written shorter and equivalently as unfoldr (\n -> listToMaybe [(x, div n x) | x<- [2..n], mod n x==0]). – Rampageous

@Kramlich "because if p | n, then p <= sqrt(n)", but number 5424684 has a prime factor 6367 which is larger than (sqrt 5424684 = ~2329). Or did I miss the point? – Oaken

@Dai A number may very well have a prime factor that is larger than its square root. But: if a number n has a prime factor f1 which is greater than its square root, it must have at least one other prime factor f2 which is smaller than its square root (since n = f1*f2*<other factors>, and if f2 was also larger than sqrt(n), then the following would hold: f1*f2 > sqrt(n) * f2 > sqrt(n) * sqrt(n) = n, which contradicts the assumption that f1 and f2 are prime factors of n. – Decoteau

@FrankSchmitt Thanks for the explanation! So we recursively divide n by the factor smaller than sqrt n, thus never need to check numbers beyond sqrt n. – Oaken

@Dai Exactly. this_comment_was_too_short – Decoteau

Answer 2

This is a good-performanced and easy-to-understand implementation, in which isPrime and primes are defined recursively, and primes will be cached by default. primeFactors definition is just a proper use of primes, the result will contains continuous-duplicated numbers, this feature makes it easy to count the number of each factor via (map (head &&& length) . group) and it's easy to unique it via (map head . group) :

isPrime :: Int -> Bool
primes :: [Int]

isPrime n | n < 2 = False
isPrime n = all (\p -> n `mod` p /= 0) . takeWhile ((<= n) . (^ 2)) $ primes
primes = 2 : filter isPrime [3..]

primeFactors :: Int -> [Int]
primeFactors n = iter n primes where
    iter n (p:_) | n < p^2 = [n | n > 1]
    iter n ps@(p:ps') =
        let (d, r) = n `divMod` p
        in if r == 0 then p : iter d ps else iter n ps'

And the usage:

> import Data.List
> import Control.Arrow

> primeFactors 12312
[2,2,2,3,3,3,3,19]

> (map (head &&& length) . group) (primeFactors 12312)
[(2,3),(3,4),(19,1)]

> (map head . group) (primeFactors 12312)
[2,3,19]

Parados

Answer 3

Until the dividend m < 2,

1. take the first divisor n from primes.
2. repeat dividing m by n while divisible.
3. take the next divisor n from primes, and go to 2.

The list of all divisors actually used are prime factors of original m.

Code:

-- | prime factors
--
-- >>> factors 13
-- [13]
-- >>> factors 16
-- [2,2,2,2]
-- >>> factors 60
-- [2,2,3,5]
--
factors :: Int -> [Int]
factors m = f m (head primes) (tail primes) where
  f m n ns
    | m < 2 = []
    | m `mod` n == 0 = n : f (m `div` n) n ns
    | otherwise = f m (head ns) (tail ns)

-- | primes
--
-- >>> take 10 primes
-- [2,3,5,7,11,13,17,19,23,29]
--
primes :: [Int]
primes = f [2..] where f (p : ns) = p : f [n | n <- ns, n `mod` p /= 0]

Update:

This replacement code improves performance by avoiding unnecessary evaluations:

factors m = f m (head primes) (tail primes) where
  f m n ns
    | m < 2 = []
    | m < n ^ 2 = [m]   -- stop early
    | m `mod` n == 0 = n : f (m `div` n) n ns
    | otherwise = f m (head ns) (tail ns)

primes can also be sped up drastically, as mentioned in Will Ness's comment:

primes = 2 : filter (\n-> head (factors n) == n) [3,5..]

Palua

this code is slow twice over: in factors it doesn't stop at sqrt, and in primes it doesn't wait until sqr -- of a prime. – Rampageous

@WillNess Thanks to your advice, I can update code. This code become faster significantly by adding only a guard. Original code takes 15 seconds to factorize 599946, while updated code takes only 6 milliseconds. – Palua

great; you've improved your factors but not the primes which still is much too slow. compare your definition's speed with that of 2:filter (\n-> head (factors n) == ...... (complete the definition). :) can you find the reason for such a difference in speed? see also "empirical orders of growth". – Rampageous

@WillNess I understood what you intend to do, but I could not complete your code snippet. I feel sorry but can you show me primes you defined? – Palua

sure, it's primes = 2:filter (\n-> head (factors n) == n) [3,5..]. :) – Rampageous

Answer 4

Haskell allows you to create infinite lists, that are mutually recursive. Let's take an advantage of this.

First let's create a helper function that divides a number by another as much as possible. We'll need it, once we find a factor, to completely eliminate it from a number.

import Data.Maybe (mapMaybe)

-- Divide the first argument as many times as possible by the second one.
divFully :: Integer -> Integer -> Integer
divFully n q | n `mod` q == 0   = divFully (n `div` q) q
             | otherwise        = n

Next, assuming we have somewhere the list of all primes, we can easily find factors of a numbers by dividing it by all primes less than the square root of the number, and if the number is divisible, noting the prime number.

-- | A lazy infinite list of non-trivial factors of all numbers.
factors :: [(Integer, [Integer])]
factors = (1, []) : (2, [2]) : map (\n -> (n, divisors primes n)) [3..]
  where
    divisors :: [Integer] -> Integer -> [Integer]
    divisors _ 1          = []   -- no more divisors
    divisors (p:ps) n
        | p^2 > n       = [n]  -- no more divisors, `n` must be prime
        | n' < n        = p : divisors ps n'    -- divides
        | otherwise     = divisors ps n'        -- doesn't divide
      where
        n' = divFully n p

Conversely, when we have the list of all factors of numbers, it's easy to find primes: They are exactly those numbers, whose only prime factor is the number itself.

-- | A lazy infinite list of primes.
primes :: [Integer]
primes = mapMaybe isPrime factors
  where
    -- |  A number is prime if it's only prime factor is the number itself.
    isPrime (n, [p]) | n == p  = Just p
    isPrime _                  = Nothing

The trick is that we start the list of factors manually, and that to determine the list of prime factors of a number we only need primes less then its square root. Let's see what happens when we consume the list of factors a bit and we're trying to compute the list of factors of 3. We're consuming the list of primes, taking 2 (which can be computed from what we've given manually). We see that it doesn't divide 3 and that since it's greater than the square root of 3, there are no more possible divisors of 3. Therefore the list of factors for 3 is [3]. From this, we can compute that 3 is another prime. Etc.

Sultry

Answer 5

I just worked on this problem. Here's my solution.

Two helping functions are

factors n = [x | x <- [1..n], mod n x == 0]
isPrime n = factors n == [1,n]

Then using a list comprehension to get all prime factors and how many are they.

prime_factors num = [(last $ takeWhile (\n -> (x^n) `elem` (factors num)) [1..], x) | x <- filter isPrime $ factors num]

where

x <- filter isPrime $ factors num

tells me what prime factors the given number has, and

last $ takeWhile (\n -> (x^n) `elem` (factors num)) [1..]

tells me how many this factor is.

Examples

> prime_factors 36    -- 36 = 4 * 9
[(2,2),(2,3)]

> prime_factors 1800  -- 1800 = 8 * 9 * 25
[(3,2),(2,3),(2,5)]

Guadiana

Answer 6

More elegant code，use 2 and odd numbers to divide the number.

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
          | otherwise    = mfac m (if odd x then x + 2 else x + 1)

Feudatory

let fac = mfac n 2 in fac : factors' (n `div` fac) is suboptimal, always starts from 2, can start from the previous fac. – Rampageous

Answer 7

Here's my version. Not as concise as the others, but I think it's very readable and easy to understand.

import Data.List

factor :: Int -> [Int]
factor n
  | n <= 1 = []
  | even n = 2 : factor(div n 2)
  | otherwise =
    let root = floor $ sqrt $ fromIntegral n
    in
      case find ((==) 0 . mod n) [3, 5.. root] of
        Nothing  -> [n]
        Just fac -> fac : factor(div n fac)

Selfdeprecating

Answer 8

I'm sure this code is ugly enough to drive a real Haskell programmer to tears, but it works in GHCI 9.0.1 to provide prime factors with a count of each prime factor.

import Data.List
factors n = [x | x <- [2..(n`div` 2)], mod n x == 0] ++ [n]
factormap n = fmap factors $ factors n
isPrime n = case factormap n of [a] -> True; _ -> False
primeList (x:xs) = filter (isPrime) (x:xs)
numPrimes n a = length $ (factors n) `intersect` (takeWhile ( <=n) $ iterate (a*) a)
primeFactors n = primeList $ factors n
result1 n = fmap (numPrimes n) (primeFactors n)
answer n =  ((primeFactors n),(result1 n))

Example: ghci> answer 504

([2,3,7],[3,2,1])

The answer is a list of prime factors and a second list showing how many times each prime factor is in the submitted number.

Legging

Answer 9

How about this piece of code? I'm afraid that I'm going through unnecessary evaluations in order to reach prime numbers of the form 6*m+5 or 6*m+7 after the first two, but that is the price of writing a recursive solution that may need repeated numbers on the result without having a prior list of primes. In other words, if I don't have their list, then I can check the natural number against their form.

primeFactors :: Integer -> [Integer]
primeFactors n | n <= 1 = []
               | otherwise = k : primeFactors (n `div` k)
  where
    k | n `mod` 2 == 0 = 2
      | n `mod` 3 == 0 = 3
      | otherwise = go 0
        where
          go m | n `mod` (6*m+5) == 0 = 6*m+5
               | n `mod` (6*m+7) == 0 = 6*m+7
               | otherwise = go (m+1)

---

Why this form though? You can write natural numbers using the first, 2*n, 6*n+3, 6*n+5 and 6*n+7, but two of these are divisble by the first two primes, while the others are only divisible by the following two under special conditions, which I used in this solution.

Scanty