Helping Nancy and Esther with a math problem

The Nancy comic strip for June 12 2018 has a math problem which Nancy and Esther seemed to be intimidated by; I figured I would try to help out a little.
The problem consists of two equations in two unknowns; we are asked to find values for the unknowns which make both equations true at the same time
x2 + y2 = 3
16x2 – 4y2 = 0
(A geometric aside: the first equation describes an ellipse, and the second equation describes a degenerate hyperbola, in particular, the two lines y = 2x and y = -2x. The common solutions, if any exist, are those points where the ellipse intersects the hyperbola.)
We can simplify the second equation by dividing throughout by 4
x2 + y2 = 3
4x2y2 = 0
Add the two equations to get a single equation in x alone
5x2 = 3
Solving this for x, there are two solutions: x = √(3/5) (about 0.77) and x = -√(3/5)
Plug either of these values for x into 4x2y2 = 0
4(√(3/5))2y2 = 0
y2 = 12/5
Solving this for y, there are two solutions: y = 2√(3/5) (about 1.55) and y = -2√(3/5)
Plugging in and checking:
3/5 + 12/5 = 15/5 = 3
16(3/5) – 4(12/5) = 48/5 – 48/5 = 0
So there are four solutions:
  • (xy) = (√(3/5), 2√(3/5))
  • (xy) = (-√(3/5), 2√(3/5))
  • (xy) = (√(3/5), -2√(3/5))
  • (xy) = (-√(3/5), -2√(3/5))

How far apart can squares be?

Consider the squares: 0, 1, 4, 9, 16, 25, 36…

“No two squares are 6 apart,” I say. After some subtractions you believe me.

Exercise: prove no two squares are 6 apart.

“Nor are any two squares 134 apart,” I say. You look at me in surprise. After some puzzlement, inspiration strikes.

Exercise: prove that any two squares differ either by an odd number or by multiple of 4.

“In fact,” I say, “the set {2, 6, 10, 14, …} = {4k + 2} completely describes how far apart two squares cannot be…”

Exercise: given any n not of the form 4k + 2, prove that it is possible to find two squares that are n apart.

“… with one exception.”

Exercise: find the exception. Find what was wrong with your previous proof. Convince yourself the proof is now correct and there are no other exceptions.

Perl scripts to encrypting and decrypt text using Rijndael

I talked about Rijndael in a few previous posts: Expressing a function f: GF(2⁸) → GF(2⁸) as a polynomial using a Lagrange polynomial, Generating the Rijndael S-box, Efficient multiplication and division in GF(28), Sieving irreducible monic polynomials over a finite field, Addition and multiplication table for GF(22).

Today I wrote a couple of Perl scripts which use off-the-shelf CPAN modules to encrypt and decrypt text using a passphrase. For ease of transport, at the cost of some storage space, the encrypted payload is Base64 encoded.

The particular CPAN modules in question are Crypt::CBC and Crypt::Rijndael.

encrypt.pl and decrypt.pl source in my GitHub repository. They operate on STDIN and produce output on STDOUT.

Note that if you pass an incorrect passphrase to decrypt.pl you get garbage output, rather than an error.

Here is a secret message which is encoded with the passphrase “hunter2”: encrypted.txt Right-click the link and save, rather than opening it in the browser.

Expressing a function f: GF(2⁸) → GF(2⁸) as a polynomial using a Lagrange polynomial

I talked about Rijndael in a couple of previous posts: Generating the Rijndael S-box, Efficient multiplication and division in GF(28), Sieving irreducible monic polynomials over a finite field, Addition and multiplication table for GF(22).

I’m going to talk some more about it today.

The Rijndael non-linear S-box S(x) is a composition of two invertible functions f(g(x)). Last time we showed how to generate these, and their inverses, as 256-entry tables.

As Daemen and Rijmen point out, any function from a finite field to itself can be expressed as a polynomial. In fact, given a tabular form of the function, it is possible to generate the Lagrange polynomial and then simplify.

They also give the polynomial for S(x):

SRD[x] = 05·x255 + 09·x253 + F9·x251 + 25·x247 + F4·x239 + 01·x223 + B5·x191 + 8F·x127 + 63

Well, let’s check this.
And while we’re at it, let’s find the polynomials for g(x), f(x) and even S-1(x) too.
First of all, let’s start with the Lagrange polynomial.

Given a table of entries { (00, y00), (01, y01), …, (ff, yff) }, there is a polynomial L(x) which gives the same output, namely:

L(x) = Σi ∈ { 00, 01, …, ff } yi pi(x)
where pi(x) = Πj ∈ { 00, 01, …, ff }, ji (xj) / (ij)

Can we simplify this?
Yes. Note that (ij)-1 term varies over all the non-zero elements of the finite field. Since this is a field, every non-zero element has an inverse, which might or might not be itself.
If the inverse is not itself, we can pair the two inverse elements together, and we get 01, which is the multiplicative identity, so we can ignore it.
What are we left with? The product of those non-zero elements which are their own inverses.
In the case of GF(28), there is only one such element, namely 01.
Great! We can ignore the (ij)-1 terms altogether:

