Primes, factorization, and Euler's totient function

Haskell

Public Domain

Download (right click, save as, rename as appropriate)

Embed

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
primes :: Integral a => [a]
primes = 2 : filter (isPrime 0) [3, 5..]
 where isPrime i n
	| squares !! i > n         = True
	| mod n (primes !! i) == 0 = False
	| otherwise                = isPrime (i + 1) n
	where squares = map (^ 2) primes

-- |@factor n@ returns a list of pairs of prime factors of @n@ and their
-- multiplicities.  If @n@ is @1@, an empty list is returned.  If @n@ is @0@,
-- @[(0, 1)]@ is returned.  If @n@ is negative, the result is equivalent to
-- @(-1, 1) : factor (-n)@.
--
-- This function is not recommended for use on large integers.
factor :: Integral a => a -> [(a, Int)]
factor = fst . factorWith primes

-- |@factorWith primes n@ decomposes @n@ into a product of powers of elements
-- of @primes@ (which are assumed to all be coprime and positive) and a
-- coprime quotient.  Specifically, it returns a pair containing:
--
-- * a list of pairs of factors of @n@ in @primes@ and their multiplicities,
--   and
--
-- * the product of any remaining factors of @n@ that do not appear in
--   @primes@ (or @1@ if there are none).
--
-- If @n@ is @1@, @([], 1)@ is always returned.  If @n@ is @0@, @([(0, 1)],
-- 1)@ is always returned.  If @n@ is negative, the result is the same as
-- @factorWith primes (-n)@ but with @(-1, 1)@ prepended to the list of
-- factors.
factorWith :: Integral a => [a] -> a -> ([(a, Int)], a)
factorWith _ 0 = ([(0, 1)], 1)
factorWith primal n
 | n < 0 = let (f, cp) = factor' (filter (\p -> mod (-n) p == 0) primal) (-n)
	   in ((-1, 1) : f, cp)
 | n > 0 = factor' (filter (\p -> mod n p == 0) primal) n
 where factor' _  1    = ([], 1)
       factor' [] n    = ([], n)
       factor' (p:q) n = let (k, m) = until (\(_, m) -> mod m p /= 0)
				       (\(k, m) -> (k+1, div m p)) (0, n)
			     (a, b) = factor' q m
			 in ((p, k) : a, b)

-- |Euler's totient phi function.  Returns the number of positive integers
-- less than or equal to @n@ that are coprime to @n@.
totient :: Integral a => a -> a
totient n | n <= 0 = undefined
totient n = product [p^k - p^(k-1) | (p, k) <- factor n]