Last week I wrote about primailty certificates. These certificates offer a way to verify that a number is prime using less computation than was used to discover than the number was prime.

This post gives a couple more examples of primality certificates using sonnet primes. As described here,

These are primes of the form ababcdcdefefgg, the rhyme scheme of an English (Shakespearean) sonnet, where the letters a through g represent digits and a is not zero.

There are 16,942 sonnet primes, ranging from 10102323454577 to 98987676505033.

Here are Pratt certificates verifying that the largest and smallest sonnet primes are indeed primes.

As before, I assume primes less than 100 are recognizable as primes, though the certificate algorithm continue until all branches of the tree end in a 2.

The quadratic reciprocity theorem addresses the question of whether a number is a square modulo a prime. For an odd prime p, the Legendre symbol

is defined to be 0 if a is a multiple of p, 1 if a is a (non-zero) square mod p, and -1 otherwise. It looks like a fraction, but it’s not. It acts something like a fraction, which was presumably the motivation for the notation.

The quadratic reciprocity theorem says that if p and q are odd primes, then

We can rewrite this in the less symmetric but more useful form

The idea behind quadratic reciprocity as an algorithm is that when p < q, we can use the equation above to turn a question about the existence of square roots mod q into a question about the existence of square roots mod p. This is progress because p is smaller, and we get to reduce q modulo p.

The quadratic reciprocity theorem only applies to the case of odd primes. So we need a couple more facts to complete our algorithm: how to deal with numbers that are not odd or not prime.

If the “numerator” of our Legendre symbol is not prime, we factor it using the theorem that

for numbers m and n that are not multiples of p. And if we get a factor of 2 in a numerator there is a theorem that says

result

Is today a square for Jenny?

As an example, we will determine whether 112023, today’s date, is a square mod 8675309, the telephone number in the 1981 pop song 867-5309/Jenny.

We begin by factoring 112023. A consequence of the factoring rule above is that we can throw away prime factors with an even exponent, and we can replace odd exponents with 1. Therefore

The full calculation is worked out in detail here. We apply the factoring rule and quadratic reciprocity over and over until we get Legendre symbols with such small numbers that we can evaluate them directly. In our case, all we need to know is that 1 is a square mod 3 and 2 is not.

The bottom line of our calculation is 1, so yes, 112023 is a square mod 8675309.

Calculating modular square roots

The quadratic reciprocity theorem gives us a way to determine whether a number has a square root mod p but it does not give us a way to calculate that root. Or rather those roots. Always two there are, like siths.

We could calculate the square roots of 112023 mod 8675309 using the Tonelli-Shanks algorithm, which is unfortunately more complicated than the quadratic reciprocity algorithm.

Mathematica can give us the roots, presumably by using some variation of the Tonelli-Shanks algorithm, or perhaps a newer algorithm. In any case,

    Reduce[x^2 == 112023, x, Modulus->8675309]

tells us the roots are 1955529 and 6719780. The latter is the negative of the former, modulo 8675309, i.e. -1955529 = 6719780 mod 8675309.

You can verify that the roots are correct by computing

1955529*1955529 – 440802*8675309 = 112023

and

6719780*6719780 – 5205053*8675309 = 112023.

Almost all binomial coefficients are divisible by their row number.

This is a theorem from [1]. What does it mean?

If you iterate through Pascal’s triangle, left-to-right and top-to-bottom, noting which entries C(m, k) are divisible by m, the proportion approaches 1 in the limit.

The author proves that the ratio converges to 1, but doesn’t say anything about the rate at which it approaches 1. By writing a little Python script we can observe how the ratio approaches 1.

First, let’s create a generator that will let us walk through Pascal’s triangle.

    def pascal():
        m, k = 0, 0
        while True:
            yield (m, k)
            k += 1
            if k > m:
                m += 1
                k = 0

We can use this generator to illustrate Pascal’s triangle mod row number.

    from math import comb

    p = pascal()
    m, k = next(p)
    while m < 10:
        ch = "*"
        if m > 0 and comb(m, k) % m == 0:
            ch = "0"
        print(ch, end="")
        if m == k:
            print()
        m, k = next(p)

This produces the following.

    *
    00
    *0*
    *00*
    *0*0*
    *0000*
    *0***0*
    *000000*
    *0*0*0*0*
    *00*00*00*

The theorem says that as we keep going, the proportion of 0’s in a diagram like the one above approaches 1.

Now let’s plot the proportion of zeros as we traverse Pascal’s triangle mod row number.

    N = 1000
    x = np.zeros(N)
    p = pascal()
    for n in range(N):
        m, k = next(p)
        if m > 0 and comb(m, k) % m == 0:
            x[n] = 1
    y = np.arange(N)
    plt.plot(y, x.cumsum()/y)
    plt.show()

Here’s what we get.

The ratio doesn’t seem to be converging to 1. If I had to guess, I might say it’s converging to 0.7 by looking at this plot. But if we go further out and switch to a log scale it seems more plausible that the sequence is converging to 1, albeit very slowly.

[1] Heiko Harborth. Divisibility of Binomial Coefficients by Their Row Number. The American Mathematical Monthly, Jan. 1977, Vol. 84, No. 1, pp. 35–37.