L(x) = Σi ∈ { 00, 01, …, ff } yi Πj ∈ { 00, 01, …, ff }, ji (xj)

What is this? A sum of 256 different polynomials.
Each of these summands is a product of 255 polynomials of degree 1, and so is a polynomial of degree 255.
But in GF(28), x255 = 01 for all x. So each of the summands is actually a polynomial of degree 254, or less.
So the sum is also a polynomial of degree 254, or less; terms with an exponent >= 255 go away and just contribute to lower-order terms.

Great. We have { (00, y00), (01, y01), …, (ff, yff) } for f, g, S, and S-1. Let’s plug them in and see what happens!

I wrote a Perl script to do this; the script is attached. Here’s the output. (It’s quite slow; it takes about two and a half minutes to run.)

g(x) = 01 x^254

f(x) =
63 + 05 x + 09 x^2 + f9 x^4 + 25 x^8 + f4 x^16 + 01 x^32 + b5 x^64 +
8f x^128

S(x) =
63 + 8f x^127 + b5 x^191 + 01 x^223 + f4 x^239 + 25 x^247 + f9 x^251 + 09 x^253 +
05 x^254

S^(-1)(x) =
52 + f3 x + 7e x^2 + 1e x^3 + 90 x^4 + bb x^5 + 2c x^6 + 8a x^7 +
1c x^8 + 85 x^9 + 6d x^10 + c0 x^11 + b2 x^12 + 1b x^13 + 40 x^14 + 23 x^15 +
f6 x^16 + 73 x^17 + 29 x^18 + d9 x^19 + 39 x^20 + 21 x^21 + cf x^22 + 3d x^23 +
9a x^24 + 8a x^25 + 2f x^26 + cf x^27 + 7b x^28 + 04 x^29 + e8 x^30 + c8 x^31 +
85 x^32 + 7b x^33 + 7c x^34 + af x^35 + 86 x^36 + 2f x^37 + 13 x^38 + 65 x^39 +
75 x^40 + d3 x^41 + 6d x^42 + d4 x^43 + 89 x^44 + 8e x^45 + 65 x^46 + 05 x^47 +
ea x^48 + 77 x^49 + 50 x^50 + a3 x^51 + c5 x^52 + 01 x^53 + 0b x^54 + 46 x^55 +
bf x^56 + a7 x^57 + 0c x^58 + c7 x^59 + 8e x^60 + f2 x^61 + b1 x^62 + cb x^63 +
e5 x^64 + e2 x^65 + 10 x^66 + d1 x^67 + 05 x^68 + b0 x^69 + f5 x^70 + 86 x^71 +
e4 x^72 + 03 x^73 + 71 x^74 + a6 x^75 + 56 x^76 + 03 x^77 + 9e x^78 + 3e x^79 +
19 x^80 + 18 x^81 + 52 x^82 + 16 x^83 + b9 x^84 + d3 x^85 + 38 x^86 + d9 x^87 +
04 x^88 + e3 x^89 + 72 x^90 + 6b x^91 + ba x^92 + e8 x^93 + bf x^94 + 9d x^95 +
1d x^96 + 5a x^97 + 55 x^98 + ff x^99 + 71 x^100 + e1 x^101 + a8 x^102 + 8e x^103 +
fe x^104 + a2 x^105 + a7 x^106 + 1f x^107 + df x^108 + b0 x^109 + 03 x^110 + cb x^111 +
08 x^112 + 53 x^113 + 6f x^114 + b0 x^115 + 7f x^116 + 87 x^117 + 8b x^118 + 02 x^119 +
b1 x^120 + 92 x^121 + 81 x^122 + 27 x^123 + 40 x^124 + 2e x^125 + 1a x^126 + ee x^127 +
10 x^128 + ca x^129 + 82 x^130 + 4f x^131 + 09 x^132 + aa x^133 + c7 x^134 + 55 x^135 +
24 x^136 + 6c x^137 + e2 x^138 + 58 x^139 + bc x^140 + e0 x^141 + 26 x^142 + 37 x^143 +
ed x^144 + 8d x^145 + 2a x^146 + d5 x^147 + ed x^148 + 45 x^149 + c3 x^150 + ec x^151 +
1c x^152 + 3e x^153 + 2a x^154 + b3 x^155 + 9e x^156 + b7 x^157 + 38 x^158 + 82 x^159 +
23 x^160 + 2d x^161 + 87 x^162 + ea x^163 + da x^164 + 45 x^165 + 24 x^166 + 03 x^167 +
e7 x^168 + 19 x^169 + e3 x^170 + d3 x^171 + 4e x^172 + dd x^173 + 11 x^174 + 4e x^175 +
81 x^176 + 91 x^177 + 91 x^178 + 59 x^179 + a3 x^180 + 80 x^181 + 92 x^182 + 7e x^183 +
db x^184 + c4 x^185 + 20 x^186 + ec x^187 + db x^188 + 55 x^189 + 7f x^190 + a8 x^191 +
c1 x^192 + 64 x^193 + ab x^194 + 1b x^195 + fd x^196 + 60 x^197 + 05 x^198 + 13 x^199 +
2c x^200 + a9 x^201 + 76 x^202 + a5 x^203 + 1d x^204 + 32 x^205 + 8e x^206 + 1e x^207 +
c0 x^208 + 65 x^209 + cb x^210 + 8b x^211 + 93 x^212 + e4 x^213 + ae x^214 + be x^215 +
5f x^216 + 2c x^217 + 3b x^218 + d2 x^219 + 0f x^220 + 9f x^221 + 42 x^222 + cc x^223 +
6c x^224 + 80 x^225 + 68 x^226 + 43 x^227 + 09 x^228 + 23 x^229 + c5 x^230 + 6d x^231 +
1d x^232 + 18 x^233 + bd x^234 + 5e x^235 + 1b x^236 + b4 x^237 + 85 x^238 + 49 x^239 +
bc x^240 + 0d x^241 + 1f x^242 + a6 x^243 + 6b x^244 + d8 x^245 + 22 x^246 + 01 x^247 +
7a x^248 + c0 x^249 + 55 x^250 + 16 x^251 + b3 x^252 + cf x^253 + 05 x^254

