Combinatorics

Cryptographic Keys and Decks of Cards

Posted on by John

The previous post looked at the idea of storing a cryptographic key in the order of a deck of cards. A deck of 52 cards can store 225 bits of data because

⌊log2(52!)⌋ = 225.

Here ⌊x⌋ is x rounded down to the nearest integer.

If we want to store bigger keys, we’re going to need a bigger deck of cards.

Bitcoin

A Bitcoin key has 256 bits, which would require a deck of 58 cards. There is a card game called Zwicker that uses a deck of 58 cards, the usual 52 cards plus six jokers. So you could store a Bitcoin key in the permutation of a Zwicker deck.

You could also use a deck of 52 cards, plus 2 jokers, if you also consider orientation. 30 cards are rotationally symmetric, 22 are not, and neither are jokers. So, including two asymmetric jokers, you could add 24 additional bits. Permutations of a 54 card deck can encode 237 bits, and with 24 orientation bits, this is a total of 261 bits.

RSA

RSA key sizes vary, but 2048 and 3072 are common. A 2048-bit key would require a deck of 301 cards. Casinos often use a shoe of 312 cards, combining six decks of 52 cards, to deal Baccarat or Blackjack. However, casinos combine identical decks. If you were to combine six unique decks, you could store a 2048-bit key.

Storing a 3072-bit key would require a deck of 422 cards. You could make a deck of 432 cards by combining 8 distinguishable packs of 54 cards (52 + 2 jokers).

ML-KEM

ML-KEM is a proposed quantum-resistant replacement for RSA. As with RSA, key sizes for ML-KEM vary, the smallest being ML-KEM-512 with a key size of 1632 bytes, which equals 13056 bits. This would require a deck of 1442 cards. You could combine 28 distinct packs of 52 cards, but that’s unwieldy.

This illustrates one of the difficult trade-offs with post-quantum cryptography: key sizes are much bigger. If you wanted to create a deck of 1442 cards, you’d probably want to make your “cards” something other than standard playing cards. You’d want to use permutations of something else.

Verification

The following Python code verifies the calculations above.

from math import log2, factorial, floor

def capacity(cards):
    return floor(log2(factorial(cards)))

def verify(bits, cards):
    return capacity(cards) >= bits and capacity(cards-1) < bits

print(verify(237, 54))
print(verify(256, 58))
print(verify(2048, 301))
print(verify(3072, 422))
print(verify(1632*8, 1442))

For more on how I came up with the deck sizes, see the next post on computing the inverse factorial.

Hiding data in permutations

Posted on by John

The latest issue of Paged Out! has an article by Stephen Hewitt “An off-line backup of your cryptographic key using playing cards.” The idea is to use a deck of 52 to store a 128-bit cryptographic key. To erase the key, shuffle the deck. Hewitt gives his algorithm for embedding a key, one that can be carried out manually but isn’t maximally efficient.

You could store a 225-bit key as a permutation of 52 cards because

log2(52!) = 225.581.

But then how would you number permutations so you could go from a number to a particular permutation and later decode the permutation to a number? Is this even practical? For a small number n, you could encode a number k < n by enumerating the first k permutations of a set of n items, and you could decode by enumerating permutations until you find the one you have. But this is completely impractical for large n, such as n = 52.

The process of mapping permutation to an integer is called ranking, and the mapping from an integer to a permutation is called unranking. How efficiently can rankings and unrankings be calculated?

Let n be the number of symbols being permuted. Then there are simple algorithms for ranking and unranking with respect to lexicographical order that have complexity O(n²) and more sophisticated algorithms that have complexity O(n log n). There are also O(n) algorithms that do not preserve lexicographical order.

The Permutations class in SymPy has methods unrank_lex and rank to unrank and rank permutations according to lexicographical order.

The notation the Permutations class uses requires a little explanation. For example, suppose we unrank 2026.

>>> from sympy.combinatorics import Permutation
>>> Permutation.unrank_lex(52, 2026)
Permutation(45, 47, 51, 48, 46, 50)

