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]