sage « Yet Another Mathblog

Sestinas and Sage

According to [1], a sestina is a highly structured poem consisting of six six-line stanzas followed by a tercet (called its envoy or tornada), for a total of thirty-nine lines. The same set of six words ends the lines of each of the six-line stanzas, but in a shuffled order each time. The shuffle used is very similar to the Mongean shuffle.

Define $f_n(k) = 2k$ , if k <= n/2 and $f_n(k) = 2n+1-2k$ , if $k > n/2.$ Let $p = (p_1,...,p_n) \in S_n$ , where $p_j = f_n(p_{j-1})$ and $S_n$ is the symmetric group of order $n$ . From [2], we have the following result.

Theorem: If p is an n-cycle then 2n+1 is a prime.

Call such a prime a “sestina prime”. Which primes are sestina primes?

Here is Python/Sage code for this permutation:


def sestina(n):
    """
    Computes the element of the symmetric group S_n associated to the shuffle above. 
  
    EXAMPLES: 
        sage: sestina(4)
        (1,2,4) 
        sage: sestina(6)
        (1,2,4,5,3,6)
        sage: sestina(8)
        (1,2,4,8)(3,6,5,7)
        sage: sestina(10)
        (1,2,4,8,5,10)(3,6,9)
        sage: sestina(12)
        (1,2,4,8,9,7,11,3,6,12)(5,10)
        sage: sestina(14)
        (1,2,4,8,13,3,6,12,5,10,9,11,7,14) 
        sage: sestina(16)
        (1,2,4,8,16)(3,6,12,9,15)(5,10,13,7,14)
        sage: sestina(18)
        (1,2,4,8,16,5,10,17,3,6,12,13,11,15,7,14,9,18)
        sage: sestina(20) (1,2,4,8,16,9,18,5,10,20)(3,6,12,17,7,14,13,15,11,19) 
        sage: sestina(22) (1,2,4,8,16,13,19,7,14,17,11,22)(3,6,12,21)(5,10,20)(9,18) 

    """ 
    def fcn(k, n):
        if k<=int(n/2): 
            return 2*k 
        else: 
            return 2*n+1-2*k 
    L = [fcn(k,n) for k in range(1,n+1)] 
    G = SymmetricGroup(n) 
    return G(L)

And here is an example due to Ezra Pound [3]:

                                  I

Damn it all! all this our South stinks peace.
You whoreson dog, Papiols, come! Let’s to music!
I have no life save when the swords clash.
But ah! when I see the standards gold, vair, purple, opposing
And the broad fields beneath them turn crimson,
Then howl I my heart nigh mad with rejoicing.

                                     II

In hot summer have I great rejoicing
When the tempests kill the earth’s foul peace,
And the light’nings from black heav’n flash crimson,
And the fierce thunders roar me their music
And the winds shriek through the clouds mad, opposing,
And through all the riven skies God’s swords clash.

                                     III

Hell grant soon we hear again the swords clash!
And the shrill neighs of destriers in battle rejoicing,
Spiked breast to spiked breast opposing!
Better one hour’s stour than a year’s peace
With fat boards, bawds, wine and frail music!
Bah! there’s no wine like the blood’s crimson!

                                     IV

And I love to see the sun rise blood-crimson.
And I watch his spears through the dark clash
And it fills all my heart with rejoicing
And prys wide my mouth with fast music
When I see him so scorn and defy peace,
His lone might ’gainst all darkness opposing.

                                     V

The man who fears war and squats opposing
My words for stour, hath no blood of crimson
But is fit only to rot in womanish peace
Far from where worth’s won and the swords clash
For the death of such sluts I go rejoicing;
Yea, I fill all the air with my music.

                                     VI

Papiols, Papiols, to the music!
There’s no sound like to swords swords opposing,
No cry like the battle’s rejoicing
When our elbows and swords drip the crimson
And our charges ’gainst “The Leopard’s” rush clash.
May God damn for ever all who cry “Peace!”

                                     VII

And let the music of the swords make them crimson
Hell grant soon we hear again the swords clash!
Hell blot black for always the thought “Peace”!

References:

[1] http://en.wikipedia.org/wiki/Sestina

[2] Richard Dore and Anton Geraschenko,”Sestinas and Primes” posted to http://stacky.net/wiki/index.php?title=Course_notes, and http://math.berkeley.edu/~anton/written/sestina.pdf

[3] Ezra Pound, “Sestina: Altaforte” (1909), (originally published int the English Review, 1909)

[4] John Bullitt, N. J. A. Sloane and J. H. Conway , http://oeis.org/A019567

2012/04/02 Posted by wdjoyner | math, sage | Mongean shuffle, poetry, prime numbers, Python, sage, sestina, symmetric group | Leave a Comment

Remarks on the Hamming [7.4.3] code and Sage

This post simply collects some very well-known facts and observations in one place, since I was having a hard time locating a convenient reference.

Let $C$ be the binary Hamming [7,4,3] code defined by the generator matrix $G = \left(\begin{array}{rrrrrrr} 1 & 0 & 0 & 0 & 1 & 1 & 1 \\ 0 & 1 & 0 & 0 & 1 & 1 & 0 \\ 0 & 0 & 1 & 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 & 1 & 1 \end{array}\right)$ and check matrix $H = \left(\begin{array}{rrrrrrr} 1 & 1 & 1 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 1 \end{array}\right)$ . In other words, this code is the row space of G and the kernel of H. We can enter these into Sage as follows:

sage: G = matrix(GF(2), [[1,0,0,0,1,1,1],[0,1,0,0,1,1,0],[0,0,1,0,1,0,1],[0,0,0,1,0,1,1]])
sage: G
[1 0 0 0 1 1 1]
[0 1 0 0 1 1 0]
[0 0 1 0 1 0 1]
[0 0 0 1 0 1 1]
sage: H = matrix(GF(2), [[1,1,1,0,1,0,0],[1,1,0,1,0,1,0],[1,0,1,1,0,0,1]])
sage: H
[1 1 1 0 1 0 0]
[1 1 0 1 0 1 0]
[1 0 1 1 0 0 1]
sage: C = LinearCode(G)
sage: C
Linear code of length 7, dimension 4 over Finite Field of size 2
sage: C = LinearCodeFromCheckMatrix(H)
sage: LinearCode(G) == LinearCodeFromCheckMatrix(H)
True

The generator matrix gives us a one-to-one onto map $G:GF(2)^4\to C$ defined by $m \longmapsto m\cdot G$ . Using this map, the codewords are easy to describe and enumerate:

$\begin{tabular}{ccc} decimal & binary & codeword \\ 0 & 0000 & 0000000 \\ 1 & 0001 & 0001011 \\ 2 & 0010 & 0010101 \\ 3 & 0011 & 0011110 \\ 4 & 0100 & 0100110 \\ 5 & 0101 & 0101101 \\ 6 & 0110 & 0110011 \\ 7 & 0111 & 0111000 \\ 8 & 1000 & 1000111 \\ 9 & 1001 & 1001100 \\ 10 & 1010 & 1010010 \\ 11 & 1011 & 1011001 \\ 12 & 1100 & 1100001 \\ 13 & 1101 & 1101010 \\ 14 & 1110 & 1110100 \\ 15 & 1111 & 1111111 \end{tabular}$ .

Using this code, we first describe a guessing game you can play with even small children.

Number Guessing game: Pick an integer from 0 to 15. I will ask you 7 yes/no questions. You may lie once.
I will tell you when you lied and what the correct number is.

Question 1: Is n in {0,1,2,3,4,5,6,7}?
(Translated: Is 1st bit of Hamming_code(n) a 0?)
Question 2: Is n in {0,1,2,3,8,9,10,11}?
(Is 2nd bit of Hamming_code(n) a 0?)
Question 3: Is n in {0,1,4,5,8,9,12,13}?
(Is 3rd bit of Hamming_code(n) a 0?)
Question 4: Is n in {0,2,4,6,8,10,12,14} (ie, is n even)?
(Is 4th bit of Hamming_code(n) a 0?)
Question 5: Is n in {0,1,6,7,10,11,12,13}?
(Is 5th bit of Hamming_code(n) a 0?)
Question 6: Is n in {0,2,5,7,9,11,12,14}?
(Is 6th bit of Hamming_code(n) a 0?)
Question 7: Is n in {0,3,4,7,9,10,13,14}?
(Is 7th bit of Hamming_code(n) a 0?)

Record the answers in a vector (0 for yes, 1 for no): $v = (v_1,v_2,...,v_7)$ . This must be a codeword (no lies) or differ from a codeword by exactly one bit (1 lie). In either case, you can find n by decoding this vector.

We discuss a few decoding algorithms next.

Venn diagram decoding:

We use a simple Venn diagram to describe a decoding method.