The output is not a full list of 52 numbers in permuted order; it is only a cycle. The notation refers to the permutation that sends 45 to 47, 47 to 51, …, 50 to 45 and leaves everything else fixed.

If we rank the permutation given above, we get 2026 back.

>>> Permutation.rank(Permutation(45, 47, 51, 48, 46, 50))
2026

Note that we didn’t say how many elements (45, 47, 51, 48, 46, 50) is a permutation of. Because of lexicographical order, the rank would be the same whether we viewed this as a permutation of 52 objects or of more objects.

Now let’s do something larger. Let’s generate a 220-bit number and encode it as a permutation.

>>> n = random.getrandbits(225)
>>> a = Permutation.unrank_lex(52, n)
>>> n
40234719030664563684489051530416964877785781669439875437823431388841
>>> a
Permutation(0, 25, 32, 15, 8, 28)(1, 48, 34, 14, 10, 51, 38, 31, 21, 5, 42, 47, 29, 26, 46, 30, 50, 49, 37, 22, 18, 23)(2, 45, 17, 20, 36, 40, 11, 4, 7, 41, 33, 3, 43, 44, 19, 16, 35, 39, 12, 6, 9)
>>> Permutation.rank(a) == n
True

Now just for fun, let’s display the permutation above applied to a standard (French) deck of 52 cards. As explained here, symbols associated with these cards have a range of Unicode values. By printing these values, we can visualize the permuted deck.

Here’s the code that made the image above.

spades = list(range(0x1F0A1, 0x1F0AF))
spades.remove(0x1F0AC) # take out the knight
cards = [s + 16*i for s in spades for i in range(4)]

a = Permutation.unrank_lex(52, n)
p = a(cards)

for i in range(4):
    for j in range(13):
        print(chr(p[13*i + j]), end="")
    print()

The code above is plenty fast, but Permutation has methods rank_nonlex and unrank_nonlex that run in O(n) time, which could be useful for n much larger than 52.

Counting permutations with roots

Posted on by John

My post from yesterday on permutation roots ends with a Mathematica code for finding the probability that a permutation of n elements has a kth root. This is done by finding the coefficient of xn in the generating function

\prod_{m=1}^\infty \exp_{\text{gcd}(m, k)} \left\frac{x^m}{m}\right)

I wanted to say more about this, and look at implementing the same code in SymPy. I was curious how well SymPy would do because I’ve noticed that LLMs often generate SymPy code since it’s an open source CAS.

Wilf [1] describes the infinite product above as the exponential generating function (egf) of f(n, k), the number of permutations of n objects that have a kth root. Since egfs have a n! term in the denominator, this is also the ordinary generating function (ogf) of the probability that a randomly chosen permutation on n objects has a kth root.

My first attempt at using Mathematica to probe the generating function was

expq[x_, q_] := MittagLefflerE[q, x^q]	 
p[n_, k_] :=  SeriesCoefficient[	 
    Product[expq[x^m/m, GCD[m, k]], {m, 1, Infinity}], {x, 0, n}]

This hung forever when I tried to use it on a small example. I realized, but apparently Mathematica did not, that Infinity could be replaced by n since terms higher than n do not contribute to the coefficient of xn. With that change, the code ran quickly.

This morning I tried converting the Mathematica code to Sympy; Claude did this in one shot. I also reproduced the table of f(n, k) values on page 150 of [1] to test the code. Since Wilf tabulated f(n, k), not f(n, k)/n!, I multiplied the results by n!.

Here is the output:

k = 2 [1, 1, 3, 12, 60, 270, 1890, 14280, 128520, 1096200]
k = 3 [1, 2, 4, 16, 80, 400, 2800, 22400, 181440, 1814400]
k = 4 [1, 1, 3, 12, 60, 270, 1890, 13020, 117180, 1039500]
k = 5 [1, 2, 6, 24, 96, 576, 4032, 32256, 290304, 2612736]
k = 6 [1, 1, 1, 4, 40, 190, 1330, 8680, 52920, 340200]
k = 7 [1, 2, 6, 24, 120, 720, 4320, 34560, 311040, 3110400]