Yesterday I wrote about the Chicken McNugget problem: if Chicken McNuggets are sold in boxes of 6, 9, and 20, what’s the largest number of nuggets you cannot buy? That post showed that the solution was 43. The technical name for these kinds of problems is numerical monoids.

The method of solution in the previous post is reusable. But before we reuse it, let’s look at a simpler problem involving two units rather than three.

Cicada cycles

Some cicadas appear every 13 years and some every 17 years. So, for example, 43 years from now would be two 13-year cicada cycles and one 17-year cicada cycle. We could ask what the largest number of years is that cannot be expressed in terms of cicada cycles. This is known as the Frobenius number of 13 and 17.

For one thing, there is a largest number of years that cannot be expressed in terms of cicada cycles. If x and y are positive integers, there is an upper bound on the numbers that cannot be expressed as sum of non-negative integer multiples of x and y if and only if x and y are relatively prime. And if x and y are relatively prime, this upper bound is

xy − x − y.

J. J. Sylvester proved this back when Chester A. Arthur was US president [1]. It’s not quite obvious that x and y being relatively prime is sufficient for there to be an upper bound—it follows from Bézout’s lemma—but it is obvious it’s necessary: if x and y are multiples of k > 1, then all sums of integer multiples of x and y are also multiples of k, and so non-multiples of k are unreachable.

So Sylvester’s theorem shows that the Frobenius number of 13 and 17 is 13×17 – 13 – 17 = 191.

Throwing in Venus

Venus goes around the sun 13 times in the time it takes Earth to go around 8 times. That is, Venus and Earth line up every 8 years (from the perspective of our planet).  What if we throw in multiples of 8 year cycles? We could, for example, express 25 years as one 17-year cicada cycle and one cycle of Venus.

We know there’s an upper bound on numbers that cannot be expressed as non-negative integer combinations of 13, 17, and 8 because we can express every number greater than 191 using 13 and 17 alone. Adding a new generator cannot increase the Frobenius number. In fact, in this case we know it lowers it because we could apply Sylvester’s theorem to 8 and 13 and show that the Frobenius number of 8, 13, and 17 is at most 83.

Algorithm

There is no analog of Sylvester’s theorem for more than two generators, but we could always compute the Frobenius number in a particular case using the approach from the McNugget post: use brute force up to a point, then use theory to argue that’s enough. In that post, the smallest generator was 6, and we used brute force until we found 6 consecutive numbers that could be produced as non-negative integer combinations of our generators. We could improve on this approach by using Sylvester’s theorem: we use brute force for all combinations that could possibly generate a result lower than an upper bound obtained by Sylvester’s theorem. There are more efficient algorithms—see [2] for references—but the following will work.

    from math import gcd, floor
    from itertools import product

    def frobenius(generators):
        "assumes the smallest two generators are relatively prime"
        g = sorted(generators)
        assert(len(g) > 1)
        assert(gcd(g[0], g[1]) == 1)
        ub = g[0]*g[1] 
        
        limits = [range(int(floor(ub/n))) for n in generators]
        reachable = set()
        for p in product(*limits):
            combo = sum(p[i]*g[i] for i in range(len(g)))
            if combo < ub:
                reachable.add(combo)
        complement = set(range(ub)) - reachable
        return(max(complement))
    
    print( frobenius([13, 17, 8]) )

This shows that the Frobenius number of 8, 13, and 17 is 44.

Just to do another example, the code can show that the smallest number of days not representable by multiples of 7, 30, and 365 is 209 days.

Notes

[1] J. J. Sylvester, Mathematical questions with their solutions, Educational Times 41 (1884) 171–178.

[2] Factoring in the Chicken McNugget monoid. arXiv:1709.01606

Let s be an integer equal to at least 5 and let n be a positive integer.

Look at all tuples of s integers, each integer being between 0 and n-1 inclusive, and look at the sum of their squares mod n. About 1/n will fall into each possible value.

So, for example, if you look at a lists of 6 digits, sum the squares of the digits in each list, then about 1/10 of the results will end in 0, about 1/10 will end in 1, about 1/10 will end in 2, etc.

This is the gist of a theorem discussed on Math Overflow here.

And here is Python code to demonstrate it.

    from itertools import product

    n = 10
    s = 6

    counts = [0]*n
    for p in product(range(n), repeat=s):
        counts[ sum(x**2 for x in p) % n ] += 1

    print(counts)

The result is

    [103200, 99200, 99200, 99200, 99200, 
     103200, 99200, 99200, 99200, 99200]

and so between 99,200 and 103,200 sums end in each digit.

Let’s run it again with n = 11 and s = 5. Now we get

    [14641, 14762, 14520, 14762, 14762, 14762, 
     14520, 14520, 14520, 14762, 14520]

which shows the counts are between 14,520 and 14,762. If they were perfectly uniformly distributed we’d expect 11⁴ = 14,641 in each bin. The maximum deviation from a uniform distribution is 0.83%.

The code becomes impractical as s or n get larger, but we can try s = 7 and keep n = 11. Then we get a maximum deviation from uniformity of 0.075%, i.e. about 10 times closer to being uniform.