## Archive for April, 2012

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]

-- 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
*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.

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.

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.