and here is the SymPy code. I edited the main but the rest is verbatim from Claude.

from sympy import symbols, gcd, factorial, Rational, S

x = symbols('x')

def expq_coeffs(m, q, n):
    """
    Truncated (degree <= n) series coefficients of
        expq(x**m/m, q) = MittagLefflerE(q, (x**m/m)**q)
    Since q is a positive integer:
        E_q(y^q) = sum_j y^(q*j) / (q*j)!
    with y = x**m/m, so the term of degree m*q*j has coefficient
        1 / ( m**(q*j) * (q*j)! ).
    Returns a list c[0..n] of coefficients.
    """
    c = [S.Zero] * (n + 1)
    j = 0
    while m * q * j <= n:
        deg = m * q * j
        c[deg] += Rational(1, m**(q * j) * factorial(q * j))
        j += 1
    return c

def poly_mult_trunc(a, b, n):
    """Multiply two series (lists of coeffs, index = degree) truncated to degree n."""
    c = [S.Zero] * (n + 1)
    for i, ai in enumerate(a):
        if ai == 0:
            continue
        max_j = n - i
        for j2 in range(max_j + 1):
            bj = b[j2]
            if bj != 0:
                c[i + j2] += ai * bj
    return c

def p(n, k):
    """
    SymPy equivalent of:
        expq[x_, q_] := MittagLefflerE[q, x^q]
        p[n_, k_] := SeriesCoefficient[
            Product[expq[x^m/m, GCD[m, k]], {m, 1, n}], {x, 0, n}]
    """
    result = [S.Zero] * (n + 1)
    result[0] = S.One
    for m in range(1, n + 1):
        q = gcd(m, k)
        factor = expq_coeffs(m, q, n)
        result = poly_mult_trunc(result, factor, n)
    return result[n]

# example
if __name__ == "__main__":
    for k in range(2, 8):
        print("k =", k, [factorial(n)*p(n, k) for n in range(1,11)])

[1] Herbert Wilf. Generatingfunctionology. Available online here.

DNA Sequence Alignment and Kings

Posted on by John

This morning I wrote a post that included the central Delannoy numbers. The nth central Delannoy number Dn counts the number of ways a king can move from one corner of a chessboard to the diagonally opposite corner without backtracking.

The more general Delannoy numbers Dm,n are the analogy for an m × n rectangular board, not necessarily square.

Dm,n is also the number of possible sequence alignments for a strand of DNA with m base pairs and a strand with n base pairs [1]. At each step in the alignment process, you can introduce a gap in the first strand, the second strand or neither, which is analogous to the king who can move N, E, or NE at each step.

The Delannoy numbers can be computed recursively:

def D(m, n):
    if m == 0 or n == 0:
        return 1
    return D(m - 1, n) + D(m, n - 1) + D(m - 1, n - 1)

The code above can be sped up tremendously by adding the decorator

@lru_cache(maxsize=None)

above the function definition to turn on memoization. I did an experiment computing D12,15 with and without memoization and the times were 77.1805 seconds and 0.000062 seconds respectively, i.e. memoization made the code over a million times faster.

Incidentally, D12,15 = 2653649025 and so there are a lot of ways to align even short sequences unless you place some restriction on the permissible alignments.

Update: Here’s a heatmap plotting log10(Dm,n). Obviously the function increases with m and n: bigger chessboards have more possible paths. Moreover, it’s larger along the diagonal (i.e. the central Delannoy numbers). If you look along northeast to southwest diagonals, the function is largest in the middle where m = n.

[1] Torres, A., Cabada, A., & Nieto, J. J. (2003). An exact formula for the number of alignments between two DNA sequences. DNA Sequence, 14(6), 427–430. https://doi.org/10.1080/10425170310001617894

Silver Rectangles and the Ways of Kings

Posted on by John

Golden rectangles

The defining property of golden rectangle is that if you stick a square on its longer side, you get another golden rectangle.

The smaller vertical rectangle is similar to the larger horizontal rectangle. This means

φ / 1 = (1 + φ) / φ

which tells us φ² = 1 + φ and so the golden ratio φ equals (1 + √5)/2.