sage: t = var('t')
sage: circle1 = parametric_plot([10*cos(t)-5,10*sin(t)+5], (t,0,2*pi))
sage: circle2 = parametric_plot([10*cos(t)+5,10*sin(t)+5], (t,0,2*pi))
sage: circle3 = parametric_plot([10*cos(t),10*sin(t)-5], (t,0,2*pi))
sage: text1 = text("$1$", (0,0))  
sage: text2 = text("$2$", (-6,-2)) 
sage: text3 = text("$3$", (0,7))
sage: text4 = text("$4$", (6,-2)) 
sage: text5 = text("$5$", (-9,9)) 
sage: text6 = text("$6$", (9,9)) 
sage: text7 = text("$7$", (0,-9))  
sage: textA = text("$A$", (-13,13)) 
sage: textB = text("$B$", (13,13)) 
sage: textC = text("$C$", (0,-17)) 
sage: text_all = text1+text2+text3+text4+text5+text6+text7+textA+textB+textC
sage: show(circle1+circle2+circle3+text_all,axes=false)

This gives us the following diagram:

Decoding algorithm:
Suppose you receive $v = ( v_1, v_2, v_3, v_4, v_5, v_6, v_7)$ .
Assume at most one error is made.
Decoding process:

Place $v_i$ in region i of the Venn diagram.
For each of the circles A, B, C, determine if the sum of the bits in four regions add up to 0 or to 1. If they add to 1, say that that circle has a “parity failure”.
The error region is determined form the following table.

$\begin{tabular}{cc} parity failure region(s) & error position \\ none & none \\ A, B, and C & 1 \\ B, C & 4 \\ A, C & 2 \\ A, B & 3 \\ A & 5 \\ B & 6 \\ C & 7 \end{tabular}$

For example, suppose v = (1,1,1,1,1,0,1). The filled in diagram looks like

This only fails in circle B, so the table says (correctly) that the error is in the 6th bit. The decoded codeword is $c = v+e_6 = (1,1,1,1,1,1,1).$

Next, we discuss a decoding method based on the Tanner graph.

Tanner graph for hamming 7,4,3 code

The above Venn diagram corresponds to a bipartite graph, where the left “bit vertices” (1,2,3,4,5,6,7) correspond to the coordinates in the codeword and the right “check vertices” (8,9,10) correspond to the parity check equations as defined by the check matrix. This graph corresponds to the above Venn diagram, where the check vertices 8, 9, 10 were represented by circles A, B, C:

sage: Gamma = Graph({8:[1,2,3,5], 9:[1,2,4,6], 10:[1,3,4,7]})
sage: B = BipartiteGraph(Gamma)
sage: B.show()
sage: B.left
set([1, 2, 3, 4, 5, 6, 7])
sage: B.right
set([8, 9, 10])
sage: B.show()

This gives us the graph in the following picture:

Decoding algorithm:
Suppose you receive $v = ( v_1, v_2, v_3, v_4, v_5, v_6, v_7)$ .
Assume at most one error is made.
Decoding process:

Place $v_i$ at the vertex i on the left side of the bipartite graph.
For each of the check vertices 8,9,10 on the right side of the graph, determine of the if the sum of the bits in the four left-hand vertices connected to it add up to 0 or to 1. If they add to 1, we say that that check vertex has a “parity failure”.
Those check vertices which do not fail are connected to bit vertices which we assume are correct. The remaining bit vertices
connected to check vertices which fail are to be determined (if possible) by solving the corresponding check equations.

check vertex 8: $v_2+v_3+v_4+v_5 = 0$

check vertex 9: $v_1+v_3+v_4+v_6 = 0$

check vertex 10: $v_1+v_2+v_4+v_7 = 0$

Warning: This method is not guaranteed to succeed in general. However, it does work very efficiently when the check matrix H is “sparse” and the number of 1′s in each row and column is “small.”

For example, suppose v = (1,1,1,1,1,0,1). The check vertex 9 fails its parity check, but vertex 8 and 10 do not. Therefore, only bit vertex 6 is unknown, since vertex 6 is the only one not connected to 8 and 10. This tells us that the decoding codeword is $c = (1,1,1,1,1,v_6,1)$ , for some unknown $v_6$ . We solve for this unknown using the check vertex equation $v_1+v_3+v_4+v_6 = 0$ , giving us $v_6 = 1$ . The decoded codeword is $c = (1,1,1,1,1,1,1).$

This last example was pretty simple, so let’s try $v=(0,1,1,1,1,1,1)$ . In this case, we know the vertices 9 and 10 fail, so $c = (v_1,1,1,1,1,v_6,v_7)$ . We solve using

$v_1+1+1+v_6 = 0$

$v_1+1+1+v_7 = 0$

This simply tells us $v_1=v_6=v_7$ . By majority vote, we get $c = (1,1,1,1,1,1,1)$ .

2012/02/28 Posted by wdjoyner | math, sage | coding theory, graph theory, Hamming codes | Leave a Comment

Digital steganography and Sage

This post shall explain how the F5 stegosystem (in a simple case) can be implemented in Sage. I thank Carlos Munuera for teaching me these ideas and for many helpful conversations on this matter. I also thank Kyle Tucker-Davis [1] for many interesting conversations on this topic.

Steganography, meaning “covered writing,” is the science of secret communication [4]. The medium used to carry the information is called the “cover” or “stego-cover” (depending on the context – these are not synonymous terms). The term “digital steganography” refers to secret communication where the cover is a digital media file.

One of the most common systems of digital steganography is the Least Significant Bit (LSB) system. In this system, we assume the “cover” image is represented as a finite sequence of binary vectors of a given length. In other words, a greyscale image is regarded as an array of pixels, where each pixel is represented by a binary vector of a certain fixed length. Each such vector represents a pixel’s brightness in the cover image. The encoder embeds (at most) one bit of information in the least significant bit of eaach vector. Care must be taken with this system to ensure that the changes made do not betray the stego-cover, while still maximizing the information hidden.

From a short note of Crandell [3] in 1998, it was realized that error-correcting codes can give rise to “stego-schemes”, ie, methods by which a message can be hidden in a digital file efficiently.

Idea in a nutshell: If $C$ is the r-th binary Hamming code (so, $n = 2^r-1$ is its length, $k = 2^r - r - 1$ is its dimension, and $d = 3$ is its minimum distance), $G$ is a generating matrix, $H$ is a check matrix, and $C^\perp$ is the dual code of $C$ , then we take an element $v \in GF(2)^n$ to be a cover, and the message $m$ we embed is an element of $GF(2)^r$ . Once we find a vector $z$ of lowest weight such that $H(v+z)=m$ , we call $v+z$ the stegocover. The stegocover looks a lot like the original cover and “contains” the message m. This will be explained in more detai below.
(This particular scheme is called the F5 stegosystem, and is due to Westfeld.)

Quick background on error-correcting, linear, block codes [5]

A linear error-correcting block code is a finite dimensional vector space over a finite field with a fixed basis. We assume the finite field is the binary field $GF(2)$ .

We shall typically think of a such a code as a subspace $C$ of $GF(2)^n$ with a fixed basis, where $n>0$ is an integer called the length of the code. Moreover, the basis for the ambient space $GF(2)^n$ will be the standard basis,

$e_1=(1,0,\dots, 0), e_2=(0,1,0,\dots, 0), \dots, e_n=(0,\dots, 0,1).$

There are two common ways to specify a linear code $C$ .

You can give $C$ as a vector subspace of $GF(2)^n$ by specifying a set of basis vectors for $C$ . This set of basis vectors is, by convention, placed as
the rows of a matrix called a generator matrix $G$ of $C$ . Obviously, the order in which the rows are presented does not
affect the code itself.

If $g_1, \dots, g_k$ are the rows of $G$ then

$C= \{c=m_1g_1+\dots +m_kg_k\ |\ {\rm some}\ m_i\in GF(2)\},$

is the set of linear combinations of the row vectors $g_i$ . The vector of coefficients, $m=(m_1,\dots, m_k)$ represents the information you want to encode and transmit.

In other words, encoding of a message can be defined via the generator matrix:

$\begin{array}{ccc} m = (m_1,\dots, m_k) & \to & c=m_1g_1+\dots +m_kg_k = m^t\cdot G,\\ GF(2)^k & \to & C. \end{array}$
You can give $C$ as a vector subspace of $GF(2)^n$ by specifying a matrix $H$ for which $C$ is the kernel of $H$ , $C={\rm ker}(C)$ . This matrix is called a check matrix of $C$ . Again, the order in which the rows are presented does not affect the code itself.

Note that if $G$ is a full rank $k\times n$ matrix then a full rank check matrix $H$ must be a $(n-k) \times n$ matrix.

These two ways of defining a code are not unrelated.