Daemen and Rijmen’s polynomial rendition of S(x) is confirmed.

Also note how simple g(x) is, and how complex S-1(x) is.

Finally, I must confess that none of this is actually useful. The tabular form is much more convenient for applications. This is just a fun theoretical exercise.

(Well, I had fun.)

EDIT: 2015-10-31 moved script to https://github.com/mvaneerde/blog/blob/master/rijndael/lagrange-interpolation.pl

Generating the Rijndael S-box

I talked about Rijndael in a couple of previous posts: Efficient multiplication and division in GF(28), Sieving irreducible monic polynomials over a finite field, Addition and multiplication table for GF(22).

I’m going to talk some more about it today.

One of the more interesting steps used in the Rijndael transformation is the non-linear S-box. This is a function S(a) that takes an element of GF(28) (which could be represented as a byte) and returns another one. It is invertible but non-linear.

The spec defines S(a) in terms of two other invertible functions g(a) and f(a). In particular S(a) = f(g(a)). It follows that S-1(a) = g-1(f-1(a)).

g(a) is the non-linear piece, and is quite simple:

g(00) is defined as 00.
g(a) = a-1 for all other a in GF(28). In particular, if a = 03b, then g(a) = 03255 – b.

This is clearly invertible; in fact, it is its own inverse.

Daemen and Rijmen seem to have been almost embarrassed by the simplicity of g so they introduced f. To quote from section 3.4.1:

By definition, g has a very simple algebraic expression. This could allow algebraic manipulations that can be used to mount attacks such as interpolation attacks. Therefore, we built the S-box as the sequence of g and an invertible affine transformation f.

I don’t know if I buy the “simplicity” argument. It seems to me that if Rijndael, without f, is robust, then it’s robust. And if you add f to a non-robust scheme, I don’t understand how that makes it robust; I do see that it complicates analysis, but that seems like a drawback rather than an advantage.

But I’ll play along for now. What is f?

b = f(a) is defined using the following matrix multiplication over GF(2) (each entry can be represented as a bit; a row or column as a byte.) Multiplication can be implemented as bitwise AND, and addition by bitwise XOR.

S-box-f

In particular, note that f(00) = 63. So S(00) = f(g(00)) = f(00) = 63. So S-1(63) = 00.

The reference implementation in the book uses hardcoded 256-byte lookup tables for S(a) and S-1(a). I wrote a Perl script which generates these lookup tables and prints them out. This script is attached.

Here’s its output, which matches the listing in the book:

>perl -w s-box.pl
S(xy)
| x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 xa xb xc xd xe xf
—+————————————————
0y | 63 7c 77 7b f2 6b 6f c5 30 01 67 2b fe d7 ab 76
1y | ca 82 c9 7d fa 59 47 f0 ad d4 a2 af 9c a4 72 c0
2y | b7 fd 93 26 36 3f f7 cc 34 a5 e5 f1 71 d8 31 15
3y | 04 c7 23 c3 18 96 05 9a 07 12 80 e2 eb 27 b2 75
4y | 09 83 2c 1a 1b 6e 5a a0 52 3b d6 b3 29 e3 2f 84
5y | 53 d1 00 ed 20 fc b1 5b 6a cb be 39 4a 4c 58 cf
6y | d0 ef aa fb 43 4d 33 85 45 f9 02 7f 50 3c 9f a8
7y | 51 a3 40 8f 92 9d 38 f5 bc b6 da 21 10 ff f3 d2
8y | cd 0c 13 ec 5f 97 44 17 c4 a7 7e 3d 64 5d 19 73
9y | 60 81 4f dc 22 2a 90 88 46 ee b8 14 de 5e 0b db
ay | e0 32 3a 0a 49 06 24 5c c2 d3 ac 62 91 95 e4 79
by | e7 c8 37 6d 8d d5 4e a9 6c 56 f4 ea 65 7a ae 08
cy | ba 78 25 2e 1c a6 b4 c6 e8 dd 74 1f 4b bd 8b 8a
dy | 70 3e b5 66 48 03 f6 0e 61 35 57 b9 86 c1 1d 9e
ey | e1 f8 98 11 69 d9 8e 94 9b 1e 87 e9 ce 55 28 df
fy | 8c a1 89 0d bf e6 42 68 41 99 2d 0f b0 54 bb 16