Silver rectangles

A silver rectangle is one that if you stick two squares on its longer side you get another rectangle with the same aspect ratio.

This tells us

σ / 1 = (1 + 2σ) / σ

and so σ² = 1 + 2σ and the silver ratio is σ = 1 + √2.

Just as you can define a golden ratio and a silver ratio, there’s an analogous way to define a sequence of metallic ratios.

Kings and Delannoy numbers

The silver ratio has several connections to the ways of ways kings. By that I mean the number of ways a king can go from one corner of a chessboard to the diagonally opposite corner without backtracking.

A king can move one space in any direction. If we start with a king in the bottom left corner of the board, the no-backtracking requirement means the king can move up, right, or up and right.

The number of paths a king can take from one corner to the opposite corner of an n × n chessboard is the nth central Delannoy number Dn. more generally Delannoy numbers are defined for an m × n chessboard, but I’ll stick to the case m = n called the central Delannoy number, or just Delannoy numbers for short.

The first Delannoy number is 1 because there’s only one way for a king to get from one corner to the other: do nothing, because the opposite corner is the same corner. The second Delannoy number is 3 because the king can move up then right, or right then up, or move diagonally up and right.

For a 3 × 3 grid things are significantly more complicated, and D3 = 13. For an 8 × 8 grid the number of paths is 48,639.

Generating function

How would you estimate the number of paths on an n × n board for large values of n without calculating it exactly? You might start by finding a generating function for the Delannoy numbers, which works out to be

(x² − 6x + 1)−1/2

The radius of convergence r for the generating function series is the distance from 0 to the closest singularity of the generating function, which is the smaller root of

x² − 6x + 1

which is

3 − √8 = (3 + √8)−1 = (1 + √2)−2 = 1/σ²

i.e. the radius of convergence is the reciprocal of the silver ratio squared.

Asymptotic estimate

The radius of convergence gives us a first approximation to the asymptotic size of the series coefficients. Since we’re working with the generating function of the Delannoy numbers, these coefficients are the Delannoy numbers. That is,

Dn ~ r−n = (σ2)n = σ2n.

That’s as good as you can do just knowing the radius of convergence. A more careful analysis would refine this estimate by dividing by a factor proportional to √n.

Related posts

Queens on a prime order board

Posted on by John

The n queens problem is to place on an n × n chessboard n queens so that none attacks any other. This means there is only one queen on every horizontal, vertical, and diagonal line.

When n is a prime number ≥ 5, it is sufficient to place the queens on a line that has slope 2, 3, 4, …, n − 2. (The slope cannot be 1 because that’s a diagonal. And it cannot be n − 1 because n − 1 = −1 mod n is also a diagonal.) [1]

Here we imagine opposite edges of the board being joined together. Geometrically, this makes the chessboard a torus (donut). Algebraically, the points on a line of slope s have the coordinates

(a + k, b + ks)

where addition is carried out mod n.

All solutions to the n queens problem have this form when n = 5. Some solutions will have this form for larger prime values of n but not all.

For example, when n = 7, here is a solution where all the queens are on a line of slope 2.

But here is another solution where the queens do not all lie on a line of constant slope.

Related posts

[1] W. H. Bussey. A Note on the Problem of the Eight Queens. The American Mathematical Monthly, Vol. 29, No. 7 (August 1922), pp. 252–253

Partitions over permutations

Posted on by John

I was thinking more about the cosine approximation to the Gaussian

exp(−z²) ≈ (1 + cos(sin(z) + z))/2

that I wrote about last week. The two expressions above are close along the real axis but not along the imaginary axis.

If z = iy, the right side grows much faster than the left, behaving like exp(exp(y)).

This led to me looking up the power series for the double exponential function exp(exp(y)). This is an interesting series because the coefficient of xn is

e Bn / n!

where Bn is the nth Bell number, which equals the number of ways to partition a set of n labeled items [1]. And of course n! is the number of ways to permute a set of n labeled items. So the nth coefficient in the power series for exp(exp(y)) is the ratio of the number of partitions to permutations for a set of n labeled things, multiplied by e.

