Archive for April, 2012

Perfect shuffles in Haskell

April 17, 2012

Consider a deck of 2 n cards. There exists a bijection from the set of cards to the finite set \{1, 2, \dots, 2 n\}, which allows us to represent the deck using a (2 n)-tuple, as follows

(1, 2, \dots, n, n+1, n+2, \dots, 2 n)

where the first element represents the card at the top of the deck, and the (2 n)-th element represents the card at the bottom of the deck. A perfect shuffle consists of cutting the deck exactly in half and then interlacing / interleaving the two halves perfectly. The in shuffle moves the original top card to the second position. The out shuffle leaves the original top card on top [1]. After a perfect in shuffle the deck will be

(n+1, 1, n+2, 2, \dots, 2 n, n),

whereas after a perfect out shuffle the deck will be

(1, n+1, 2, n+2, \dots, n, 2 n).

Is it possible to restore the deck to its original order after a certain number of perfect shuffles? Amazingly (or perhaps not), it is indeed possible!

__________

Example 1

Suppose we are given a deck of 8 cards (i.e., n = 4). A perfect out shuffle of this deck is depicted below

[ source ]

Using pencil and paper, it is easy to see that after three perfect out shuffles the deck will be restored to its original order. Here is a Haskell script:

-- create Card data type
data Card = Ace | Two | Three | Four | 
            Five | Six | Seven | Eight 
            deriving (Show, Eq)

-- create type synonym
type Deck = [Card]

-- interlace two decks of the same size
interlace :: Deck -> Deck -> Deck
interlace [] [] = []
interlace (c1:d1) (c2:d2) = [c1,c2] ++ interlace d1 d2

-- perfect out shuffle
outShuffle :: Deck -> Deck
outShuffle deck = interlace deck1 deck2
                  where deck1 = fst (splitAt 4 deck)
                        deck2 = snd (splitAt 4 deck)

-- perfect in shuffle
inShuffle :: Deck -> Deck
inShuffle deck = interlace deck2 deck1
                 where deck1 = fst (splitAt 4 deck)
                       deck2 = snd (splitAt 4 deck)

where I used function splitAt to cut the deck into two halves of equal size. Since I am not acquainted with any function that performs interleaving, I did create my own function. We load this script into GHCi. Here is a brief GHCi session:

*Main> -- create deck
*Main> let deck = [Ace,Two,Three,Four,Five,Six,Seven,Eight]
*Main> -- perform 3 out shuffles
*Main> take 4 (iterate outShuffle deck)
[[Ace,Two,Three,Four,Five,Six,Seven,Eight],
 [Ace,Five,Two,Six,Three,Seven,Four,Eight],
 [Ace,Three,Five,Seven,Two,Four,Six,Eight],
 [Ace,Two,Three,Four,Five,Six,Seven,Eight]]
*Main>

Indeed, after three perfect out shuffles, we do obtain a deck in the original order! One is tempted to wonder what the order of the deck will be after three perfect in shuffles. Let us see:

*Main> -- create deck
*Main> let deck = [Ace,Two,Three,Four,Five,Six,Seven,Eight]
*Main> -- perform 3 in shuffles
*Main> take 4 (iterate inShuffle deck)
[[Ace,Two,Three,Four,Five,Six,Seven,Eight],
 [Five,Ace,Six,Two,Seven,Three,Eight,Four],
 [Seven,Five,Three,Ace,Eight,Six,Four,Two],
 [Eight,Seven,Six,Five,Four,Three,Two,Ace]]

After three perfect in shuffles we do not obtain the original deck. Instead, we obtain a reversal of the original deck.

__________

Example 2

Let us now work with the standard deck of 52 cards. People say that eight consecutive perfect out shuffles will restore the deck to its original order. Let us check whether this is true or not.

Here is a new Haskell script:

-- create Value data type
data Value = Ace | Two | Three | Four | Five | Six | Seven | Eight | 
             Nine | Ten | Jack | Queen | King deriving (Show, Eq)

-- create Suit data type
data Suit  = Club | Diamond | Heart | Spade deriving (Show, Eq)

-- create type synonyms
type Card = (Value,Suit)
type Deck = [Card]

-- create list of card values
values :: [Value]
values = [Ace,Two,Three,Four,Five,Six,Seven,
          Eight,Nine,Ten,Jack,Queen,King]

-- create list of card suits
suits :: [Suit]
suits = [Club,Diamond,Heart,Spade]

-- create deck
deck :: Deck
deck = [(v,s) | v <- values, s <- suits] 

