John D. CookModify

2026-07-28 22:11:30

A couple years ago I wrote about how to compute the inverse of factorial. I used that code in writing the previous post because the post required solving the equation

⌊log₂(n!)⌋ ≥ b

given b. That is, given a number of bits b, find the smallest value of n such that n! ≥ 2^b.

What the code got right

Looking back on the code in that post, there are a few changes I’d like to make. But first of all, I’d like to point out something the post does right: instead of trying to solve

Γ(y) = x

it solves

log Γ(y) = log x.

That’s why the argument to inverse_log_gamma is logarg. That makes the code useful for values of x that would far exceed the maximum floating point value, such as in the calculations for the previous post.

What I’d change

Rounding

The function inverse_factorial from the old post solves finds the closest integer solution. It would be better for it to return the solution without rounding and then let the user round result if they want to. In my calculations in the previous post, I wanted to take the floor, not round.

Newton’s method

The code in the previous post uses the bisection method. This method is very safe, and fast enough for my purposes, but it could be made faster. Newton’s method is faster, but it can be ill-behaved if you don’t start close enough to the solution.

It’s safe to use Newton’s method to invert log Γ for two reasons. First, you can get a good starting point based on Stirling’s approximation. Second, and more importantly, log Γ is convex. Newton’s method will converge from any starting point when applied to a convex function. A little caution is necessary because log Γ is not convex everywhere, but it is convex on the positive real axis.

Another difficulty with Newton’s method is that you need to supply the derivative of the function whose root you’re trying to find. But this isn’t an issue here because the derivative of log Γ is the digamma function, which is implemented in SciPy.

Tolerance

Finally, the previous code used the default tolerance for deciding when to stop refining the solution. The revised method lets the user specify tolerance. It provides a default value, but that default is visible in the function call, not hidden down in SciPy.

Revised code

Here’s the revised code.

from scipy.special import gammaln, digamma
from scipy.optimize import newton

def inverse_log_gamma(logarg, tol=1e-12):
    assert(logarg > 0)    
    x0 = logarg / log(logarg + 1) + 1 if logarg > 1 else 2.0
    def f(z): return gammaln(z) - logarg
    return newton(f, x0, fprime=digamma, tol=tol)

def inverse_factorial(logarg):
    g = inverse_log_gamma(logarg)
    return g - 1 

2026-07-28 20:13:17

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

⌊log₂(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.

2026-07-28 07:08:59

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

log₂(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.

2026-07-28 00:03:44

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 xⁿ in the generating function

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 xⁿ. 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.

2026-07-27 23:34:05

It’s well known that you can convert the base 16 (hex) representation of an integer to the base 2 (binary) representation by simply converting each digit from hex to binary. For example,

CAFE_hex = 1100 1010 1111 1110_two

I imagine it’s less well known that you can do the same thing with floating point numbers.

I wanted to find the binary representation of a floating point number using Python, and discovered that it has no function to do this. However, there is a method on floats to show a hex representation. For example, here’s the hex representation of π.

>>> import math
>>> (math.pi).hex()
'0x1.921fb54442d18p+1'

Curiously, the p+k part at the end is an exponent of 2, not an exponent of 16. So after we convert 1.921fb54442d18 to binary, we’ll need to multiply by 2, i.e. move the fractional point one space to the right.

So first we convert 1.921fb54442d18_hex to binary by converting 1, 9, 2, etc. each to binary.

1.1001 0010 0001 1111 1011 0101 0100 0100 0100 0010 1101 0001 1000_two

Then after shifting the fraction point to account for the p+1 part we have

π = 11.001001000011111101101010100010001000010110100011000_two

You could use Python’s bin() function to convert the fractional part, interpreted as an integer, to hex, though you may need to pad with 0 bits. For example,

>>> (1.03).hex()
'0x1.07ae147ae147bp+0

>>> bin(0x7ae147ae147)
'0b1111010111000010100011110101110000101000111'

The binary representation of 1.03_ten is

1.000001111010111000010100011110101110000101000111_two

We added a total of five zero bits, four for the 0 after the fractional point and one for converting 7 to 0111_two.

2026-07-27 04:32:02

Let σ be a permutation on n elements. If there is a permutation τ such that applying τ twice has the same effect on the list of elements as applying σ once, we say σ = τ² and τ is a square root of σ.

If we let our n elements be the integers 0 through n − 1, then we can represent permutations by what they do to this list of numbers. In Python as a tuple of length n and compose permutations with the following function:

import itertools

def compose(sigma, tau):
    "Return the composition σ ∘ τ (apply τ first, then σ)."
    return tuple(sigma[j] for j in tau)

We can always construct permutations that have square roots by squaring a permutation. If we run the following code

tau = (3, 1, 4, 5, 2, 0)
sigma = compose(tau, tau)

we find σ = (5, 1, 2, 0, 4, 3), and by construction (3, 1, 4, 5, 2, 0) is a square root of &sigma, though it’s not the only one.

The following code shows that σ has four roots.

import itertools

def numroots(sigma):
    n = len(sigma)
    c = 0
    for tau in itertools.permutations(range(n)):
        if sigma == compose(tau, tau):
            c += 1
    return c

print( numroots(sigma) )
print( numroots( (1, 2, 3, 4, 5, 0) ) )

It also shows that the rotation (1, 2, 3, 4, 5. 0) has no roots.

The function numroots has runtime proportional to n! and so it’s not practical for large permutations. There is a theorem that says a permutation σ has a square root if and only if the number of cycles it has of every even length is even. See [1].

We can also define cubes and cube roots of permutations, and higher powers and roots.

How common is it for permutations to have square roots, or cube roots, etc.? If you pick a random permutation on n elements, what is the probability that it has a kth root?

This is a hard question in general, but it is equivalent to finding the coefficient of x^k in the infinite product

This is theorem 4.8.3 in [1]. This theorem was the motivation for writing about exp_q in the previous post.

Although the product is infinite, there’s no need to compute terms in the product that only contribute powers of x higher than you’re interested in. The following Mathematica code will compute the probability that a permutation on n elements has a kth root.

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}]

So, for example, the probability that a permutation of 10 elements has a square root is 29/96.

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