The number of ways to partition a set of n things grows quickly as n increases, almost as quickly as the number of permutations, and so the series for the double exponential function converges very slowly.

Computing

SymPy has a function bell for computing Bell numbers, so you could compute the ratio of partitions to permutations as follows.

from sympy import bell, factorial
f = lambda n: bell(n)/factorial(n)

This returns a number of type sympy.core.numbers.Rational and so the result is exact. You can cast it to float for convenience.

Asymptotics

If we look at only the terms in the asymptotic series for log Bn and log n! that grow with n we have

log Bn ~ n log n − n log log n
log n! ~ n log n − ½ log n

and so

log( Bn / n! ) ~ ½ log n − n log log n

There’s also an asymptotic series for log( Bn / n! ) involving the Lambert W function:

log( Bn / n! ) ~ n/r − 1 − n log r

where r = W(n).

Related posts

[1] It’s important that the items are labeled. Partition numbers are the number of partitions of an unlabeled set. Partition numbers are much smaller than Bell numbers.

Tone row operations

Posted on by John

The previous post introduced the idea of a twelve-tone row, a permutation of the twelve pitch classes of a chromatic scale.

The earlier post also introduced a group of operations on a tone row with elements P, R, I, and RI. Here P, which stands for “prime”, is the identity operator: it leaves the tone row alone. R stands for retrograde, which means to reverse the sequence. I stands for inversion, inverting each of the intervals in the row. If you apply R then I, you get the same result as if you first applied I then R. See the previous post for an example of each.

Do these operations always return a new tone row? In other words, are P(t), R(t), I(t), and RI(t) always distinct for all tone rows t?

It’s easy to see that P(t) and R(t) are always different. The first and last elements of t are different, and so the first element of P(t) does not equal the first element of R(t).

There is one interval that remains the same when inverted, namely the tritone. It’s possible to create a two-tone row that does not change under inversion, such as the sequence [B, F]. A tone row with more than two notes must have an interval that is not a tritone, and so inverting the tone row changes it.

If a tone row t has at least three notes, then P(t), R(t), I(t), and RI(t) are all distinct.

The pitch classes in a tone row are repeated in a cycle, so we might consider all rotations of a tone row to be in a sense the same. In mathematical terms, we can mod out by cyclic permutations.

If we apply R, I, and RI to equivalence classes of tone rows, the results are not always unique. For example, let t be the chromatic scale.

C C♯ D D♯ E F F♯ G G♯ A A♯ B

Then R(t) is

B A♯ A G♯ G F♯ F E D♯ D C♯ C

and I(t) is

C B A♯ A G♯ G F♯ F E D♯ D C♯

which is a rotation of R(t).

The chromatic scale is a boring tone row. For a random tone row t, all the variations P(t), R(t), I(t), and RI(t) will very likely be distinct, including rotations. The following Python code will illustrate the points above and estimate the probability that manipulations of a tone row are distinct even when equating rotations.

import random
import itertools

def inversion(x, m = 12):
    return [(2*x[0] - x[i]) % m for i in range(m)]

def rotate_left(x, n):
    return x[n:] + x[:n]

def rotdiff(x, y):
    for i in range(len(x)):
        if rotate_left(x, i) == y:
            return False
    return True

c = 0
for _ in range(1_000_000):
    P = list(range(12))
    random.shuffle(P) # shuffle in place
    R = list(reversed(P))
    I = list(inversion(P))
    RI = list(reversed(I))
    for pair in combinations([P, R, I, RI], 2):
        assert(pair[0] != pair[1])
        if not rotdiff(pair[0], pair[1]):
            c += 1
print(c)

This generates and tests a million random tone rows. When I ran it, it found 226 tone rows were one manipulation of the tone row was equal to a rotation of another. So when I said the variations are “very likely” to be distinct, you could quantify that as having over a 99.9% chance.

Combining in-shuffles and out-shuffles

Posted on by John

A few days ago I wrote two posts about perfect shuffles. Once you’ve cut a deck of cards in half, an in-shuffle lets a card from the top half fall first, and an out-shuffle lets a card from the bottom half fall first.