Fact:
If $G=(I_k\ \vert\ A)$ is the generating matrix for $C$ then $H=(-A^t\ \vert\ I_{n-k})$ is a parity check matrix.

Exaample:
Let $r>0$ be an integer and let $H$ be a $r\times (2^r-1)$ matrix whose columns are all the distinct non-zero vectors of $GF(2)^r$ . Then, the code having $H$ as its check matrix is called a binary Hamming code, denoted Ham(r,2).

Let $r=3$ , and let

$H= \begin{pmatrix} 0 & 1 & 1 & 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 & 1 & 0 \\ 1 & 1 & 0 & 1 & 0 & 0 & 1 \end{pmatrix}.$

In this form, namely when the columns of $H$ are arranged in standard form such that the rightmost $k\times k$ entries is the identity matrix $I_k$ , the generator matrix $G$ can be quickly found to be

$G= \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 1 & 1 \\ 0 & 1 & 0 & 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 1 & 1 & 1 & 1 \end{pmatrix}.$

A coset of $GF(2)^n/C$ is a subset of $GF(2)^n$ of the form $C+v$ for some $v\in GF(2)^n$ . Let $S$ be a coset of $C$ . Then,
a coset leader of $S$ is an element of $S$ having smallest Hamming weight.

Fact:
The coset leaders of a Hamming code are those vectors of $wt \leq 1$ .

Steganographic systems from error-correcting codes

This section basically describes Crandell’s idea [3] in a more formalized language.

Following Munuera’s notation in [6], a steganographic system $S$ can be formally defined as

$S= \{\mathcal{C}, \mathcal{M}, \mathcal{K}, emb, rec \},$
where

$\mathcal{C}$ is a set of all possible covers
$\mathcal{M}$ is a set of all possible messages
$\mathcal{K}$ is a set of all possible keys
$emb:\mathcal{C}\times\mathcal{M}\times\mathcal{K}\to\mathcal{C}$ is an embedding function
$rec:\mathcal{C}\times\mathcal{K}\to\mathcal{M}$ is a recovery function

and these all satisfy

$rec(emb(c,m,k),k)=m,$

for all $m\in\mathcal{M}$ , $c\in\mathcal{C}$ , $k\in\mathcal{K}$ . We will assume that a fixed key $k\in\mathcal{K}$ is used, and therefore,
the dependence on an element in the keyspace $\mathcal{K}$ can be ignored. The original cover $c$ is called the plain cover, $m$ is called the message, and $emb(c,m,k)=emb(c,m)$ is called the stegocover. Let $\mathcal{C}$ be $GF(2)^n$ , representing both plain covers and stegocovers. Also, let $\mathcal{M}$ be $GF(2)^k$ , where $k$ is a fixed integer such that $0 \leq k \leq n$ .

Sage examples

How do we compute the emb map?

We need the following Sage function.

def hamming_code_coset_leader(C, y):
      """
      Finds the coset leader of a binary Hamming code C.
      EXAMPLES:
      """
      F = C.base_ring()
      n = C.length()
      k = C.dimension()
      r = n-k
      V0 = F^r
      if not(y in V0):
          RaiseError, "Input vector is not a syndrome."
      H = C.check_mat()
      colsH = H.columns()
      i = colsH.index(y)
      V = F^n
      v = V(0)
      v[i] = F(1)
      return v

Let $V = GF(2)^{63}$ and consider the cover $v = ([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 1,1,1,1,0,0,1,1,1, 1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,1,0,0,1,1,0,0,1,1)$ . Regarded as a $7\times 9$ matrix,

V = GF(2)^(63)
rhino = V([1,1,1,1,1,1,1,1,1,
           1,1,1,1,1,1,1,1,1,
           1,1,1,1,0,0,1,1,1,
           1,0,0,0,0,0,0,1,1,
           1,0,0,0,0,0,0,0,1,
           1,0,0,0,0,0,0,1,1,
           1,0,0,1,1,0,0,1,1])
A = matrix(GF(2),7,rhino.list())
matrix_plot(A)

this looks like an elephant or a rhino:

Now we embed the message $m = (1,0,1,0,1,0)$ . First we compute the stegocover:

C = HammingCode(6,GF(2))
H = C.check_mat()
V0 = GF(2)^6
m = V0([1,0,1,0,1,0])
z = hamming_code_coset_leader(C, m)
stegocover = rhino+z
A = matrix(GF(2),7,stegocover.list())
matrix_plot(A)

It looks like another rhino/elephant:

Note only one bit is changed since the Hamming weight of z is at most 1. To recover the message m, just multiply the vector “stegocover” by H.

That’s how you can use Sage to understand the F5 stegosystem, at least in a very simple case.

REFERENCES:
[1] Kyle Tucker-Davis, “An analysis of the F5 steganographic system”, Honors Project 2010-2011
http://www.usna.edu/Users/math/wdj/tucker-davis/

[2] J. Bierbrauer and J. Fridrich. “Constructing good covering codes for applications in Steganography}. 2006.
http://www.math.mtu.edu/jbierbra/.

[3] Crandall, Ron. “Some Notes on Steganography”. Posted on a Steganography Mailing List, 1998.
http://os.inf.tu-dresden.de/westfeld/crandall.pdf

[4] Wayner, Peter. Disappearing Cryptography. Morgan Kauffman Publishers. 2009.

[5] Hill, Raymond. A First Course in Coding Theory. Oxford University Press. 1997.

[6] Munuera, Carlos. “Steganography from a Coding Theory Point of View”. 2010.
http://www.singacom.uva.es/oldsite/Actividad/s3cm/s3cm10/10courses.html

[7] Zhang, Weiming and S. Li. “Steganographic Codes- a New Problem of Coding Theory”. preprint, 2005.
http://arxiv.org/abs/cs/0505072.

2012/02/22 Posted by wdjoyner | math, sage | error-correcting codes, F5 stegosystem, Hamming codes, steganography | Leave a Comment

Applications of graphs to Boolean functions

Let f be a Boolean function on $GF(2)^n$ . The Cayley graph of f is defined to be the graph

$\Gamma_f = (GF(2)^n, E_f )$ ,

whose vertex set is $GF(2)^n$ and the set of edges is defined by

$E_f =\{ (u,v) \in GF(2)^n \times GF(2)^n \ |\ f(u+v)=1\}$ .

The adjacency matrix $A_f$ is the matrix whose entries are

$A_{i,j} = f(b(i) + b(j))$ ,

where b(k) is the binary representation of the integer k.
Note $\Gamma_f$ is a regular graph of degree wt(f), where wt denotes the Hamming weight of f when regarded as a vector of values (of length $2^n$ ).

Recall that, given a graph $\Gamma$ and its adjacency matrix A, the spectrum Spec( $\Gamma$ ) is the multi-set of eigenvalues of A.

The Walsh transform of a Boolean function f is an integer-valued function over $GF(2)^n$ that can be defined as

$W_f(u) = \sum_{x in GF(2)^n} (-1)^{f(x)+ \langle u,x\rangle}.$
A Boolean function f is bent if $|W_f(a)| = 2^{n/2}$ (this only makes sense if n is even). The Hadamard transform of a integer-valued function f is an integer-valued function over $GF(2)^n$ that can be defined as

$H_f(u) = \sum_{x in GF(2)^n} f(x)(-1)^{\langle u,x\rangle}.$
It turns out that the spectrum of $\Gamma_f$ is equal to the Hadamard transform of f when regarded as a vector of (integer) 0,1-values. (This nice fact seems to have first appeared in [2], [3].)

A graph is regular of degree r (or r-regular) if every vertex has degree r (number of edges incident to it).   We say that an r-regular graph $\Gamma$ is a strongly regular graph with parameters (v, r, d, e) (for nonnegative integers e, d) provided, for all vertices u, v the number of vertices adjacent to both u, v is equal to

          e, if u, v are adjacent,
          d, if u, v are nonadjacent.

It turns out tht f is bent iff $\Gamma_f$ is strongly regular and e = d (see [3] and [4]).

The following Sage computations illustrate these and other theorems in [1], [2], [3], [4].

Consider the Boolean function $f: GF(2)^4 \to GF(2)$ given by $f(x_0,x_1,x_2) = x_0x_1+x_2x_3$ .

sage: V = GF(2)^4
sage: f = lambda x: x[0]*x[1]+x[2]*x[3]
sage: CartesianProduct(range(16), range(16))
Cartesian product of [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], 
                     [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]