-- interlace two decks of the same size
interlace :: Deck -> Deck -> Deck
interlace [] [] = []
interlace (c1:d1) (c2:d2) = [c1,c2] ++ interlace d1 d2

-- perfect out shuffle
outShuffle :: Deck -> Deck
outShuffle deck = interlace deck1 deck2
                  where n = div (length deck) 2
                        deck1 = fst (splitAt n deck)
                        deck2 = snd (splitAt n deck)

-- perfect in shuffle
inShuffle :: Deck -> Deck
inShuffle deck = interlace deck2 deck1
                 where n = div (length deck) 2
                       deck1 = fst (splitAt n deck)
                       deck2 = snd (splitAt n deck)

We load this script into GHCi. Here is a GHCi session:

*Main> -- print deck
*Main> deck
[(Ace,Club),(Ace,Diamond),(Ace,Heart),(Ace,Spade),
 (Two,Club),(Two,Diamond),(Two,Heart),(Two,Spade),
 (Three,Club),(Three,Diamond),(Three,Heart),(Three,Spade),
 (Four,Club),(Four,Diamond),(Four,Heart),(Four,Spade),
 (Five,Club),(Five,Diamond),(Five,Heart),(Five,Spade),
 (Six,Club),(Six,Diamond),(Six,Heart),(Six,Spade),
 (Seven,Club),(Seven,Diamond),(Seven,Heart),(Seven,Spade),
 (Eight,Club),(Eight,Diamond),(Eight,Heart),(Eight,Spade),
 (Nine,Club),(Nine,Diamond),(Nine,Heart),(Nine,Spade),
 (Ten,Club),(Ten,Diamond),(Ten,Heart),(Ten,Spade),
 (Jack,Club),(Jack,Diamond),(Jack,Heart),(Jack,Spade),
 (Queen,Club),(Queen,Diamond),(Queen,Heart),(Queen,Spade),
 (King,Club),(King,Diamond),(King,Heart),(King,Spade)]
*Main> -- perform 8 out shuffles
*Main> let deck8 = (iterate outShuffle deck) !! 8
*Main> -- is the new deck the same as the original one?
*Main> deck8 == deck
True

Indeed, it is true! Eight perfect out shuffles suffice!

__________

References

[1] Persi Diaconis, R. L. Graham, William M. Kantor, The Mathematics of Perfect Shuffles, Advances in Applied Mathematics, Volume 4, Issue 2, June 1983.

Grabenkampf

April 7, 2012

Ernst Jünger on the brutality of trench warfare:

Auch das moderne Gefecht hat seine großen Augenblicke. Man hört so oft die irrige Ansicht, daß der Infanteriekampf zu einer uninteressanten Massenschlächterei herabgesunken ist. Im Gegenteil, heute mehr denn je entscheidet der einzelne. Das weiß jeder, der sie in ihrem Reich gesehen hat, die Fürsten des Grabens mit den harten, entschlossenen Gesichtern, tollkühn, so sehnig, geschmeidig vor- und zurückspringend, mit scharfen, blutdürstigen Augen, Helden, die kein Bericht nennt. Der Grabenkampf ist der blutigste, wildeste, brutalste von allen, doch auch er hat seine Männer gehabt, Männer, die ihrer Stunde gewachsen waren, unbekannte, verwegene Kämpfer. Unter allen nervenerregenden Momenten des Krieges ist keiner so stark, wie die Begegnung zweier Stoßtruppführer zwischen den engen Lehmwänden des Grabens. Da gibt es kein Zurück und kein Erbarmen. Blut klingt aus dem schrillen Erkennungsschrei, der sich wie Alpdruck von der Brust ringt.

[ source ]

__________

Here is Basil Creighton‘s 1929 translation:

Even modern battle has its great moments. One hears it said very often and very mistakenly that the infantry battle has degenerated to an uninteresting butchery. On the contrary, to-day more than ever it is the individual that counts. Every one knows that who has seen them in their own realm, these princes of the trenches, with their hard, set faces, brave to madness, tough and agile to leap forward or back, with keen bloodthirsty nerves, whom no despatch ever mentions. Trench warfare is the bloodiest, wildest, and most brutal of all warfare, yet it too has had its men, men whom the call of the hour has raised up, unknown foolhardy fighters. Of all the nerve-racking moments of war none is so formidable as the meeting of two storm-troop leaders between the narrow walls of the trench. There is no retreat and no mercy then. Blood sounds in the shrill cry that is wrung like a nightmare from the breast.

__________

Source:

Ernst Jünger, In Stahlgewittern, Berlin 1922.

Digraphs in Haskell

April 3, 2012