Suppose we have a deck of 52 cards. We said in the earlier posts that the order of an in-shuffle I is 52. That is, after 52 in-shuffles, a deck returns to its initial order. And the order of an out-shuffle O is 8.

We can think of I and O as generators of subgroups of order 52 and 8 respectively in the group S of all permutations of 52 cards. I was curious when I wrote the earlier posts how large the group generated by I and O together would be. Is it possible to reach all 52! permutations of the deck by some combination of applying I and O? If not, how many permutations can be generated?

I’ve since found the answer in [1] in a theorem by Diaconis, Graham, and Kantor. I don’t know who Kantor is, but it’s no surprise that a theorem on card shuffles would come from Persi Diaconis and Ron Graham. The theorem covers the case for decks of size N = 2n, which branches into different results depending on the size of n and the value of n mod 4.

For N = 52, the group generated by I and O has

26! × 226

elements.

On the one hand, that’s a big number, approximately 2.7 × 1034. On the other hand, it’s quite small compared to 52! = 8 × 1067. So while there are a lot of permutations reachable by a combination of in-shuffles and out-shuffles, your chances of selecting such a permutation from the set of all such permutations is vanishingly small.

To put it yet another way, the number of arrangements is on the order of the square root of 52!, a big number, but not big relative to 52!. (Does this pattern

√52! ≈ 26! × 226

generalize? See the next post.)

Not only does the theorem of Diaconis et al give the order of the group, it gives the group itself: the group of permutations generated by I and O is isomorphic to the group of symmetries of a 26-dimensional octahedron.

[1] S. Brent Morris. Magic Tricks, Card Shuffling and Dynamic Computer Memories. MAA 1998.

In-shuffles and out-shuffles

Posted on by John

The previous post talked about doing perfect shuffles: divide a deck in half, and alternately let one card from each half fall.

It matters which half lets a card fall first. If the top half’s bottom card falls first, this is called an in-shuffle. If the bottom half’s bottom card falls first, it’s called an out-shuffle.

With an out-shuffle, the top and bottom cards don’t move. Presumably it’s called an out-shuffle because the outside cards remain in place.

An out-shuffle amounts to an in-shuffle of the inner cards, i.e. the rest of the deck not including the top and bottom card.

The previous post had a Python function for doing an in-shuffle. Here we generalize the function to do either an in-shuffle or an out-shuffle. We also get rid of the list comprehension, making the code longer but easier to understand.

def shuffle2(deck, inside = True):
    n = len(deck)
    top = deck[: n//2]
    bottom = deck[n//2 :]
    if inside:
        first, second = bottom, top
    else:
        first, second = top, bottom
    newdeck = []
    for p in zip(first, second):
        newdeck.extend(p)
    return newdeck

Let’s use this code to demonstrate that an out-shuffle amounts to an in-shuffle of the inner cards.

deck = list(range(10))
d1 = shuffle2(deck, False) 
d2 = [deck[0]] + shuffle2(deck[1:9], True) + [deck[9]]
print(d1)
print(d2)

Both print statements produce [0, 5, 1, 6, 2, 7, 3, 8, 4, 9].

I said in the previous post that k perfect in-shuffles will restore the order of a deck of n cards if

2k = 1 (mod n + 1).

It follows that k perfect out-shuffles will restore the order of a deck of n cards if

2k = 1 (mod n − 1)

since an out-shuffle of n cards is essentially an in-shuffle of the n − 2 cards in the middle.

So, for example, it only takes 8 out-shuffles to return a deck of 52 cards to its original order. In the previous post we said it takes 52 in-shuffles, so it takes a lot fewer out-shuffles than in-shuffles.

It’s plausible to conjecture that it takes fewer out-shuffles than in-shuffles to return a deck to its initial order, since the former leaves the two outside cards in place. But that’s not always true. It’s true for a deck of 52 cards, but not for a deck of 14, for example. For a deck of 14 cards, it takes 4 in-shuffles or 12 out-shuffles to restore the deck.