S^(-1)(xy)
| x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 xa xb xc xd xe xf
—+————————————————
0y | 52 09 6a d5 30 36 a5 38 bf 40 a3 9e 81 f3 d7 fb
1y | 7c e3 39 82 9b 2f ff 87 34 8e 43 44 c4 de e9 cb
2y | 54 7b 94 32 a6 c2 23 3d ee 4c 95 0b 42 fa c3 4e
3y | 08 2e a1 66 28 d9 24 b2 76 5b a2 49 6d 8b d1 25
4y | 72 f8 f6 64 86 68 98 16 d4 a4 5c cc 5d 65 b6 92
5y | 6c 70 48 50 fd ed b9 da 5e 15 46 57 a7 8d 9d 84
6y | 90 d8 ab 00 8c bc d3 0a f7 e4 58 05 b8 b3 45 06
7y | d0 2c 1e 8f ca 3f 0f 02 c1 af bd 03 01 13 8a 6b
8y | 3a 91 11 41 4f 67 dc ea 97 f2 cf ce f0 b4 e6 73
9y | 96 ac 74 22 e7 ad 35 85 e2 f9 37 e8 1c 75 df 6e
ay | 47 f1 1a 71 1d 29 c5 89 6f b7 62 0e aa 18 be 1b
by | fc 56 3e 4b c6 d2 79 20 9a db c0 fe 78 cd 5a f4
cy | 1f dd a8 33 88 07 c7 31 b1 12 10 59 27 80 ec 5f
dy | 60 51 7f a9 19 b5 4a 0d 2d e5 7a 9f 93 c9 9c ef
ey | a0 e0 3b 4d ae 2a f5 b0 c8 eb bb 3c 83 53 99 61
fy | 17 2b 04 7e ba 77 d6 26 e1 69 14 63 55 21 0c 7d

And here’s the interesting part, the implementation of g(a) and f(a):

sub g($) {
my ($a) = shift;

# g(0) we define to be 0
return 0 unless $a;

# otherwise g(a) = a^(-1)
# a = (x + 1)^loga
# so a^(-1) = (x + 1)^(-loga) = (x + 1)^(255 – loga)
my $loga = $log_xplusone_of[$a];
my $logb = 255 – $loga;
$logb -= 255 if $logb >= 255;

return $xplusone_to[$logb];
}

# f(a) = b is defined as follows:
#
# [ b7 ]   ( [ 1 1 1 1 1 0 0 0 ] [ a7 ] )   [ 0 ]
# [ b6 ]   ( [ 0 1 1 1 1 1 0 0 ] [ a6 ] )   [ 1 ]
# [ b5 ]   ( [ 0 0 1 1 1 1 1 0 ] [ a5 ] )   [ 1 ]
# [ b4 ] = ( [ 0 0 0 1 1 1 1 1 ] [ a4 ] ) + [ 0 ]
# [ b3 ]   ( [ 1 0 0 0 1 1 1 1 ] [ a3 ] )   [ 0 ]
# [ b2 ]   ( [ 1 1 0 0 0 1 1 1 ] [ a2 ] )   [ 0 ]
# [ b1 ]   ( [ 1 1 1 0 0 0 1 1 ] [ a1 ] )   [ 1 ]
# [ b0 ]   ( [ 1 1 1 1 0 0 0 1 ] [ a0 ] )   [ 1 ]
#
# where the + is XOR
sub f($) {
my ($a) = @_;

# start with the addition
my $b = 0x63; # 0b01100011;

# do the matrix multiplication
# one matrix column at a time
for (
my ($c, $i) = (0x8f, 0x80); # 0b10001111, 0b10000000
$i;
$c = rotate_byte_right($c), $i >>= 1
) {
# i is used to select a bit out of the a column
# c is the matrix column which is multiplied by that bit
# the resulting product influences the eventual b column

# printf(“%02x %02xn”, $c, $i);

# if this bit in the a column is 0, all of the products are 0, so don’t bother
next unless $a & $i;

# this bit in the a column is 1
# so XOR b with the matrix column
$b ^= $c;
}

return $b;
}

EDIT 2015-10-31: moved script to https://github.com/mvaneerde/blog/blob/master/rijndael/s-box.pl

Efficient multiplication and division in GF(2⁸)