sage: C = CartesianProduct(range(16), range(16))
sage: Vlist = V.list()                  
sage: E = [(x[0],x[1]) for x in C if f(Vlist[x[0]]+Vlist[x[1]])==1]
sage: len(E)
96
sage: E = Set([Set(s) for s in E])
sage: E = [tuple(s) for s in E] 
sage: Gamma = Graph(E)
sage: Gamma
Graph on 16 vertices
sage: VG = Gamma.vertices()
sage: L1 = []
sage: L2 = []
sage: for v1 in VG:
....:         for v2 in VG:
....:             N1 = Gamma.neighbors(v1)
....:         N2 = Gamma.neighbors(v2)
....:         if v1 in N2:
....:                 L1 = L1+[len([x for x in N1 if x in N2])]
....:         if not(v1 in N2) and v1!=v2:
....:                 L2 = L2+[len([x for x in N1 if x in N2])]
....: 
....: 
sage: L1; L2
[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
 2, 2, 2, 2]

This implies the graph is strongly regular with d=e=2.

sage: Gamma.spectrum()
[6, 2, 2, 2, 2, 2, 2, -2, -2, -2, -2, -2, -2, -2, -2, -2]
sage: [walsh_transform(f, a) for a in V]
[4, 4, 4, -4, 4, 4, 4, -4, 4, 4, 4, -4, -4, -4, -4, 4]
sage: Omega_f = [v for v in V if f(v)==1] 
sage: len(Omega_f)
6
sage: Gamma.is_bipartite()
False
sage: Gamma.is_hamiltonian()
True
sage: Gamma.is_planar()     
False
sage: Gamma.is_regular()
True
sage: Gamma.is_eulerian()
True
sage: Gamma.is_connected()
True
sage: Gamma.is_triangle_free()
False
sage: Gamma.diameter()
2
sage: Gamma.degree_sequence()
[6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6]
sage: show(Gamma)
# bent-fcns-cayley-graphs1.png

Here is the picture of the graph:

sage: H = matrix(QQ, 16, 16, [(-1)^(Vlist[x[0]]).dot_product(Vlist[x[1]]) for x in C])
sage: H
[ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1]
[ 1 -1  1 -1  1 -1  1 -1  1 -1  1 -1  1 -1  1 -1]
[ 1  1 -1 -1  1  1 -1 -1  1  1 -1 -1  1  1 -1 -1]
[ 1 -1 -1  1  1 -1 -1  1  1 -1 -1  1  1 -1 -1  1]
[ 1  1  1  1 -1 -1 -1 -1  1  1  1  1 -1 -1 -1 -1]
[ 1 -1  1 -1 -1  1 -1  1  1 -1  1 -1 -1  1 -1  1]
[ 1  1 -1 -1 -1 -1  1  1  1  1 -1 -1 -1 -1  1  1]
[ 1 -1 -1  1 -1  1  1 -1  1 -1 -1  1 -1  1  1 -1]
[ 1  1  1  1  1  1  1  1 -1 -1 -1 -1 -1 -1 -1 -1]
[ 1 -1  1 -1  1 -1  1 -1 -1  1 -1  1 -1  1 -1  1]
[ 1  1 -1 -1  1  1 -1 -1 -1 -1  1  1 -1 -1  1  1]
[ 1 -1 -1  1  1 -1 -1  1 -1  1  1 -1 -1  1  1 -1]
[ 1  1  1  1 -1 -1 -1 -1 -1 -1 -1 -1  1  1  1  1]
[ 1 -1  1 -1 -1  1 -1  1 -1  1 -1  1  1 -1  1 -1]
[ 1  1 -1 -1 -1 -1  1  1 -1 -1  1  1  1  1 -1 -1]
[ 1 -1 -1  1 -1  1  1 -1 -1  1  1 -1  1 -1 -1  1]
sage: flist = vector(QQ, [int(f(v)) for v in V])
sage: H*flist  
(6, -2, -2, 2, -2, -2, -2, 2, -2, -2, -2, 2, 2, 2, 2, -2)
sage: A = matrix(QQ, 16, 16, [f(Vlist[x[0]]+Vlist[x[1]]) for x in C])
sage: A.eigenvalues()
[6, 2, 2, 2, 2, 2, 2, -2, -2, -2, -2, -2, -2, -2, -2, -2]

Here is another example: $f: GF(2)^3 \to GF(2)$ given by $f(x_0,x_1,x_2) = x_0x_1+x_2$ .

sage: V = GF(2)^3
sage: f = lambda x: x[0]*x[1]+x[2]
sage: Omega_f = [v for v in V if f(v)==1] 
sage: len(Omega_f)
4
sage: C = CartesianProduct(range(8), range(8))
sage: Vlist = V.list()    
sage: E = [(x[0],x[1]) for x in C if f(Vlist[x[0]]+Vlist[x[1]])==1]
sage: E = Set([Set(s) for s in E])
sage: E = [tuple(s) for s in E] 
sage: Gamma = Graph(E)
sage: Gamma
Graph on 8 vertices
sage: 
sage: VG = Gamma.vertices()
sage: L1 = []
sage: L2 = []
sage: for v1 in VG:
....:         for v2 in VG:
....:             N1 = Gamma.neighbors(v1)
....:         N2 = Gamma.neighbors(v2)
....:         if v1 in N2:
....:                 L1 = L1+[len([x for x in N1 if x in N2])]
....:         if not(v1 in N2) and v1!=v2:
....:                 L2 = L2+[len([x for x in N1 if x in N2])]
....: 
sage: L1; L2
[2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 0, 2, 2, 2, 2, 0, 0, 2, 2, 2, 2, 0, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2]
[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]

This implies that the graph is not strongly regular, therefore f is not bent.

sage: Gamma.spectrum()
[4, 2, 0, 0, 0, -2, -2, -2]
sage: 
sage: Gamma.is_bipartite()
False
sage: Gamma.is_hamiltonian()
True
sage: Gamma.is_planar()     
False
sage: Gamma.is_regular()
True
sage: Gamma.is_eulerian()
True
sage: Gamma.is_connected()
True
sage: Gamma.is_triangle_free()
False
sage: Gamma.diameter()
2
sage: Gamma.degree_sequence()
[4, 4, 4, 4, 4, 4, 4, 4]
sage: H = matrix(QQ, 8, 8, [(-1)^(Vlist[x[0]]).dot_product(Vlist[x[1]]) for x in C])
sage: H
[ 1  1  1  1  1  1  1  1]
[ 1 -1  1 -1  1 -1  1 -1]
[ 1  1 -1 -1  1  1 -1 -1]
[ 1 -1 -1  1  1 -1 -1  1]
[ 1  1  1  1 -1 -1 -1 -1]
[ 1 -1  1 -1 -1  1 -1  1]
[ 1  1 -1 -1 -1 -1  1  1]
[ 1 -1 -1  1 -1  1  1 -1]
sage: flist = vector(QQ, [int(f(v)) for v in V])
sage: H*flist  
(4, 0, 0, 0, -2, -2, -2, 2)
sage: Gamma.spectrum()
[4, 2, 0, 0, 0, -2, -2, -2]
sage: A = matrix(QQ, 8, 8, [f(Vlist[x[0]]+Vlist[x[1]]) for x in C])
sage: A.eigenvalues()
[4, 2, 0, 0, 0, -2, -2, -2]

sage: show(Gamma)
# bent-fcns-cayley-graphs2.png

Here is the picture:

REFERENCES:
[1] Pantelimon Stanica, Graph eigenvalues and Walsh spectrum of Boolean functions, INTEGERS 7(2) (2007), #A32.

[2] Anna Bernasconi, Mathematical techniques for the analysis of Boolean functions, Ph. D. dissertation TD-2/98, Universit di Pisa-Udine, 1998.

[3] Anna Bernasconi and Bruno Codenotti, Spectral Analysis of Boolean Functions as a Graph Eigenvalue Problem, IEEE TRANSACTIONS ON COMPUTERS, VOL. 48, NO. 3, MARCH 1999.

[4] A. Bernasconi, B. Codenotti, J.M. VanderKam. A Characterization of Bent Functions in terms of Strongly Regular Graphs, IEEE Transactions on Computers, 50:9 (2001), 984-985.

[5] Michel Mitton, Minimal polynomial of Cayley graph adjacency matrix for Boolean functions, preprint, 2007.

[6] ——, On the Walsh-Fourier analysis of Boolean functions, preprint, 2006.

2012/02/16 Posted by wdjoyner | math, sage | adjacency matrix, boolean function, Boolean functions, cayley graph, graph theory, regular graph | Leave a Comment

The Vigenère cipher and Sage

