Say that we want to compute the 3-smooth numbers,
i.e., the integers that can be written in the form
It's easy enough to compute the powers of 2 and 3:
powersOf2 :: [Integer]
powersOf2 = iterate (*2) 1
powersOf3 :: [Integer]
powersOf3 = iterate (*3) 1But how do we implement the following?
-- take 10 smooth3 == [1, 2, 3, 4, 6, 9, 12, 16, 18, 24]
smooth3 :: [Integer]
smooth3 = ..If we limited powersOf2 and powersOf3, say to the first 10 elements, we
could write:
powersOf2Bounded :: [Integer]
powersOf2Bounded = take 10 (iterate (*2) 1)
powersOf3Bounded :: [Integer]
powersOf3Bounded = take 10 (iterate (*3) 1)
smooth3Bounded :: [Integer]
smooth3Bounded = sort (liftA2 (*) powersOf2Bounded powersOf3Bounded)And smooth3Bounded would consist of powersOf2 and powersOf3 are infinite. We need a
way to lazily generate smooth3 given powersOf2 and powersOf3.
Let's lay out the elements of smooth3 in a 2D grid:
| 1 2 4 8 ..
---------------------------------
1 | 1 2 4 8 ..
|
3 | 3 6 12 24 ..
|
9 | 9 18 36 72 ..
|
27 | 27 54 108 216 ..
|
: | : : : : ॱ.
Since (*) is monotonically increasing in both arguments, each element in the
grid is ≥ its neighbors above and to the left. Therefore, when an element
appears in smooth3, its up- and left-neighbors must already have appeared.
Let's think about smooth3 after 3 elements have been produced:
| 1 2 4 8 .. --------------------------------- 1 |124 8 .. | 3 |36 12 24 .. | 9 | 9 18 36 72 .. | 27 | 27 54 108 216 .. | : | : : : : ॱ.
After producing 1, 2, 3, the next element in smooth3 can only be one of
{4, 6, 9}. We know this without performing any comparisons, just by the
positions of these elements in the grid, as these are the only elements whose
up- and left-neighbors have already been produced.
After producing the next element, 4, our grid becomes
| 1 2 4 8 .. --------------------------------- 1 |1248 .. | 3 |36 12 24 .. | 9 | 9 18 36 72 .. | 27 | 27 54 108 216 .. | : | : : : : ॱ.
Now, the "frontier" of possible next elements is {6, 8, 9}. We added 8 to
the frontier as it is no longer "blocked" by 4. Note that we cannot yet add
12 to the frontier even though it's no longer blocked by 4, as it's still
blocked by 6.
After producing the next element, 6, our grid becomes
| 1 2 4 8 .. --------------------------------- 1 |1248 .. | 3 |3612 24 .. | 9 | 9 18 36 72 .. | 27 | 27 54 108 216 .. | : | : : : : ॱ.
Now, the frontier consists of {8, 9, 12}. We added 12 as it is now unblocked
by 4 and 6. We cannot yet add 18, as it is still blocked by 9.
This investigation suggests the following algorithm for producing smooth3:
1. Initialize a "frontier" with the top left element
2. Delete the smallest element z from the frontier and yield it as the next element of smooth3.
3. If the down- or right-neighbors of z are unblocked, add them to the frontier.
4. Go to step 2.
The idea in the algorithm above is quite general, and instead of using it on
(*), powersOf2, and powersOf3, we can apply it to any binary function f
that is non-decreasing in both arguments, and any (potentially infinite) ordered
lists xs and ys. Then we can define a general
applyMerge :: (Ord c) => (a -> b -> c) -> [a] -> [b] -> [c]where applyMerge f xs ys is an ordered list of all f x y, for each x in xs
and y in ys. In particular, smooth3 = applyMerge (*) powersOf2 powersOf3.
We can use a priority queue to maintain the frontier of elements. To determine
in step 3 whether a down- or right-neighbor is unblocked, we maintain two
IntSets, for the x- and y- indices of all elements in the frontier. Then
in step 3, if the x- and y-index of a neighbor are not in the respective
IntSets, the neighbor is unblocked and can be added to the frontier.
After producing IntSets to produce the next element takes
applyMerge f xs ys takes f and compare take
Note that applyMerge retains the input lists in memory, which could cause
unexpected memory usage when the input lists are lazily generated. For example,
sum (take n (applyMerge const [1 :: Int ..] [1 :: Int ..]))
requires retaining the first
sum (take n (applyMerge (+) [1 :: Int ..] [1 :: Int ..]))
which requires retaining the first
With applyMerge, we can implement a variety of complex algorithms succinctly.
For example, the Sieve of Erastosthenes1 to generate prime numbers:
primes :: [Int]
primes = 2 : ([3..] `minus` composites) -- `minus` from data-ordlist
composites :: [Int]
composites = applyMerge (*) primes [2..]Gaussian integers, ordered by norm (Wikipedia):
zs :: [Integer]
zs = 0 : concatMap (\i -> [i, -i]) [1..]
gaussianIntegers :: [GaussianInteger] -- `GaussianInteger` from arithmoi
gaussianIntegers = map snd (applyMerge (\x y -> (norm (x :+ y), x :+ y)) zs zs)Square-free integers (Wikipedia):
squarefrees :: [Int]
squarefrees = [1..] `minus` applyMerge (*) (map (^2) primes) [1..]In data-ordlist,
there is mergeAll :: Ord a => [[a]] -> [a],
which merges a potentially infinite list of ordered lists, where the heads of
the inner lists are sorted. This is more general than applyMerge, in the sense
that we can implement applyMerge in terms of mergeAll:
applyMerge :: (Ord c) => (a -> b -> c) -> [a] -> [b] -> [c]
applyMerge f xs ys =
mergeAll (map (\y -> map (\x -> f x y) xs) ys)However, mergeAll uses applyMerge uses just
In The Algorithm Design Manual, Steven Skiena describes an algorithm for minimizing the sum of two airline ticket fares:
“Got it!,” I said. “We will keep track of index pairs in a priority queue, with the sum of the fare costs as the key for the pair. Initially we put only pair (1, 1) on the queue. If it proves it is not feasible, we put its two successors on—namely (1, 2) and (2, 1). In general, we enqueue pairs (i + 1, j) and (i, j + 1) after evaluating/rejecting pair (i, j). We will get through all the pairs in the right order if we do so.”
The gang caught on quickly. “Sure. But what about duplicates? We will construct pair (x, y) two different ways, both when expanding (x − 1, y) and (x, y −1).”
“You are right. We need an extra data structure to guard against duplicates. The simplest might be a hash table to tell us whether a given pair exists in the priority queue before we insert a duplicate. In fact, we will never have more than n active pairs in our data structure, since there can only be one pair for each distinct value of the first coordinate.”
This is similar to the applyMerge algorithm, except that applyMerge has an
optimization to check that we don’t add (x, y) to the priority queue when there
is already an (x′, y) with x < x′ or (x, y′) with y < y′ in the queue.
Footnotes
-
Note that this is really the Sieve of Erastosthenes, as defined in the classic The Genuine Sieve of Eratosthenes. Constrast this to other simple prime generation implementations, such as
primes = sieve [2..] where sieve (p : xs) = p : sieve [x | x <- xs, x `rem` p > 0]
which is actually trial division and not a faithful implementation of the Sieve of Erastosthenes. ↩