I talked about Rijndael in a couple of previous posts: Sieving irreducible monic polynomials over a finite field, Addition and multiplication table for GF(2²)

I’m going to talk some more about it today.

Rijndael does a lot of arithmetic operations (addition and multiplication) on elements of GF(28)4.

These are represented as polynomials b0 + b1 x + b2 x2 + b3 x3.

The individual coefficients bi are elements of GF(28).

These are themselves polynomials a0 + a1 x + … + a7 x7 where the individual coefficients ai are bits (0 or 1).

Multiplication bi bj is made into a closed operation by taking the modulus of the product under the reduction polynomial m = x8 + x4 + x3 + x + 1.

The reduction polynomial is prime, which is important in establishing that this representation of the Galois field really is a field.

For convenience we will represent the elements of GF(28) as hexadecimal numbers 0x00, 0x01, … 0xfe, 0xff.

By contrast, we will represent the elements of GF(28)4 as vectors (b0 b1 b2 b3); for example, (0x00 0x01 0x02 0x03).

Multiplying vectors in GF(28)4 induces a bunch of multiplications in GF(28). It would be good to come up with a way to make this fast. Doing a polynomial reduction every time is not fast.

Of course, one way to do this is to create a multiplication table with 256 rows and 256 columns, do each multiplication the slow way once, and compile in the results. This would require on the order of 64 kB. But there’s a better way.

The trick that is used in Daemen and Rijmen’s implementation is to find a primitive element e in GF(28) so that the orbit of e (that is to say, the set {e0, e1, …, e254}) spans all of GF(28) – {0x00}.

If we find such an e, then we can save a cheat sheet which is almost as fast as the multiplication table, but requires storing only on the order of 256 bytes.

Well, let’s calculate the orbits of all the elements in order until we find a primitive element.

Orbit(0x00) = { 0x000, 0x001, 0x002, … } = { 0x01, 0x00, 0x00, … } = { 0x00, 0x01 }. Nope.
Orbit(0x01) = { 0x010, 0x011, 0x012, … } = { 0x01, 0x01, 0x01, … } = { 0x01 }. Even worse.
Orbit(0x02) = { 0x020, 0x021, 0x022, … } = { 0x01, 0x02, 0x04, …, 0x8d }. This looks promising at first but we end up with an orbit of only 51 out of the necessary 255 elements (note that 51 divides 255:)

0x020 = 0x01
0x021 = 0x02
0x022 = 0x04
0x023 = 0x08
0x024 = 0x10
0x025 = 0x20
0x026 = 0x40
0x027 = 0x80
0x028 = 0x1b
0x029 = 0x36
0x0210 = 0x6c
0x0211 = 0xd8
0x0212 = 0xab
0x0213 = 0x4d
0x0214 = 0x9a
0x0215 = 0x2f
0x0216 = 0x5e
0x0217 = 0xbc
0x0218 = 0x63
0x0219 = 0xc6
0x0220 = 0x97
0x0221 = 0x35
0x0222 = 0x6a
0x0223 = 0xd4
0x0224 = 0xb3
0x0225 = 0x7d
0x0226 = 0xfa
0x0227 = 0xef
0x0228 = 0xc5
0x0229 = 0x91
0x0230 = 0x39
0x0231 = 0x72
0x0232 = 0xe4
0x0233 = 0xd3
0x0234 = 0xbd
0x0235 = 0x61
0x0236 = 0xc2
0x0237 = 0x9f
0x0238 = 0x25
0x0239 = 0x4a
0x0240 = 0x94
0x0241 = 0x33
0x0242 = 0x66
0x0243 = 0xcc
0x0244 = 0x83
0x0245 = 0x1d
0x0246 = 0x3a
0x0247 = 0x74
0x0248 = 0xe8
0x0249 = 0xcb
0x0250 = 0x8d
And then we circle back around to…
0x0251 = 0x01

Well, let’s keep looking. What about 0x03?

Orbit(0x03) = { 0x030, 0x031, 0x032, … } = { 0x01, 0x03, 0x05, …, 0xf6 }. Bingo, we hit all 255 non-0x00 elements!