The Vigenère cipher is named after Blaise de Vigenère, a sixteenth century diplomat and cryptographer, by a historical accident. Vigene`re actually invented a different and more complicated cipher. The so-called “Vigenère cipher” cipher was actually invented by Giovan Batista Belaso in 1553. In any case, it is this cipher which we shall discuss here first.

This cipher has been re-invented by several authors, such as author and mathematician Charles Lutwidge Dodgson (Lewis Carroll) who claimed his 1868 “The Alphabet Cipher” was unbreakable. Several others claimed the so-called Vigenère cipher was unbreakable (e.g., the Scientific American magazine in 1917). However, Friedrich Kasiski and Charles Babbage broke the cipher in the 1800′s [1]. This cipher was used in the 1700′s, for example, during the American Civil War. The Confederacy used a brass cipher disk to implement the Vigenère cipher (now on display in the NSA Museum in Fort Meade) [1].

The so-called Vigenère cipher is a generalization of the Cesaer shift cipher. Whereas the shift cipher shifts each letter by the same amount (that amount being the key of the shift cipher) the so-called Vigenère cipher shifts a letter by an amount determined by the key, which is a word or phrase known only to the sender and receiver).

For example, if the key was a single letter, such as “C”, then the so-called Vigenère cipher is actually a shift cipher with a shift of 2 (since “C” is the 2-nd letter of the alphabet, if you start counting at 0). If the key was a word with two letters, such as “CA”, then the so-called Vigenère cipher will shift letters in even positions by 2 and letters in odd positions are left alone (or shifted by 0, since “A” is the 0-th letter, if you start counting at 0).

REFERENCES:
[1] Wikipedia article on the Vigenere cipher.

Let’s look at a message (a chant at football games between rivals USNA and West Point):

sage: AS = AlphabeticStrings()           
sage: A = AS.alphabet()
sage: VC = VigenereCryptosystem(AS, 1) # sets the key length to be = 1
sage: m = VC.encoding("Beat Army!"); m  # trivial example
BEATARMY

Now, we choose for illustration a simple key of length 1, and encipher this message:

sage: key = AS("D")
sage: c = VC.enciphering(key, m)
sage: c  # a trivial example
EHDWDUPB

You see here that in this case the cipher boils down to the Caesar/shift cipher (shifting by 3).

Deciphering is easy:

sage: VC.deciphering(key, c)
BEATARMY

Next, we choose for illustration a simple key of length 2, and encipher the same message:

sage: VC = VigenereCryptosystem(AS, 2)
sage: key = AS("CA")
sage: m = VC.encoding("Beat Army!"); m
BEATARMY
sage: c = VC.enciphering(key, m); c
DECTCROY

Since one of the key letters is “A” (which shifts by 0), half the plaintext is unchanged in going to the ciphertext.

Here is the algorithmic description of the above (so-called) Vigenère cipher:

    ALGORITHM:
    INPUT: 
      key - a string of upper-case letters (the secret "key")
      m - string of upper-case letters (the "plaintext" message)
    OUTPUT:
      c - string of upper-case letters (the "ciphertext" message)

  Identify the alphabet A, ..., Z with the integers 0, ..., 25. 
    Step 1: Compute from the string key a list L1 of corresponding
            integers. Let n1 = len(L1).
    Step 2: Compute from the string m a list L2 of corresponding
            integers. Let n2 = len(L2).
    Step 3: Break L2 up sequencially into sublists of size n1, and one sublist
            at the end of size <=n1. 
    Step 4: For each of these sublists L of L2, compute a new list C given by 
            C[i] = L[i]+L1[i] (mod 26) to the i-th element in the sublist, 
            for each i.
    Step 5: Assemble these lists C by concatenation into a new list of length n2.
    Step 6: Compute from the new list a string c of corresponding letters.

Once it is known that the key is, say, n characters long, frequency analysis can be applied to every n-th letter of the ciphertext to determine the plaintext. This method is called “Kasiski examination”, or the “Kasiski attack” (although it was first discovered by Charles Babbage).

The cipher Vigenère actually discovered is an “auto-key cipher” described as follows.

ALGORITHM:
    INPUT: 
      key - a string of upper-case letters (the secret "key")
      m - string of upper-case letters (the "plaintext" message)
    OUTPUT:
      c - string of upper-case letters (the "ciphertext" message)

  Identify the alphabet A, ..., Z with the integers 0, ..., 25. 
    Step 1: Compute from the string m a list L2 of corresponding
            integers. Let n2 = len(L2). 
    Step 2: Let n1 be the length of the key. Concatenate the string 
            key with the first n2-n1 characters of the plaintext message.
            Compute from this string of length n2 a list L1 of corresponding
            integers. Note n2 = len(L1).
    Step 3: Compute a new list C given by C[i] = L1[i]+L2[i] (mod 26), for each i.
            Note n2 = len(C).
    Step 5: Compute from the new list a string c of corresponding letters.

Note how the key is mixed with the plaintext to create the enciphering of the plaintext to ciphertext in Steps 2 and 3.

A screencast describing this has been posted to vimeo.

2012/01/20 Posted by wdjoyner | math, sage | cipher, math, sage | Leave a Comment

Bifid cipher and Sage

The Bifid cipher was invented around 1901 by Felix Delastelle. It is a “fractional substitution” cipher, where letters are replaced by pairs of symbols from a smaller alphabet. The cipher uses a 5×5 square filled with some ordering of the alphabet, except that “i”‘s and “j”‘s are identified (this is a so-called Polybius square; there is a 6×6 analog if you add back in “j” and also append onto the usual 26 letter alphabet, the digits 0, 1, …, 9). According to Helen Gaines’ book “Cryptanalysis”, this type of cipher was used in the field by the German army during World War I.

The Bifid cipher was discusses in Alasdair McAndrew’s book on Cryptography and Sage. We shall follow his discussion. As an example of a Polybius square for the Bifid cipher, pick the key to be “encrypt” (as Alasdair does). In that case, the Polybius square is $\left(\begin{array}{rrrrr} E & N & C & R & Y \\ P & T & A & B & C \\ D & E & F & G & H \\ I & K & L & M & N \\ O & P & Q & R & S \end{array}\right).$ BTW, the $6\times 6$ analog is: $\left(\begin{array}{rrrrrr} E & N & C & R & Y & P \\ T & A & B & C & D & E \\ F & G & H & I & J & K \\ L & M & N & O & P & Q \\ R & S & T & U & V & W \\ X & Y & Z & 0 & 1 & 2 \end{array}\right)$ .

Here is Sage code to produce the $6\times 6$ case (the $5\times 5$ case is in Alasdair’s book):

def bifid(pt, key):
    """
    INPUT:
        pt - plaintext string      (digits okay)
        key - short string for key (no repetitions, digits okay)
    
    OUTPUT:
        ciphertext from Bifid cipher (all caps, no spaces)

    This is the version of the Bifid cipher that uses the 6x6
    Polybius square.

    AUTHOR:
        Alasdair McAndrew

    EXAMPLES:
        sage: key = "encrypt"
        sage: pt = "meet me on monday at 8am"
        sage: bifid(pt, key)
        [[2, 5], [0, 0], [0, 0], [1, 0], [2, 5], [0, 0], [3, 0], 
         [0, 1], [2, 5], [3, 0], [0, 1], [1, 3], [1, 1], [0, 4], 
         [1, 1], [1, 0], [5, 4], [1, 1], [2, 5]]
        'HNHOKNTA5MEPEGNQZYG'

    """
    AS = AlphabeticStrings()
    A = list(AS.alphabet())+[str(x) for x in range(10)]
    # first make sure the letters are capitalized
    # and text has no spaces
    key0 = [x.capitalize() for x in key if not(x.isspace())]
    pt0 = [x.capitalize() for x in pt if not(x.isspace())]
    # create long key
    long_key = key0+[x for x in A if not(x in key0)]
    n = len(pt0)
    # the fractionalization
    pairs = [[long_key.index(x)//6, long_key.index(x)%6] for x in pt0]
    print pairs
    tmp_cipher = flatten([x[0] for x in pairs]+[x[1] for x in pairs])
    ct = join([long_key[6*tmp_cipher[2*i]+tmp_cipher[2*i+1]] for i in range(n)], sep="")
    return ct

def bifid_square(key):
    """
    Produced the Polybius square for the 6x6 Bifid cipher.

    EXAMPLE:
        sage: key = "encrypt"
        sage: bifid_square(key)
        [E N C R Y P]
        [T A B C D E]
        [F G H I J K]
        [L M N O P Q]
        [R S T U V W]
        [X Y Z 0 1 2]

    """
    AS = AlphabeticStrings()
    A = list(AS.alphabet())+[str(x) for x in range(10)]
    # first make sure the letters are capitalized
    # and text has no spaces
    key0 = [SR(x.capitalize()) for x in key if not(x.isspace())]
    # create long key
    long_key = key0+[SR(x) for x in A if not(x in key0)]
    # the fractionalization
    pairs = [[long_key.index(SR(x))//6, long_key.index(SR(x))%6] for x in A]
    f = lambda i,j: long_key[6*i+j] 
    M = matrix(SR, 6, 6, f)
    return M

Have fun!

2012/01/04 Posted by wdjoyner | math, sage, software | Bifid cipher, cipher, math, mathematics, sage, sage math | Leave a Comment

Zombies and Mathematics

What do you do if there is a Zombie attack? Can mathematics help? This post is (humorously) dedicated to collecting links to papers or blog posted related to the mathematical models of Zombies.

George Romero‘s 1968 Night of the Living Dead, now in the public domain, introduced reanimated ghouls, otherwise known as zombies, which craved live human flesh. Romero’s script was insired on Richard Matheson’s I Am Legend. In Romero’s version, the zombies could be killed by destroying the zombie’s brain. A dead human could, in some cases be “reanimated,” turning into a zombie. These conditions are modeled mathematically in several papers, given below.

When Zombies Attack! Mathematical Modelling of a Zombie Outbreak!, paper by Mathematicians at the University of Ottawa, Philip Munz, Ioan Hudea, Joe Imad and Robert J. Smith? (yes, his last name is spelled “Smith?”).
More Zombie publications by Univ. Ottawa Math Dept Prof. Robert Smith?.
Epidemics in the presence of social attraction and repulsion, Oct 2010 Zombie paper by Evelyn Sander and Chad M. Topaz.
Statistical Inference in a Zombie Outbreak Model, slides for a talk given by Des Higman, May 2010.
Mathematics kills zombies dead!, 08/17/2009 blog post by “Zombie Research Society Staff”.
The Mathematics of Zombies, August 18, 2009 blog post by Mike Elliot.
Love, War and Zombies – Systems of Differential Equations using Sage, April 2011 slides by David Joyner.
Sage commands for Love, War and Zombies talk.
This was given on April 22, 2011, as a Project Mosaic/M-cast broadcast.

Public domain 1968 film Night of the Living Dead by George Romero.

2011/05/05 Posted by wdjoyner | cimema, math, sage | Leave a Comment

uniform matroids and MDS codes

It is known that uniform (resp. paving) matroids correspond to MDS (resp. “almost MDS” codes). This post explains this connection.

An MDS code is an $[n,k,d]$ linear error correcting block code $C$ which meets the Singleton bound, $d+k=n+1$ . A uniform matroid is a matroid for which all circuits are of size $\geq r(M)+1$ , where $r(M)$ is the rank of M. Recall, a circuit in a matroid M=(E,J) is a minimal dependent subset of E — that is, a dependent set whose proper subsets are all independent (i.e., all in J).

Consider a linear code $C$ whose check matrix is an $(n-k)\times n$ matrix $H=(\vec{h}_1,\dots , \vec{h}_n)$ . The vector matroid M=M[H] is a matroid for which the smallest sized dependency relation among the columns of H is determined by the check relations $c_1\vec{h}_1 + \dots + c_n \vec{h}_n = H\vec{c}=\vec{0}$ , where $\vec{c}=(c_1,\dots,c_n)$ is a codeword (in C which has minimum dimension d). Such a minimum dependency relation of H corresponds to a circuit of M=M[H].

"""

Implements a Matroid class and some methods.
A finite matroid M is a pair (E, I), where E is a finite set

and I is a collection of subsets of E (called the independent sets)

with the following properties:
 * The empty set is independent.

 * Every subset of an independent set is independent.

 * If A and B are two independent sets of I and A has

   more elements than B, then there exists an element in A

   which is not in B that when added to B still gives an independent set.
A subset of E that is not independent is called dependent.

A maximal independent set is called a basis for the matroid.

Any two bases of a matroid M have the same number of elements.

This number is called the rank of M. A circuit in a matroid M is

a minimal dependent subset of E.
REFERENCES:
http://en.wikipedia.org/wiki/Matroid
"""
class Matroid_generic():

    """

    Implements a generic matroid.
    """

    def __init__(self, X, IndepSets = None, Bases = None):

        """

        A matroid can either be specified by its independent

        sets, its bases or its circuits (probably not both).

        This class provides methods for converting one type

        into the others.
        TODO:

           Add "Circuits = None" option.
        EXAMPLES:

            sage: A = matrix(GF(2), [[1,0,0,1,1,0,1],[0,1,0,1,0,1,1],[0,0,1,0,1,1,1]])

            sage: Fano = vector_matroid(A)

            sage: Fano

            Matroid([0, 1, 2, 3, 4, 5, 6], Bases = [[0, 2, 3],

            [1, 3, 6], [1, 2, 4], [0, 3, 5], [1, 5, 6], [0, 1, 5],

            [2, 3, 4], [0, 1, 4], [0, 1, 2], [3, 5, 6], [3, 4, 6],

            [0, 1, 6], [1, 4, 5], [1, 3, 4], [1, 3, 5], [0, 3, 6],

            [0, 2, 6], [0, 4, 6], [0, 2, 5], [2, 5, 6], [2, 4, 6],

            [1, 2, 3], [0, 3, 4], [0, 4, 5], [2, 3, 5], [1, 2, 6],

            [2, 4, 5], [4, 5, 6]])

            sage: Fano.circuits()

            [[0, 1, 3], [0, 2, 4], [0, 5, 6], [1, 2, 5], [1, 4, 6],

             [2, 3, 6], [3, 4, 5], [0, 1, 2, 6], [0, 1, 4, 5],

             [0, 2, 3, 5], [0, 3, 4, 6], [1, 2, 3, 4], [1, 3, 5, 6],

             [2, 4, 5, 6]]

        """

        self.ground_set = X

        self.bases = Bases

        if Bases != None:

            S = {}

            for B in Bases:

                SB = Set(B).subsets()

                S.union(SB)

            self.indep_sets = [list(x) for x in Set(S)]

            J = Set([Set(x) for x in self.indep_sets])

            D = [list(x) for x in SX if not(x in J)]

            self.dep_sets = D

        if IndepSets != None:

            self.indep_sets = IndepSets

            J = Set([Set(x) for x in self.indep_sets])

            M = max([x.cardinality() for x in J])

            B = [list(x) for x in J if x.cardinality() == M]

            self.bases = B

            SX = Set(X).subsets()

            D = [list(x) for x in SX if not(x in J)]

            self.dep_sets = D

        # now create the circuits - BUG

        C = []

        for x in D:

            addx = [Set([b for b in x if b != a]) in J for a in x]

            if not(False in addx):

                C.append(x)

        self.circuts = [list(x) for x in C]
    def __repr__(self):

        return "Matroid(%s, Bases = %s)"%(self.base_set(), self.bases())
    def __str__(self):

        return "Matroid on base set %s of rank %s"%(self.base_set(), self.rank())
    def basis(self):

        """

        EXAMPLES:

            sage: A = matrix(GF(2), [[1,0,0,1,1,0],[0,1,0,1,0,1],[0,0,1,0,1,1]])

            sage: M = vector_matroid(A)

            sage: M.basis()

            [[0, 2, 3], [1, 2, 4], [0, 3, 5], [0, 1, 5], [2, 3, 4],

             [0, 1, 4], [0, 1, 2], [0, 2, 5], [1, 4, 5], [1, 3, 5],

             [0, 3, 4], [0, 4, 5], [2, 3, 5], [1, 2, 3], [1, 3, 4],

             [2, 4, 5]]

        """

        return self.bases
    def circuits(self):

        """

        Returns the minimal dependent sets.
        EXAMPLES:

            sage: A = matrix(GF(2), [[1,0,0,1,1,0],[0,1,0,1,0,1],[0,0,1,0,1,1]])

            sage: M = vector_matroid(A)

            sage: M.circuits()

            [[0, 1, 3], [0, 2, 4], [1, 2, 5], [3, 4, 5], [0, 1, 4, 5],

             [0, 2, 3, 5], [1, 2, 3, 4]]

        """

        return self.circuts
    def dependent_sets(self):

        return self.dep_sets
    def base_set(self):

        return self.ground_set
    def independent_sets(self):

        S = self.indep_sets

        S.sort()

        return S
    def rank(self):

        r = len(self.bases()[0])

        return r
def vector_matroid(mat):

    """

    Returns the vector matroid defined by the matrix M
    EXAMPLES:

        sage: A = matrix(GF(2), [[1,0,0,1,1,0],[0,1,0,1,0,1],[0,0,1,0,1,1]])

        sage: M = vector_matroid(A)

        sage: M.rank()

        3

        sage: print M

        Matroid on base set [0, 1, 2, 3, 4, 5] of rank 3

        sage: M.base_set()

        [0, 1, 2, 3, 4, 5]

        sage: M.basis()

        [[0, 2, 3], [1, 2, 4], [0, 3, 5], [0, 1, 5], [2, 3, 4],

         [0, 1, 4], [0, 1, 2], [0, 2, 5], [1, 4, 5], [1, 3, 5],

         [0, 3, 4], [0, 4, 5], [2, 3, 5], [1, 2, 3], [1, 3, 4], [2, 4, 5]]

        sage: M.circuits()

        [[1], [2], [2], [4], [1, 3], [1, 4], [1, 5], [1, 4], [1, 5],

         [1, 5], [2, 4], [2, 5], [2, 5], [3, 5], [2, 4], [2, 5], [2, 5],

         [3, 5], [3, 5], [1, 3], [1, 3], [1, 4], [1, 4], [2, 4], [2, 4], [1, 3, 5]]

    """

    F = mat.base_ring()

    n = len(mat.columns())

    k = len(mat.rows())

    J = Combinations(n,k)

    indE = []

    for x in J:

        M = matrix([mat.column(i) for i in x])

        if k == M.rank():     # all indep sets of max size

            indE.append(x)

            for y in powerset(x): # all smaller indep sets

                if not(y in indE):

                    indE.append(y)

    return Matroid_generic(range(n), IndepSets=indE)
def bases(mat, invariant = False):

    """

    Find all bases in a vector/representable matroid.
    EXAMPLES:

        sage: M = matrix(GF(2), [[1,0,0,1,1,0],[0,1,0,1,0,1],[0,0,1,0,1,1]])

        sage: L = [x[1] for x in bases(M, invariant=True)]; L.sort(); L

        [5, 16, 17, 33, 37, 44, 45, 48, 81, 84, 92, 93, 96, 112, 113, 116]

        sage: C = ReedSolomonCode(6,4,GF(7))

        sage: G = C.gen_mat(); G

        [1 1 1 1 1 1]

        [0 1 2 3 4 5]

        [0 1 4 2 2 4]

        [0 1 1 6 1 6]

        sage: B = bases(G); B

        [[0, 1, 2, 3], [0, 1, 2, 4], [0, 1, 2, 5], [0, 1, 3, 4],

         [0, 1, 3, 5], [0, 1, 4, 5], [0, 2, 3, 4], [0, 2, 3, 5],

         [0, 2, 4, 5], [0, 3, 4, 5], [1, 2, 3, 4], [1, 2, 3, 5],

         [1, 2, 4, 5], [1, 3, 4, 5], [2, 3, 4, 5]]

        sage: len(B)

        15

        sage: binomial(6,4)

        15

    """

    r = mat.rank()

    ind = independent_sets(mat)

    B = []

    for x in ind:

        if len(x) == r:

            B.append(x)

    B.sort()

    if invariant == True:

        B = [(b,sum([k*2^(k-1) for k in b])) for b in B]

    return B
def circuits(mat):

    """

    Finds all circuits in a vector matroid.
    EXAMPLES:

        sage: M = matrix(GF(2), [[1,0,0,1,1,0],[0,1,0,1,0,1],[0,0,1,0,1,1]])

        sage: M

        [1 0 0 1 1 0]

        [0 1 0 1 0 1]

        [0 0 1 0 1 1]

        sage: circuits(M)

        [[0, 1, 3], [0, 2, 4], [1, 2, 3, 4], [1, 2, 5], [0, 2, 3, 5],

        [0, 1, 4, 5], [3, 4, 5]]

        sage: M = matrix(GF(3), [[1,0,0,1,1,0],[0,1,0,1,0,1],[0,0,1,0,1,1]])

        sage: circuits(M)

        [[0, 1, 3], [0, 2, 4], [1, 2, 3, 4], [1, 2, 5], [0, 2, 3, 5],

        [0, 1, 4, 5], [0, 3, 4, 5], [1, 3, 4, 5], [2, 3, 4, 5]]
    """

    J = independent_sets(mat)

    n = len(mat.columns())

    E = set(range(n))

    subsetsE = powerset(E)

    C = []

    for S in subsetsE:

        is_circuit = True

        if not(S in J):

            for x in powerset(S):

                if xS and not(x in J):

                    is_circuit = False

        if is_circuit == True:

            if not(S in J):

                C.append(S)

    return C

                

                2010/03/10 

								
                  Posted by               	  wdjoyner |
                  math, sage, software				   |
                  3 Comments

Floyd-Warshall-Roy, 3

As in the previous post, let $G=(V,E)$ be a weighted digraph having $n$ vertices and let $A=A_G$ denote its $n\times n$ adjacency matrix. We identify the vertices with the set $\{1,2,\dots, n\}$ .

The previous post discussed the following result, due to Sturmfels et al.

Theorem: The entry of the matrix A^n-1 in row i and column j equals the length of a shortest path from vertex i to vertex j in G. (Here A^n-1 denotes the n-1-st tropical power of A.)

This post discusses an implementation in Python/Sage.

Consider the following class definition.

class TropicalNumbers:
    """
    Implements the tropical semiring.

    EXAMPLES:
        sage: T = TropicalNumbers()
        sage: print T
        Tropical Semiring

    """
    def __init__(self):
        self.identity = Infinity

    def __repr__(self):
        """
        Called to compute the "official" string representation of an object.
        If at all possible, this should look like a valid Python expression
        that could be used to recreate an object with the same value.

        EXAMPLES:
            sage: TropicalNumbers()
            TropicalNumbers()

        """
        return "TropicalNumbers()"

    def __str__(self):
        """
        Called to compute the "informal" string description of an object. 

        EXAMPLES:
            sage: T = TropicalNumbers()
            sage: print T
            Tropical Semiring

        """
        return "Tropical Semiring"

    def __call__(self, a):
        """
        Coerces a into the tropical semiring.

        EXAMPLES:
            sage: T(10)
            TropicalNumber(10)
            sage: print T(10)
            Tropical element 10 in Tropical Semiring
        """
        return TropicalNumber(a)

    def __contains__(self, a):
        """
        Implements "in".

        EXAMPLES:
            sage: T = TropicalNumbers()
            sage: a = T(10)
            sage: a in T
            True

        """
        if a in RR or a == Infinity:
            return a==Infinity or (RR(a) in RR)
        else:
            return a==Infinity or (RR(a.element) in RR)

class TropicalNumber:
    def __init__(self, a):
        self.element = a
        self.base_ring = TropicalNumbers()

    def __repr__(self):
        """
        Called to compute the "official" string representation of an object.
        If at all possible, this should look like a valid Python expression
        that could be used to recreate an object with the same value.

        EXAMPLES:

        """
        return "TropicalNumber(%s)"%self.element

    def __str__(self):
        """
        Called to compute the "informal" string description of an object. 

        EXAMPLES:
            sage: T = TropicalNumbers()
            sage: print T(10)
            Tropical element 10 in Tropical Semiring
        """
        return "%s"%(self.number())

    def number(self):
        return self.element

    def __add__(self, other):
        """
        Implements +. Assumes both self and other are instances of
        TropicalNumber class.

        EXAMPLES:
            sage: T = TropicalNumbers()
            sage: a = T(10)
            sage: a in T
            True
            sage: b = T(15)
            sage: a+b
            10

        """
        T = TropicalNumbers()
        return T(min(self.element,other.element))

    def __mul__(self, other):
        """
        Implements multiplication *.

        EXAMPLES:
            sage: T = TropicalNumbers()
            sage: a = T(10)
            sage: a in T
            True
            sage: b = T(15)
            sage: a*b
            25
        """
        T = TropicalNumbers()
        return T(self.element+other.element)

class TropicalMatrix:
    def __init__(self, A):
        T = TropicalNumbers()
        self.base_ring = T
        self.row_dimen = len(A)
        self.column_dimen = len(A[0])
        # now we coerce the entries into T
        A0 = A
        m = self.row_dimen
        n = self.column_dimen
        for i in range(m):
            for j in range(n):
                A0[i][j] = T(A[i][j])
        self.array = A0

    def matrix(self):
        """
        Returns the entries (as ordinary numbers).

        EXAMPLES:
            sage: A = [[0,1,3,7],[2,0,1,3],[4,5,0,1],[6,3,1,0]]
            sage: AT = TropicalMatrix(A)
            sage: AT.matrix()
            [[0, 1, 3, 7], [2, 0, 1, 3], [4, 5, 0, 1], [6, 3, 1, 0]]
        """
        m = self.row_dim()
        n = self.column_dim()
        A0 = [[0 for i in range(n)] for j in range(m)]
        for i in range(m):
            for j in range(n):
                A0[i][j] = (self.array[i][j]).number()
        return A0

    def row_dim(self):
        return self.row_dimen

    def column_dim(self):
        return self.column_dimen

    def __repr__(self):
        """
        Called to compute the "official" string representation of an object.
        If at all possible, this should look like a valid Python expression
        that could be used to recreate an object with the same value.

        EXAMPLES:

        """
        return "TropicalMatrix(%s)"%self.array

    def __str__(self):
        """
        Called to compute the "informal" string description of an object. 

        EXAMPLES:

        """
        return "Tropical matrix %s"%(self.matrix())

    def __add__(self, other):
        """
        Implements +. Assumes both self and other are instances of
        TropicalMatrix class.

        EXAMPLES:
            sage: A = [[1,2,Infinity],[3,Infinity,0]]
            sage: B = [[2,Infinity,1],[3,-1,1]]
            sage: AT = TropicalMatrix(A)
            sage: BT = TropicalMatrix(B)
            sage: AT
            TropicalMatrix([[TropicalNumber(1), TropicalNumber(2), TropicalNumber(+Infinity)],
              [TropicalNumber(3), TropicalNumber(+Infinity), TropicalNumber(0)]])
            sage: AT+BT
            [[TropicalNumber(1), TropicalNumber(2), TropicalNumber(1)],
             [TropicalNumber(3), TropicalNumber(-1), TropicalNumber(0)]]
        """
        A = self.array
        B = other.array
        C = []
        m = self.row_dim()
        n = self.column_dim()
        if m != other.row_dim:
            raise ValueError, "Row dimensions must be equal."
        if n != other.column_dim:
            raise ValueError, "Column dimensions must be equal."
        for i in range(m):
            row = [A[i][j]+B[i][j] for j in range(n)] # + as tropical numbers
            C.append(row)
        return C

    def __mul__(self, other):
        """
        Implements multiplication *.

        EXAMPLES:
            sage: A = [[1,2,Infinity],[3,Infinity,0]]
            sage: AT = TropicalMatrix(A)
            sage: B = [[2,Infinity],[-1,1],[Infinity,0]]
            sage: BT = TropicalMatrix(B)
            sage: AT*BT
             [[TropicalNumber(1), TropicalNumber(3)],
              [TropicalNumber(5), TropicalNumber(0)]]
            sage: A = [[0,1,3,7],[2,0,1,3],[4,5,0,1],[6,3,1,0]]
            sage: AT = TropicalMatrix(A)
            sage: A = [[0,1,3,7],[2,0,1,3],[4,5,0,1],[6,3,1,0]]
            sage: AT = TropicalMatrix(A)
            sage: print AT*AT*AT
            Tropical matrix [[0, 1, 2, 3], [2, 0, 1, 2], [4, 4, 0, 1], [5, 3, 1, 0]]
        """
        T = TropicalNumbers()
        A = self.matrix()
        B = other.matrix()
        C = []
        mA = self.row_dim()
        nA = self.column_dim()
        mB = other.row_dim()
        nB = other.column_dim()
        if nA != mB:
            raise ValueError, "Column dimension of A and row dimension of B must be equal."
        for i in range(mA):
            row = []
            for j in range(nB):
                c = T(Infinity)
                for k in range(nA):
                    c = c+T(A[i][k])*T(B[k][j])
                row.append(c.number())
            C.append(row)
        return TropicalMatrix(C)

This shows that the shortest distances of digraph with adjacency matrix $\begin{pmatrix} 0&1&3&7\\ 2&0&1&3\\ 4&5&0&1\\ 6&3&1&0\end{pmatrix},$ is equal to A³, which is equal to $\begin{pmatrix} 0&1&2&3\\ 2&0&1&2\\ 4&4&0&1\\ 5&3&1&0\end{pmatrix}.$ This verifies an example given in chapter 1 of the book by Maclagan and Sturmfels, Introduction to Tropical Geometry .

2010/03/02 Posted by wdjoyner | math, sage, software | graph theory | Leave a Comment

Floyd-Warshall-Roy

The Floyd-Warshall-Roy algorithm is an algorithm for finding shortest paths in a weighted, directed graph. It allows for negative edge weights and detects a negative weight cycle if one exists. Assuming that there are no negative weight cycles, a single execution of the FWR algorithm will find the shortest paths between all pairs of vertices.

This algorithm is an example of dynamic programming, which allows one to break the computation down to simpler steps using some sort of recursive procedure. If $n=|V|$ then this is a $O(n^3)$ -time algorithm. For comparison, the Bellman-Ford algorithm has complexity $O(|V|\cdot |E|)$ , which is $O(n^3)$ -time for “dense” graphs. However, Bellman-Ford only yields the shortest paths eminating from a single vertex. To achieve comparable output, we would need to iterate Bellman-Ford over all vertices, which would be an $O(n^4)$ -time algorithm for “dense” graphs. Except for “sparse” graphs, Floyd-Warshall-Roy is better than an interated implementation of Bellman-Ford.

Here is a Sage implementation

def floyd_warshall_roy(A):
    """
    Shortest paths

    INPUT: 
        A - weighted adjacency matrix 

    OUTPUT
        dist  - a matrix of distances of shortest paths
        paths - a matrix of shortest paths

    EXAMPLES:
        sage: A = matrix([[0,1,2,3],[0,0,2,1],[20,10,0,3],[11,12,13,0]]); A
        sage: floyd_warshall_roy(A)
        ([[0, 1, 2, 2], [12, 0, 2, 1], [14, 10, 0, 3], [11, 12, 13, 0]], 
          [-1  1  2  1]
          [ 3 -1  2  3]
          [ 3 -1 -1  3]
          [-1 -1 -1 -1])
        sage: A = matrix([[0,1,2,4],[0,0,2,1],[0,0,0,5],[0,0,0,0]])
        sage: floyd_warshall_roy(A)
        ([[0, 1, 2, 2], [+Infinity, 0, 2, 1], [+Infinity, +Infinity, 0, 5],
          [+Infinity, +Infinity, +Infinity, 0]], 
          [-1  1  2  1]
          [-1 -1  2  3]
          [-1 -1 -1  3]
          [-1 -1 -1 -1])
        sage: A = matrix([[0,1,2,3],[0,0,2,1],[-5,0,0,3],[1,0,1,0]]); A
        sage: floyd_warshall_roy(A)
        Traceback (click to the left of this block for traceback)
        ...
        ValueError: A negative edge weight cycle exists.
        sage: A = matrix([[0,1,2,3],[0,0,2,1],[-1/2,0,0,3],[1,0,1,0]]); A
        sage: floyd_warshall_roy(A)
        ([[0, 1, 2, 2], [3/2, 0, 2, 1], [-1/2, 1/2, 0, 3/2], [1/2, 3/2, 1, 0]],
          [-1  1  2  1]
          [ 2 -1  2  3]
          [-1  0 -1  1]
          [ 2  2 -1 -1])
    """
    G = Graph(A, format="weighted_adjacency_matrix")
    V = G.vertices()
    E = [(e[0],e[1]) for e in G.edges()]
    n = len(V)
    dist = [[0]*n for i in range(n)]
    paths = [[-1]*n for i in range(n)]
    # initialization step
    for i in range(n):
        for j in range(n):
            if (i,j) in E:
                paths[i][j] = j
            if i == j:
                dist[i][j] = 0
            elif A[i][j] != 0:
                dist[i][j] = A[i][j]
            else:
                dist[i][j] = infinity
    # iteratively finding the shortest path
    for j in range(n):
        for i in range(n):
            if i != j:
                for k in range(n):
                    if k != j:
                        if dist[i][k]>dist[i][j]+dist[j][k]:
                            paths[i][k] = V[j]
                        dist[i][k] = min(dist[i][k], dist[i][j] +dist[j][k])
    for i in range(n):
        if dist[i][i] < 0:
            raise ValueError, "A negative edge weight cycle exists."
    return dist, matrix(paths)

2010/02/08 Posted by wdjoyner | math, sage | "shortest path", graphs, sage | 2 Comments

« Previous Entries

May 2012

M T W T F S S

« Apr

1 2 3 4 5 6

7 8 9 10 11 12 13

14 15 16 17 18 19 20

21 22 23 24 25 26 27

28 29 30 31
Blogroll
math
sage
Recent Posts

May 2012
M	T	W	T	F	S	S
« Apr
	1	2	3	4	5	6
7	8	9	10	11	12	13
14	15	16	17	18	19	20
21	22	23	24	25	26	27
28	29	30	31

Archives
- April 2012 (3)
- February 2012 (3)
- January 2012 (3)
- November 2011 (1)
- August 2011 (1)
- July 2011 (1)
- June 2011 (1)
- May 2011 (1)
- April 2011 (1)
- August 2010 (1)
- March 2010 (2)
- February 2010 (4)
Categories
- cimema
- gsoc
- math
- sage
- software
- sympy
- Uncategorized
RSS
Entries RSS
Comments RSS

Yet Another Mathblog

Sestinas and Sage

Remarks on the Hamming [7.4.3] code and Sage

Digital steganography and Sage

Applications of graphs to Boolean functions

The Vigenère cipher and Sage

Bifid cipher and Sage

Zombies and Mathematics

uniform matroids and MDS codes

Floyd-Warshall-Roy, 3

Floyd-Warshall-Roy

About

Blogroll

math

sage

Recent Posts

Archives

Categories

RSS

Site info

Follow “Yet Another Mathblog”