A directed graph (or digraph) is an ordered pair D = (V, A), where V is a non-empty finite set of elements called vertices, and A is a finite set of elements called arcs, where an arc is an ordered pair of distinct vertices [1]. We call V and A the vertex set and the arc set, respectively. Note that A \subseteq V \times V. For an arc (v_i, v_j), we call v_i the tail, and v_j the head, and say that the arc (v_i, v_j) leaves vertex v_i and enters vertex v_j.

For instance, consider the digraph depicted below

This digraph has 10 vertices and 15 arcs. How would we build this digraph in Haskell? A vertex could be represented by an integer, and an arc could be represented by an ordered pair of integers. However, in Mathematics we work with sets, whereas in Haskell we work with lists. Hence, we could represent the vertex set by a list of vertices, and the arc set by a list of edges.

Here is a simplistic implementation:

-- type synonyms
type Vertex  = Int
type Arc     = (Vertex,Vertex)
type Digraph = ([Vertex],[Arc])

-- create list of vertices
vs :: [Vertex]
vs = [1..10]

-- create list of arcs
as :: [Arc]
as = [(1,4),(1,5),(1,6),(2,4),(3,4),
      (4,5),(5,6),(7,6),(8,1),(8,2),
      (8,7),(9,8),(9,10),(10,2),(10,3)]

-- create digraph
digraph :: Digraph
digraph = (vs,as)

We load this script in GHCi and perform some simple testing:

*Main> -- list of vertices
*Main> vs
[1,2,3,4,5,6,7,8,9,10]
*Main> length vs
10
*Main> -- list of arcs
*Main> as
[(1,4),(1,5),(1,6),(2,4),(3,4),(4,5),(5,6),(7,6),
(8,1),(8,2),(8,7),(9,8),(9,10),(10,2),(10,3)]
*Main> length as
15
*Main> -- digraph
*Main> digraph
([1,2,3,4,5,6,7,8,9,10],[(1,4),(1,5),(1,6),
(2,4),(3,4),(4,5),(5,6),(7,6),(8,1),(8,2),
(8,7),(9,8),(9,10),(10,2),(10,3)])

We have thus created a digraph. What can we do with it? Is there anything of interest that we could compute?

__________

Neighborhoods

Given a digraph D = (V, A), we call the following sets

N_D^{+} (v) = \{ u \in V \mid (v, u) \in A\}

N_D^{-} (v) = \{ u \in V \mid (u, v) \in A\}

the out-neighborhood of vertex v and the in-neighborhood of vertex v, respectively [1]. The vertices in N_D^{+} (v) and N_D^{-} (v) are called the out-neighbors and the in-neighbors of vertex v, respectively.

In Haskell, we compute the out- and in-neighborhoods as follows:

-- compute out-neighborbood of a vertex
outNeighbors :: Digraph -> Vertex -> [Vertex]
outNeighbors (vs,as) v = map snd (filter (p v) as)
                         where p v (t,h) | t == v = True
                                         | t /= v = False

-- compute in-neighborbood of a vertex
inNeighbors :: Digraph -> Vertex -> [Vertex]
inNeighbors (vs,as) v = map fst (filter (p v) as)
                        where p v (t,h) | h == v = True
                                        | h /= v = False

Let us do some testing in GHCi:

*Main> -- list of all out-neighborhoods
*Main> map (outNeighbors digraph) vs
[[4,5,6],[4],[4],[5],[6],[],[6],[1,2,7],[8,10],[2,3]]
*Main> -- list of all in-neighborhoods
*Main> map (inNeighbors digraph) vs
[[8],[8,10],[10],[1,2,3],[1,4],[1,5,7],[8],[9],[],[9]]

 __________

Degrees

Given a digraph D = (V, A), we call the following cardinalities

d_D^{+} (v) = |N_D^{+} (v)|

d_D^{-} (v) = |N_D^{-} (v)|

the out-degree and the in-degree of vertex v, respectively [1]. We compute the out- and in-degrees in Haskell as follows:

-- compute out-degree of a vertex
outDegree :: Digraph -> Vertex -> Int
outDegree (vs,as) v = length (outNeighbors (vs,as) v)

-- compute in-degree of a vertex
inDegree :: Digraph -> Vertex -> Int
inDegree (vs,as) v = length (inNeighbors (vs,as) v)

Let us do some testing in GHCi:

*Main> -- out-degrees of all vertices
*Main> map (outDegree digraph) vs
[3,1,1,1,1,0,1,3,2,2]
*Main> -- in-degrees of all vertices
*Main> map (inDegree digraph) vs
[1,2,1,3,2,3,1,1,0,1]

__________

Concluding remarks

Here is the whole Haskell script:

