General-purpose algorithms and data structures, illustrated in Haskell. Part II
We show a simple example of achieving all the benefits of an imperative data structure -- including sharing and the efficiency of updates -- in a pure functional program. Our data structure is a doubly-linked, possibly cyclic list, with the standard operations of adding, deleting and updating elements; traversing the list in both directions; iterating over the list, with cycle detection. The code:
IORef
, STRef
, TVars
, or any other destructive updates;The algorithm is essentially imperative -- hence supporting the identity check and the in-place `update' -- but implemented purely functionally. Although the code relies on many local, type safe `heaps', there is emphatically no global heap and no global state.
It is not for nothing that Haskell has been called the best imperative language. Imperative algorithms are implementable as they are -- yet genuinely functionally, without resorting to the monadic sub-language but taking the full advantage of clausal definitions, pattern guards, laziness.
Haskell-Cafe discussion ``Updating doubly linked lists''. January 2009.
[0,1]
. Mart'n Escardo's technique can tell, in finite
time, if a given total computable predicate is satisfied
over all possible infinite bit strings. Furthermore, for
so-called sparse predicates, Marti'n Escardo's technique is very
fast.We re-formulate the problem in terms of streams and depth-limited depth-first search, and thus cast off the mystery of deciding the satisfiability of a total computable predicate over the set of all Cantor numerals, which are uncountable.
As an additional contribution, we show how to write functions over Cantor numerals in a `natural' monadic style so that those functions become self-partially evaluating. The instantiation of the functions in an appropriate pure monad gives us transparent memoization, without any changes to the functions themselves. The monad in question is pure and involves no reference cells.
On `dense' functions on numerals (i.e., those that need to examine most of the bits of their argument, up to a limit), our technique performs about 9 times faster than the most sophisticated one by Marti'n Escardo'.
StreamPEval.hs [12K]
Extensively commented Haskell98 code
Marti'n Escardo': Seemingly impossible functional programs.
< http://math.andrej.com/2007/09/28/seemingly-impossible-functional-programs/ >
We limit attention to closed terms since all interesting terms mentioned earlier are closed. De Bruijn indices let us avoid wasting time on building alpha-equivalent terms. (The terms will be pretty-printed in the conventional, named-variable notation, for clarity.) Thus we will be generating terms described by the following grammar (data type declaration), with the side condition that the terms are closed.
data Exp = V Int | A Exp Exp | L Exp
The generator is a truly straightforward application of the simplest
non-deterministic operations in the MonadPlus
interface.
a_term_naive :: MonadPlus m => m Exp a_term_naive = L `liftM` go 0 where -- go n: generate a lambda term with free variables 0..n go n = choose_var n `mplus` choose_lam n `mplus` choose_app n choose_var 0 = return (V 0) choose_var n = return (V n) `mplus` choose_var (n-1) choose_app n = liftM2 A (go n) (go n) choose_lam n = L `liftM` go (n+1)
Absent constants, a closed lambda-term must start with a lambda-binder. The rest could be either a variable, an application of two non-deterministically chosen terms, or an abstraction with a random body.
The equally straightforward implementation of MonadPlus
, the List monad,
fails to generate anything interesting within the first 10 000 terms. It is
not difficult to see why. Here are the first 8 generated terms:
Lx.x, Lx.Ly.x, Lx.Ly.y, Lx.Ly.Lz.x, Lx.Ly.Lz.y, Lx.Ly.Lz.z, Lx.Ly.Lz.Lu.x, Lx.Ly.Lz.Lu.y
The List monad realizes an incomplete search strategy: there is no guarantee that any given term will ever come. No applications are indeed forthcoming with the List monad. A better implementation of non-determinism is needed, with a complete search strategy: for example, iterative deepening or FBackTrack. With iterative deepening (which in our implementation produces the same sequence as breadth-first search but without taking Gbytes of memory), the first 10 generated terms
Lx.x, Lx.Ly.x, Lx.Ly.y, Lx.x x, Lx.Ly.Lz.x, Lx.Ly.Lz.y, Lx.Ly.Lz.z, Lx.x (Ly.x), Lx.x (Ly.y), Lx.x (x x)
look quite hopeful. Within the first 10 000 generated terms,
the self-application Lx.x x
, or delta
, comes 4th; the third Church
numeral comes 716th and omega
comes 3344th. Alas, there is no
successor or the S combinator, let alone the Y combinator. With
the FBackTrack implementation, the first few terms may look even more hopeful
Lx.x, Lx.x x, Lx.Ly.x, Lx.x x x, Lx.Ly.x x, Lx.x (x x), Lx.Ly.Lz.x, Lx.(Ly.x) x, Lx.x (Ly.x)The term
delta
comes the second, the third Church numeral comes
695th. Alas, within the first 99 977 terms nothing else interesting
comes. Generating interesting terms is not at all as straightforward
as it seems. Even sophisticated MonadPlus implementations did not
help.
It is crucial to recognize that the straightforward search for lambda-terms
is biased, and that bias is not in favor of interesting terms. Generally,
the fewer
non-deterministic choices it takes to produce a term, the sooner the term
comes. The straightforward generator clearly
favors selectors (terms of the form Lx.Ly...Lu.x
) and simple
applications such as Lx. (x x) (x x)
. We should put the brake on
generating abstractions: each new L
adds a new variable and hence
dramatically many more possible terms. We should encourage
the generator to play more with the variables it already has.
To prevent the string of applications of a variable to itself,
we should produce terms in the general order of their size. Again the goal
is to explore the search space more uniformly. It is tempting to generate
only normal forms. Alas, the Y combinator and omega
do not have normal
forms. Therefore, we do generate redices, but only of interesting kinds
(whose argument is an abstraction), and only occasionally. The following code
incorporates these ideas:
a_term :: MonadPlus m => m Exp a_term = L `liftM` go 0 where -- go n: generate a lambda term with free variables 0..n go n = do size <- iota 1 gen n True size -- gen n l s: generate a lambda term with free variables 0..n and -- of the size exactly s. The second argument tells if to generate -- abstractions gen _ _ s | s <= 0 = mzero gen n _ 1 = choose_var n gen n True 2 = choose_lam n 2 gen n False 2 = mzero gen n True s = choose_app n s `mplus` choose_lam n s gen n False s = choose_app n s choose_var n = msum . map (return . V) $ [0..n] choose_lam n s = penalty (40*n) $ L `liftM` gen (n+1) True (s-1) choose_app n s = do let s' = s - 1 -- Account for the 'A' constructor let gen_redex = do lefts <- range 4 (s' - 3) liftM2 A (choose_lam n lefts) (choose_lam n (s' - lefts)) let gen_noredex = do lefts <- range 1 (s' - 1) liftM2 A (gen n False lefts) (gen n True (s'-lefts)) gen_noredex `mplus` penalty 4 gen_redex
The auxiliary iota i
produces an integer greater or equal to i
; range i j
chooses an integer between i
and j
, inclusive. The
operation yield
-- Lowers the priority of m, so choices of m will be tried less often yield :: MonadPlus m => m a -> m a yield m = mzero `mplus` mand its n-th iterate
penalty n
produce a string of failures before trying
the argument computation m
. When a complete search strategy sees
many failures it tends to turn away and to pay more attention to other parts
of the search space.Here are the first 10 terms produced by the sophisticated generator:
Lx.x, Lx.Ly.y, Lx.Ly.x, Lx.x x, Lx.x (Ly.y), Lx.Ly.y y, Lx.Ly.x y, Lx.x (x x), Lx.Ly.y x, Lx.x (Ly.x)The term
delta
comes fourth, omega
comes 54th, the Y combinator 303d,
the third Church numeral 393d, the S combinator 1763d and the successor
4865th.The conclusion, although obvious in hindsight, is still thought-provoking: interesting lambda-terms are really hard to encounter by accident. They are exquisitely rare in the space of possible lambda-terms and distributed non-uniformly. A monkey banging on even the sophisticated lambda-typewriter may have printed the Y combinator, but would unlikely to print even the addition combinator within its lifetime.
FBackTrack.hs [3K]
Simple fair and terminating backtracking monad
Searches.hs [9K]
Iterative deepening search
n
in base b
is the integer l
such that b^l <= n < b^(l+1)
. The number one greater than the integer
logarithm base b
is the size of n
, that is, the number of digits,
in its base-b
representation. We present
the Haskell98 code that is just as efficient as the internal GHC.Integer.Logarithms.integerLogBase#
function but uses no unboxed
data, no optimizations, and is not even compiled.
Naively, the integer logarithm is computed by repeated divisions of n
by b
, until the result reaches 1. This procedure requires l
divisions, where l
is the logarithm. We present two efficient algorithms:
the first uses only multiplications, no more than 2*log_2 l
of them, and (log_2 l) * sizeof(n)
extra memory for the powers of b
. The second
algorithm does log_2 l
multiplications and no more than log_2 l
integer divisions (and the same amount of extra memory) to compute l
.
Here is the multiplication-only procedure that returns l+1
where l
is the integer logarithm of n
in base b
. It is a composition
of two functions, which together compute the bits of l
: major_bit
determines
the upper bound and other_bits
improves it.
data BaseP = BaseP !Integer -- b^k !Int -- k mul_bp :: BaseP -> BaseP -> BaseP mul_bp (BaseP bk1 k1) (BaseP bk2 k2) = BaseP (bk1*bk2) (k1+k2) basewidth :: Integer -> Integer -> Int basewidth b _ | b < 1 = error "basewidth: base must be greater than 1" basewidth b n | n < b = 1 basewidth b n | n == b = 2 basewidth b n = major_bit [BaseP b 1] where major_bit :: [BaseP] -> Int major_bit bases@(bp:bps) = let bpnext@(BaseP bk2 k2) = bp `mul_bp` bp in case compare bk2 n of EQ -> k2 + 1 -- n == b^(2k) GT -> other_bits bp bps -- b^(2k) > n LT -> major_bit (bpnext : bases) -- b^(2k) < n other_bits (BaseP _ i) [] = i+1 -- b^i < n < b^(i+1) other_bits bp (bphalf:bps) = let bpup@(BaseP bik ik) = bp `mul_bp` bphalf in case compare bik n of EQ -> ik + 1 -- n == b^(i+k) GT -> other_bits bp bps -- b^i < n < b^(i+k) LT -> other_bits bpup bps -- b^(i+k) < n < b^(i+2k)The correctness and the complexity analysis follow from the invariants of the two auxiliary functions. In any call
major_bit bases
, the list bases
is [BaseP b^k k | j <- [d,d-1..0], k=2^j]
for some d
and n > b^(2^d)
. In the first
invocation, d
is 0 and progressively increases until such d>0
is
found that b^(2^(d-1)) < n <= b^(2^d)
. In any call other_bits
(BaseP bi i) bases
, the list bases
is [BaseP b^k k | j <- [d,d-1..0], k=2^j]
for some d
, k=2^d
and b^i < n <
b^(i+2k)
. Each invocation of other_bits
halfs k
until it reaches
one. The other, more efficient as it turns out, algorithm modifies other_bits
to divide n
by the candidate lower bound b^i
.
That integer logarithm function computes that the 48th Mersenne prime 2^57885161-1
has 17425170 digits in its decimal representation --
in 1.2
seconds.
r^n
for all r>=2
and n>1
.
Thanks to lazy evaluation, the table is automatically sized as needed.
There is no need to guess the maximal size of the table so to allocate it.powers = [2..] : map (zipWith (*) (head powers)) powersThe rest of the algorithm exhibits similar stream-wise processing and computations of tables in terms of themselves.
Messages speedup help by Damien R. Sullivan, Bill Wood, Andrew J Bromage, Mark P Jones, and many others posted on the Haskell-Cafe mailing list on March 6-8, 2003
< http://www.haskell.org/pipermail/haskell-cafe/2003-March/004065.html >
< http://www.haskell.org/pipermail/haskell-cafe/2003-March/004075.html >
s f g x = f x (g x) puzzle = (!!) $ iterate (s (lzw (+)) (0:)) [1] where lzw op [] ys = ys lzw op (x:xs) (y:ys) = op x y : lzw op xs ys
Incidentally, a small change gives a different series:
puzzle1 = (!!) $ iterate (s ((lzw (+)).(0:)) (1:)) [] where lzw op [] ys = ys lzw op (x:xs) (y:ys) = op x y : lzw op xs ys
Finally, how can we possibly live without the following:
puzzle2 = (!!) $ iterate (s ((lzw (+)).(1:).(0:)) (0:)) [1,1] where lzw op xs [] = [] lzw op (x:xs) (y:ys) = op x y : lzw op xs ys
Hints:
*Main> puzzle 5 [1,5,10,10,5,1] *Main> puzzle1 5 [1,2,4,8,16] *Main> puzzle2 5 [1,1,2,3,5,8,13]
Discussion thread, started on Haskell-Cafe by Thomas Bevan on Aug 22, 2003.
< http://www.haskell.org/pipermail/haskell-cafe/2003-August/004977.html >
oleg-at-pobox.com or oleg-at-okmij.org
Your comments, problem reports, questions are very welcome!
Converted from HSXML by HSXML->HTML