0x030 = 0x01
0x031 = 0x03
0x032 = 0x05
0x033 = 0x0f
0x034 = 0x11
0x035 = 0x33
0x036 = 0x55
0x037 = 0xff
0x038 = 0x1a
0x039 = 0x2e
0x0310 = 0x72
0x0311 = 0x96
0x0312 = 0xa1
0x0313 = 0xf8
0x0314 = 0x13
0x0315 = 0x35
0x0316 = 0x5f
0x0317 = 0xe1
0x0318 = 0x38
0x0319 = 0x48
0x0320 = 0xd8
0x0321 = 0x73
0x0322 = 0x95
0x0323 = 0xa4
0x0324 = 0xf7
0x0325 = 0x02
0x0326 = 0x06
0x0327 = 0x0a
0x0328 = 0x1e
0x0329 = 0x22
0x0330 = 0x66
0x0331 = 0xaa
0x0332 = 0xe5
0x0333 = 0x34
0x0334 = 0x5c
0x0335 = 0xe4
0x0336 = 0x37
0x0337 = 0x59
0x0338 = 0xeb
0x0339 = 0x26
0x0340 = 0x6a
0x0341 = 0xbe
0x0342 = 0xd9
0x0343 = 0x70
0x0344 = 0x90
0x0345 = 0xab
0x0346 = 0xe6
0x0347 = 0x31
0x0348 = 0x53
0x0349 = 0xf5
0x0350 = 0x04
0x0351 = 0x0c
0x0352 = 0x14
0x0353 = 0x3c
0x0354 = 0x44
0x0355 = 0xcc
0x0356 = 0x4f
0x0357 = 0xd1
0x0358 = 0x68
0x0359 = 0xb8
0x0360 = 0xd3
0x0361 = 0x6e
0x0362 = 0xb2
0x0363 = 0xcd
0x0364 = 0x4c
0x0365 = 0xd4
0x0366 = 0x67
0x0367 = 0xa9
0x0368 = 0xe0
0x0369 = 0x3b
0x0370 = 0x4d
0x0371 = 0xd7
0x0372 = 0x62
0x0373 = 0xa6
0x0374 = 0xf1
0x0375 = 0x08
0x0376 = 0x18
0x0377 = 0x28
0x0378 = 0x78
0x0379 = 0x88
0x0380 = 0x83
0x0381 = 0x9e
0x0382 = 0xb9
0x0383 = 0xd0
0x0384 = 0x6b
0x0385 = 0xbd
0x0386 = 0xdc
0x0387 = 0x7f
0x0388 = 0x81
0x0389 = 0x98
0x0390 = 0xb3
0x0391 = 0xce
0x0392 = 0x49
0x0393 = 0xdb
0x0394 = 0x76
0x0395 = 0x9a
0x0396 = 0xb5
0x0397 = 0xc4
0x0398 = 0x57
0x0399 = 0xf9
0x03100 = 0x10
0x03101 = 0x30
0x03102 = 0x50
0x03103 = 0xf0
0x03104 = 0x0b
0x03105 = 0x1d
0x03106 = 0x27
0x03107 = 0x69
0x03108 = 0xbb
0x03109 = 0xd6
0x03110 = 0x61
0x03111 = 0xa3
0x03112 = 0xfe
0x03113 = 0x19
0x03114 = 0x2b
0x03115 = 0x7d
0x03116 = 0x87
0x03117 = 0x92
0x03118 = 0xad
0x03119 = 0xec
0x03120 = 0x2f
0x03121 = 0x71
0x03122 = 0x93
0x03123 = 0xae
0x03124 = 0xe9
0x03125 = 0x20
0x03126 = 0x60
0x03127 = 0xa0
0x03128 = 0xfb
0x03129 = 0x16
0x03130 = 0x3a
0x03131 = 0x4e
0x03132 = 0xd2
0x03133 = 0x6d
0x03134 = 0xb7
0x03135 = 0xc2
0x03136 = 0x5d
0x03137 = 0xe7
0x03138 = 0x32
0x03139 = 0x56
0x03140 = 0xfa
0x03141 = 0x15
0x03142 = 0x3f
0x03143 = 0x41
0x03144 = 0xc3
0x03145 = 0x5e
0x03146 = 0xe2
0x03147 = 0x3d
0x03148 = 0x47
0x03149 = 0xc9
0x03150 = 0x40
0x03151 = 0xc0
0x03152 = 0x5b
0x03153 = 0xed
0x03154 = 0x2c
0x03155 = 0x74
0x03156 = 0x9c
0x03157 = 0xbf
0x03158 = 0xda
0x03159 = 0x75
0x03160 = 0x9f
0x03161 = 0xba
0x03162 = 0xd5
0x03163 = 0x64
0x03164 = 0xac
0x03165 = 0xef
0x03166 = 0x2a
0x03167 = 0x7e
0x03168 = 0x82
0x03169 = 0x9d
0x03170 = 0xbc
0x03171 = 0xdf
0x03172 = 0x7a
0x03173 = 0x8e
0x03174 = 0x89
0x03175 = 0x80
0x03176 = 0x9b
0x03177 = 0xb6
0x03178 = 0xc1
0x03179 = 0x58
0x03180 = 0xe8
0x03181 = 0x23
0x03182 = 0x65
0x03183 = 0xaf
0x03184 = 0xea
0x03185 = 0x25
0x03186 = 0x6f
0x03187 = 0xb1
0x03188 = 0xc8
0x03189 = 0x43
0x03190 = 0xc5
0x03191 = 0x54
0x03192 = 0xfc
0x03193 = 0x1f
0x03194 = 0x21
0x03195 = 0x63
0x03196 = 0xa5
0x03197 = 0xf4
0x03198 = 0x07
0x03199 = 0x09
0x03200 = 0x1b
0x03201 = 0x2d
0x03202 = 0x77
0x03203 = 0x99
0x03204 = 0xb0
0x03205 = 0xcb
0x03206 = 0x46
0x03207 = 0xca
0x03208 = 0x45
0x03209 = 0xcf
0x03210 = 0x4a
0x03211 = 0xde
0x03212 = 0x79
0x03213 = 0x8b
0x03214 = 0x86
0x03215 = 0x91
0x03216 = 0xa8
0x03217 = 0xe3
0x03218 = 0x3e
0x03219 = 0x42
0x03220 = 0xc6
0x03221 = 0x51
0x03222 = 0xf3
0x03223 = 0x0e
0x03224 = 0x12
0x03225 = 0x36
0x03226 = 0x5a
0x03227 = 0xee
0x03228 = 0x29
0x03229 = 0x7b
0x03230 = 0x8d
0x03231 = 0x8c
0x03232 = 0x8f
0x03233 = 0x8a
0x03234 = 0x85
0x03235 = 0x94
0x03236 = 0xa7
0x03237 = 0xf2
0x03238 = 0x0d
0x03239 = 0x17
0x03240 = 0x39
0x03241 = 0x4b
0x03242 = 0xdd
0x03243 = 0x7c
0x03244 = 0x84
0x03245 = 0x97
0x03246 = 0xa2
0x03247 = 0xfd
0x03248 = 0x1c
0x03249 = 0x24
0x03250 = 0x6c
0x03251 = 0xb4
0x03252 = 0xc7
0x03253 = 0x52
0x03254 = 0xf6
… and only now do we circle around to…
0x03255 = 0x01