-- type synonyms
type Vertex  = Int
type Arc     = (Vertex,Vertex)
type Digraph = ([Vertex],[Arc])

-- create list of vertices
vs :: [Vertex]
vs = [1..10]

-- create list of arcs
as :: [Arc]
as = [(1,4),(1,5),(1,6),(2,4),(3,4),
      (4,5),(5,6),(7,6),(8,1),(8,2),
      (8,7),(9,8),(9,10),(10,2),(10,3)]

-- create digraph
digraph :: Digraph
digraph = (vs,as)

-- compute out-neighborbood of a vertex
outNeighbors :: Digraph -> Vertex -> [Vertex]
outNeighbors (vs,as) v = map snd (filter (p v) as)
                         where p v (t,h) | t == v = True
                                         | t /= v = False

-- compute in-neighborbood of a vertex
inNeighbors :: Digraph -> Vertex -> [Vertex]
inNeighbors (vs,as) v = map fst (filter (p v) as)
                        where p v (t,h) | h == v = True
                                        | h /= v = False

-- compute out-degree of a vertex
outDegree :: Digraph -> Vertex -> Int
outDegree (vs,as) v = length (outNeighbors (vs,as) v)

-- compute in-degree of a vertex
inDegree :: Digraph -> Vertex -> Int
inDegree (vs,as) v = length (inNeighbors (vs,as) v)

We now know how to create digraphs in Haskell and compute simple things, like neighborhoods and degrees. What else would it be of interest to compute?

__________

References

[1] Jørgen Bang-Jensen, Gregory Z. Gutin, Digraphs: Theory, Algorithms and Applications (2nd edition), Springer, 2009.

Newton’s method using rational arithmetic

April 2, 2012

The floating-point number system puts further restraints on what can be expected. We can at best hope to keep the relative error under control, and we cannot expect to find zeros far from the origin with great absolute accuracy.

Richard W. Hamming [1]

I am somewhat disenchanted with floating-point arithmetic. We all have implemented numerical methods using floating-point numbers. But, how many of us have implemented numerical methods using rational numbers of arbitrary precision? I certainly have not. Hence, let us give it a try. In this post I will implement Newton’s method [1] (also known as the Newton-Raphson method) in Haskell using arbitrary-precision rational numbers.

Given a continuous (real) function f : \mathbb{R} \to \mathbb{R} and x_0, an initial approximation of the zero of function f, we then obtain the next approximation using the following recurrence relation

x_{k+1} = x_k - \displaystyle\frac{f (x_k)}{f' (x_k)}

where f' is the first derivative of f. We will use arbitrary-precision rational numbers (represented as a ratio of two Integer values)

type Rational = Ratio Integer

which are defined in the Data.Ratio library.

__________

Example

Suppose we would like to find a rational approximation of \sqrt{2}. We create a real function defined by f (x) = x^2 - 2, whose zeros are  \sqrt{2} and -\sqrt{2}. Note that the first derivative of f is f' (x) = 2 x. The initial approximation is, say, x_0 = 1 / 1, which is closer to \sqrt{2} than to -\sqrt{2}. Hopefully, the sequence \{x_k\} will converge to \sqrt{2} instead of converging to -\sqrt{2}. We are performing a numerical experiment, and success is not guaranteed!

The following Haskell script implements Newton’s method:

import Data.Ratio

-- define function f
f :: Rational -> Rational
f x = x*x - 2

-- define 1st derivative of function f
f' :: Rational -> Rational
f' x = 2*x

-- define Newton map
g :: Rational -> Rational
g x = x - (f x / f' x)

-- initial approximation
x0 :: Rational
x0 = 1 % 1

-- list of rational approximations
xs :: [Rational]
xs = iterate g x0

where we used (higher-order) function iterate to create the list of rational approximations. We run the script above in GHCi and use take to compute the first few rational approximations in the list:

*Main> take 6 xs
[1 % 1,3 % 2,17 % 12,577 % 408,665857 % 470832,
886731088897 % 627013566048]

Hence, we have, for instance, the following rational approximation

x_5 = \displaystyle\frac{886731088897}{627013566048} \approx \sqrt{2}.

How good is this rational approximation? Let us see:

*Main> let x5 = xs !! 5
*Main> x5*x5 - 2
1 % 393146012008229658338304
*Main> 886731088897 / 627013566048
1.4142135623730951
*Main> sqrt 2
1.4142135623730951

Fairly good indeed! And all it took was five iterations!

__________

References

[1] Richard W. Hamming, Numerical Methods for Scientists and Engineers, Dover Publications, 1973.


Follow

Get every new post delivered to your Inbox.

Join 76 other followers