This gives us a one-to-one mapping between the 255 powers 0 through 254 and the 255 non-zero bytes 0x01 through 0xff.

Here are the perl scripts I used to generate these:

use strict;

my $m = 0x11b;

my $a = 1;
for (my $i = 0; ; $i++) {
    printf("0x02^%2d = 0x%02x\n", $i, $a);
    if ($i > 0 and $a == 1) { last; }
    $a <<= 1; # * 0x02
    if ($a > 0xff) { $a ^= $m; }
}

print "\n";

$a = 1;
for (my $i = 0; ; $i++) {
    printf("0x03^%3d = 0x%02x\n", $i, $a);
    if ($i > 0 and $a == 1) { last; }
    $a = ($a << 1) ^ $a; # * 0x03
    if ($a > 0xff) { $a ^= $m; }
}

Now that we have this orbit, we use it to construct two lookup tables:

// 255 entries xPlusOne_ToThe[0] = 0x01 to xPlusOne_ToThe[254] = 0xf6
xPlusOne_ToThe[] = { 0x01, 0x03, 0x05, ..., 0x52, 0xf6 };

And the inverse mapping:

// 255 entries logXPlusOne_Of[0x01] = 0 to logXPlusOne_Of[0xff] = 7
logXPlusOne_Of[] = { UNDEFINED, 0, 25, 1, 50, ..., 112, 7 };

Now that we have paid for these 512-ish bytes, we can do multiplication very cheaply:

multiply(b1, b2) {
    if (0x00 == b1 || 0x00 == b2) { return 0x00; }

    logAns = logXPlusOne_Of[b1] + logXPlusOne_Of[b2];
    if (logAns >= 255) { logAns -= 255; }
    return xPlusOne_ToThe[logAns];
}

Sieving irreducible monic polynomials over a finite field

Last time we talked about the addition and multiplication tables for GF(2²); GF(28) is relevant for the Rijndael encryption scheme.

Joan Daemen and Vincent Rijmen use a representation of GF(28) where each element is a polynomial of the form b7x7 + b6x6 + b5x5 + b4x4 + b3x3 + b2x2 + b1x + b0, with the bi being bits. A polynomial is represented by a single byte. Addition is component-wise (modulo 2); two polynomials can be added with a single XOR operation. Multiplication is a little more complicated.

Rijndael defines multiplication as normal polynomial multiplication, followed by taking a modulus using the reduction polynomial mx8 + x4 + x3 + x + 1. Last time we showed that, in GF(2²), x2 + x + 1 was used to reduce polynomial multiplication, and also that a necessary quality for a reduction polynomial is that it be irreducible/prime.

Last time I hinted at the question: how do we show that m is irreducible? Well, one way is to do trial division of all polynomials up to degree 4.

But that’s no fun. Instead, let’s write a script to generate all irreducible polynomials, and see if m is on it! The approach is very similar to the Sieve of Eratosthenes: generate a list of all the polynomials, then traverse it from the low end to the high end; circle each element that hasn’t been crossed out, then cross out all multiples of that element.

The first argument q is the modulus of the coefficients (in this case, 2) and the second argument d (in this case, 8) is how high the degree can go. Here is the output of the script:

>perl polynomial-sieve.pl 2 8
Finding monic irreducible polynomials of degree up to 8 and coefficients mod 2
Generating all monic polynomials…
Sieving out all reducible polynomials…
1
x
x + 1
x^2 + x + 1
x^3 + x + 1
x^3 + x^2 + 1
x^4 + x + 1
x^4 + x^3 + 1
x^4 + x^3 + x^2 + x + 1
x^5 + x^2 + 1
x^5 + x^3 + 1
x^5 + x^3 + x^2 + x + 1
x^5 + x^4 + x^2 + x + 1
x^5 + x^4 + x^3 + x + 1
x^5 + x^4 + x^3 + x^2 + 1
x^6 + x + 1
x^6 + x^3 + 1
x^6 + x^4 + x^2 + x + 1
x^6 + x^4 + x^3 + x + 1
x^6 + x^5 + 1
x^6 + x^5 + x^2 + x + 1
x^6 + x^5 + x^3 + x^2 + 1
x^6 + x^5 + x^4 + x + 1
x^6 + x^5 + x^4 + x^2 + 1
x^7 + x + 1
x^7 + x^3 + 1
x^7 + x^3 + x^2 + x + 1
x^7 + x^4 + 1
x^7 + x^4 + x^3 + x^2 + 1
x^7 + x^5 + x^2 + x + 1
x^7 + x^5 + x^3 + x + 1
x^7 + x^5 + x^4 + x^3 + 1
x^7 + x^5 + x^4 + x^3 + x^2 + x + 1
x^7 + x^6 + 1
x^7 + x^6 + x^3 + x + 1
x^7 + x^6 + x^4 + x + 1
x^7 + x^6 + x^4 + x^2 + 1
x^7 + x^6 + x^5 + x^2 + 1
x^7 + x^6 + x^5 + x^3 + x^2 + x + 1
x^7 + x^6 + x^5 + x^4 + 1
x^7 + x^6 + x^5 + x^4 + x^2 + x + 1
x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + 1
x^8 + x^4 + x^3 + x + 1
x^8 + x^4 + x^3 + x^2 + 1
x^8 + x^5 + x^3 + x + 1
x^8 + x^5 + x^3 + x^2 + 1
x^8 + x^5 + x^4 + x^3 + 1
x^8 + x^5 + x^4 + x^3 + x^2 + x + 1
x^8 + x^6 + x^3 + x^2 + 1
x^8 + x^6 + x^4 + x^3 + x^2 + x + 1
x^8 + x^6 + x^5 + x + 1
x^8 + x^6 + x^5 + x^2 + 1
x^8 + x^6 + x^5 + x^3 + 1
x^8 + x^6 + x^5 + x^4 + 1
x^8 + x^6 + x^5 + x^4 + x^2 + x + 1
x^8 + x^6 + x^5 + x^4 + x^3 + x + 1
x^8 + x^7 + x^2 + x + 1
x^8 + x^7 + x^3 + x + 1
x^8 + x^7 + x^3 + x^2 + 1
x^8 + x^7 + x^4 + x^3 + x^2 + x + 1
x^8 + x^7 + x^5 + x + 1
x^8 + x^7 + x^5 + x^3 + 1
x^8 + x^7 + x^5 + x^4 + 1
x^8 + x^7 + x^5 + x^4 + x^3 + x^2 + 1
x^8 + x^7 + x^6 + x + 1
x^8 + x^7 + x^6 + x^3 + x^2 + x + 1
x^8 + x^7 + x^6 + x^4 + x^2 + x + 1
x^8 + x^7 + x^6 + x^4 + x^3 + x^2 + 1
x^8 + x^7 + x^6 + x^5 + x^2 + x + 1
x^8 + x^7 + x^6 + x^5 + x^4 + x + 1
x^8 + x^7 + x^6 + x^5 + x^4 + x^2 + 1
x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + 1
Done!

Note that m is not just prime, it is the lowest of the monic 8-degree polynomials with coefficients mod 2.

The script itself (in Perl) is attached. It makes no claims to being pretty or efficient.

As a check, there is an elegant formula for the number of irreducible monic polynomials with coefficients in a finite field:

N(q, n) = (Σd|n μ(d) qn/d)/n

where μ(x) is the Möbius function.

In particular:

N(2, 1) = ((1)21/1)/1 = 2
N(2, 2) = ((1)22/1 – (1)22/2)/2 = 1
N(2, 3) = ((1)23/1 – (1)23/3)/3 = 2
N(2, 4) = ((1)24/1 – (1)24/2 + (0)24/4)/4 = 3
N(2, 5) = ((1)25/1 – (1)25/5)/5 = 6
N(2, 6) = ((1)26/1 – (1)26/2 – (1)26/3 + (1)26/6)/6 = 9
N(2, 7) = ((1)27/1 – (1)27/7)/7 = 18
N(2, 8) = ((1)28/1 – (1)28/2 + (0)28/4 + (0)28/8)/8 = 30

This checks out with the script.

EDIT 2015-10-31: move script to https://github.com/mvaneerde/blog/blob/master/rijndael/polynomial-sieve.pl