Pythagorean triangles with side 2023

Can a Pythagorean triangle have one size of length 2023? Yes, one possibility is a triangle with sides (2023, 6936, 7225).

right triangle with sides 2023, 6936, and 7225

Where did that come from? And can we be more systematic, listing all Pythagorean triangles with a side of length 2023?

Euclid’s formula generates Pythagorean triples by sticking integers m and n into the formula

(m² – n², 2mn, m² + n²).

If one of these three numbers equals 2023, it’s going to have to be the first one. Why? The second number is even, and 2023 is odd, so that one is ruled out.

It’s less obvious but true that the hypotenuse m² + n² cannot be 2023. We could verify this by brute force, trying integer values of m between 1 and √2023, subtracting m² from 2023, and seeing whether the result is a square. But there is a more elegant approach.

There is a theorem that says an integer N is the sum of two squares exactly when every prime that divides N is either congruent to 1 mod 4 or appears an even number of times in the factorization of N. The prime factorization of N is 7 × 17 × 17. The factor 7 appears to an odd power and 7 = 3 mod 4.

Is 2023 a difference of squares? Yes, every odd number is a difference of squares.

To write an number N as a difference of squares, we have to solve

m² – n² = (m + n)(m – n) = N.

If N is odd, we can always set mn = 1. In the case N = 2023, this means m = 1012 and n = 1011.

The solution above corresponds to factoring 2023 as 2023 × 1. We could also factor 2023 as 289 × 7 or 119 × 17. (The first factor has to be larger than the second for mn to be positive.) To find more possible value of m and n we solve the equations

m + n = 289
mn = 7

and

m + n = 119
mn = 17.

These have solution (148, 141) and (68, 51).

To summarize, we can find three Pythagorean triples with one component equal to 2023 using Euclid’s formula. These are

(2023, 2046264, 2046265)
(2023, 41736, 41785)
(2023, 6936, 7225)

corresponding to (m, n) = (1013, 1012), (148, 141), and (68, 51) respectively. The first two are primitive Pythagorean triples because the three numbers are relative prime in both cases. The last triple is not primitive because each number is divisible by 289.

Euclid’s formula generates all primitive Pythagorean triples, but it doesn’t generate all possible Pythagorean triples. To find all Pythagorean triples we have to consider solutions of the form

(km² – kn², 2kmn, km² + kn²)

for some integer k. If 2023 is k times a difference of squares, k must be a factor of 2023, and so k = 1, 7, 17, 119, 289, or 2023. We’ve already considered the case k = 1.

For k = 7, we express 289 as a difference of squares. The only solutions with positive integers is (m, n) = (145, 144). This yields the Pythagorean triple

(2023, 292320, 292327).

For k = 17, we express 119 as a difference of squares. The solutions are (m, n) = (60, 59) and (12, 5). These yield

(2023, 120360, 120377)

and

(2023, 2040, 2873).

For k = 119, we express 17 as a difference of squares. This can only be done one way, (m, n) = (9, 8). This yields

(2023, 17136, 17255).

For k = 289, we express 7 as a difference of squares. The unique solution (4, 3) yields the triple

(2023, 6936, 7225)

that we had found previously.

For k = 2023, we have to write 1 as a difference of squares. This has no solution if we require both squares to be positive.

So we’re done. We’ve found all Pythagorean triples containing 2023.

Related posts

Factoring b^n + 1

The previous post illustrated a technique for finding factors of number of the form bn – 1. This post will look at an analogous, though slightly less general, technique for numbers of the form bn + 1.

There is a theorem that says that if m divides n then bm + 1 divides bn + 1 if n is odd. To see that the exponent must be odd, let b = 10, n = 4, and m = 2: 10,0001 is not divisible by 101.

Last time we looked at 22023 – 1, so this time let’s look at J = 22023 + 1.

Since 2023 = 7×17×17, we know J is divisible by 2n + 1 for n = 7, 17, 7×17, and 17×17. Factoring each of these factors tells us that J is divisible by

  • 3
  • 43691
  • 72251
  • 79187
  • 1077971
  • 823679683
  • 18360250452977
  • 197766803208315851
  • 143162553165560959297
  • 338858733065598401355195539629373089

(If it seems this post is moving too fast, you might want to go back and read the previous post and then come back.)

Let P be the product of the factors above, and let F = J/P.  Now F is a 493-digit number, so we don’t expect to be able to get much further unless we’re lucky. We can find three factors of F

  • 43
  • 4382432993
  • 1809032819104754041

but then we get stuck. We know there are more factors, but it seems doubtful that I’m going to find any more anytime soon.

***

As in the previous post, we’ll conclude with Python code to verify the claims above.

from sympy import isprime

def verify_factors(N, factors, full=True):
    prod = 1
    for f in factors:
        assert(isprime(f))
        prod *= f
    if full:
        assert(N == prod)
    else:
        assert(N % prod == 0)
        assert(not isprime(N // prod))
    return N // prod

factors = [
    3,
    43691,
    72251,
    79187,
    1077971,
    823679683,
    18360250452977,
    197766803208315851,
    143162553165560959297,
    338858733065598401355195539629373089
]

J = 2**2023 + 1
F = verify_factors(J, factors, False)

factors = [43, 4382432993, 1809032819104754041]
verify_factors(F, factors, False)

print(len(str(F)))

Factoring b^n – 1

Suppose you want to factor a number of the form bn – 1.

There is a theorem that says that if m divides n then bm – 1 divides bn – 1.

Let’s use this theorem to try to factor J = 22023 – 1, a 609-digit number. Factoring such a large number would be more difficult if it didn’t have a special form that we can exploit.

Our theorem tells us J is divisible by 27 – 1 and 2289 – 1 because 2023 = 7×17×17.

Is J divisible by (27 – 1)(2289 – 1)? In fact it is, but this doesn’t necessarily follow from the theorem above because (bm – 1) and (bn/m -1) could share factors, though in this case they do not.

So we have

J = (27 – 1) (2289 – 1) F

for some factor F.  Now 27 – 1 = 127, a prime. What about 2289 – 1?

By the theorem above we know that 2289 – 1 is divisible by 217 – 1 = 131071, which turns out to be a prime. We can get a couple more factors of 2289 – 1 by consulting The Cunningham Project, a project that catalogs known factors of bn ± 1 for small values of b. We learn from their site that 2289 – 1 is also divisible by 12761663 and 179058312604392742511009. All together

2289 – 1 = 131071 × 12761663 × 179058312604392742511009 × R

where

R = 3320934994356628805321733520790947608989420068445023

and R turns out to be prime.

So now we have five prime factors of J:

  • 127
  • 131071
  • 12761663
  • 179058312604392742511009
  • R.

That leaves us with F above, a 520-digit number. It would seem we’ve gotten all the mileage we can out of our factoring trick. But there’s something we haven’t tried: We know that J is divisible by 2119 – 1 because 7 × 17 = 119 is a factor of 2023.

Now

2119 – 1 = 127 × 239 × 20231 × 131071 × 62983048367 × 131105292137

and so these prime factors divide J. However, two of these, 127 and 131071, we’ve already discovered. But we do learn 4 more prime factors. So

F = 239 × 20231 × 62983048367 × 131105292137 × G

where G is a 492-digit number. We can tell by Fermat’s test that G is not prime, but I’m unaware of any clever technique for easily finding any of the factors of G.

***

In general factoring a 492-digit number is hard. There are RSA challenge numbers smaller than this that have not yet been factored, such as RSA-260, a 260-digit number. On the other hand, the RSA numbers are designed to be hard to factor. RSA-260 has two prime factors, presumably both the same order of magnitude. We get a little luckier with G. It has three relatively small factors that I was able to find:

G = 166684901665026193 × 3845059207282831441 × 153641005986537578671 × H

where H is a 436-digit number. I know from Fermat’s test that H is composite but I could not find any factors.

Update: From the comments, 2605053667526976413491923719 is also a factor of G. I’ve updated the code below accordingly. Now the unfactored part H is a 408-digit number.

***

Here’s Python code to verify the claims made above.

from sympy import isprime

def verify_factors(N, factors, full=True):
    prod = 1
    for f in factors:
        assert(isprime(f))
        prod *= f
    assert(N % prod == 0)
    if full:
        assert(N == prod)

R = 3320934994356628805321733520790947608989420068445023
factors = [131071, 12761663, 179058312604392742511009, R]
verify_factors(2**289 - 1, factors)

J = 2**2023 - 1
prod = 127*(2**289 - 1)
F = J // prod
assert(J == prod*F)

factors = [127, 239, 20231, 131071, 62983048367, 131105292137]
verify_factors(2**119 - 1, factors)

prod = 239*20231*62983048367*131105292137
G = F // prod
assert(F == G*prod)

factors = [166684901665026193, 3845059207282831441, 153641005986537578671, 2605053667526976413491923719]
verify_factors(G, factors, False)

prod = 1
for f in factors:
    prod *= f
H = G // prod
assert(G == H*prod)
assert(not isprime(H))

assert(len(str(J)) == 609)
assert(len(str(F)) == 520)
assert(len(str(G)) == 492)
assert(len(str(H)) == 408)

Related post: Primality certificates

Update: See the next post for the case of bn + 1.

Special primality proofs

I’ve written lately about two general ways to prove that a number is prime: Pratt certificates for moderately-large primes and elliptic curve certificates for very large primes.

If you can say more about the prime you wish to certify, there may be special forms of certificates that are more efficient. In particular, there are efficient tests to determine whether Fermat number or a Mersenne number is prime.

Pepin’s test, implemented in four lines of Python here, determines whether or not a Fermat number is prime. It can instantly prove that the known Fermat primes are indeed prime. It can quickly show that the next several Fermat numbers are composite, but the time required to run Pepin’s test increases rapidly for larger numbers.

The Lucas-Lehmer test, implemented in six lines of Python here, tests whether a Mersenne number is prime. The largest known primes are Mersenne primes because the Lucas-Lehmer test is relatively efficient, though it still takes a lot of effort to search for new Mersenne primes.

Certifying that a number is composite is potentially much easier than certifying that a number is prime. A famous example of a proof that a number is composite was Euler’s announcement in 1732 that

232 + 1 = 641 × 6700417,

thereby proving that the fifth Fermat number is not prime, disproving Fermat’s conjecture.

Several Fermat numbers have been fully factored, concretely proving that they are composite, but some Fermat numbers have been proven composite via Pepin’s test that have not been factored. The analogous statement is true for Mersenne primes and the Lucas-Lehmer test.

Elliptic curve primality certificates

I’ve written recently about a simple kind of primality certificates, Pratt certificates. These certificates are easy to understand, and easy to verify, but they’re expensive to produce. In order to produce a Pratt certificate that n is a prime you have to factor n-1, and that can take a long time if n is large [1]. So Pratt certificates are only used in practice for moderately large primes. “Moderate” could mean a 50-digit number, which is large for some purposes, but would not be considered large in the context of cryptography.

High-level description

Elliptic curve primality proofs are used to certify primes that are too large for Pratt certificates. The Goldwasser-Kilian theorem says that a number p is prime if there exists an elliptic curve with certain properties that depend on another prime q. As with Pratt certificates, this process is applied recursively: p is prime if q is prime, and q is prime if r is prime, … down some the base case of primes so small that they are recognized as prime.

The Goldwasser-Kilian theorem gives a lower bound on the size of q:

q > p1/2 + p1/4 + 1.

Ideally you find a q that’s not too much bigger than it needs to be. So to prove a 100-digit number is prime, maybe you reduce the problem to proving than an 80-digit number is prime, then reduce that to the problem of proving a 64-digit number is prime, etc.

Details

What is this elliptic curve in the theorem? You have to find a smooth curve mod p and a point M on the curve such that M is not the group identity, but qM is the group identity. How in the world would you go about finding such a curve? That’s a more complicated story, but the main idea of certificates is that they’re easy to verify, not easy to produce. Coming up with an elliptic curve primality proof is harder than verifying one [2].

Primality certificates need to be computationally easy to verify, not necessarily conceptually easy to verify. But they tend to be conceptually easier to verify than to produce as well. The rest of this post will outline how you could write your own code to verify an elliptic curve primality certificate.

How to verify from scratch

To verify an elliptic curve primality proof, you first [3] have to verify that someone has given you a smooth elliptic curve. They give you two numbers, a and b, and claim that

y² = x³ + ax + b

is a smooth elliptic curve over the integers mod p. To verify this you have to compute

4a³ + 27b²

and show that it is not divisible by n. Then you have to verify that the point M = (x, y) they gave you is on the curve, which requires sticking x and y into the cubic equation above and verifying that the equation holds mod p.

The last and most difficult step is to verify that qM gives the group identity for the elliptic curve. You have to add M to itself q times using the group law of the curve [4].

Implementing the group law for an elliptic curve modulo a prime is a little tedious, but the only difficult part is finding a multiplicative inverse mod p. Presumably your computing environment would have a library routine for modular inverses. For example, Mathematica has the function ModularInverse.

Related posts

[1] If n is a safe prime then q = (p-1)/2 is also a prime, and so proving that q is prime is just about as hard as proving p is prime. Such hard luck can’t keep going. As you recurse you eventually get to primes p such that p-1 factors into pieces much smaller than p, but the process can start of being difficult.

[2] Elliptic curve primality certificates are sometimes called Atkin-Goldwasser-Kilian-Morain certificates or Atkin-Morain certificates because Atkin and Morain did a lot of work to make the idea of Goldwasser and Kilian practical.

[3] The Goldwasser-Kilian theorem also requires that p is not divisible by 2 or 3, but that’s trivial to verify.

[4] In practice you don’t literally add M to itself q times, but you produce a result that is the same as if you did. You’d use repeated doubling, analogous to using repeated squaring in fast exponentiation.

Primes with two non-zero bits

Suppose a number n written in binary has two 1s and all the rest of its bits are zeros. If n is prime, then the 1s must be the first and last bits of n. The first bit is 1 because the first bit of every positive integer is 1. The last bit is 1 because if the second 1 appeared anywhere else then n would be divisible by 2.

Here are the first four examples:

11two = 3ten
101two = 5ten
10001two = 17ten
100000001two = 257ten

Let’s count how many 0 bits each number contains. We get 0, 1, 3, and 7 if we count in base 10, or 0, 1, 11, and 111 if we count in binary. It seems that the number of 0 bits is either zero or a number whose binary representation is all 1s. In fact this turns out to be a theorem:

Suppose a number n is written in binary as a 1, followed by m 0’s, and finally a 1. Then if n is prime, then either m = 0 or the binary representation of m consists entirely of 1’s.

To prove this theorem, let’s move away from talking about bits and move toward more standard math terminology.

Our description of n in terms of bits is another way of saying

n = 2k + 1.

And our description of m is another way of saying

m = 2j – 1.

And a 1 followed by m 0’s and then a 1 is a number of the form

2m+1 + 1.

So the theorem above is saying that if 2k + 1 is a prime number, then k is a power of 2.

Here’s a proof. Suppose k is not a power of 2. Then k has an odd factor, i.e. k = ab where b > 1 is an odd number. In that case we can factor 2k + 1 as follows:

2^k + 1 = (2^a)^b + 1 = (2^a + 1)(2^{a(b-1)} - 2^{a(b-2)} + 2^{a(b-3)} \cdots + 1)

We’ve shown above that 2k + 1 is prime when k = 2j and j = 0, 1, 2, and 3.

What about j = 4? Yes, 216 + 1 = 65537 is prime.

We’re on a roll. What about j = 5? Is 232 + 1 prime? No, it’s not.

What about larger values of j? Do any of them produce primes? Nobody knows.

This post has secretly been about Fermat primes. The nth Fermat number is given by

F_n = 2^{2^n} + 1

Fermat observed that Fn was prime for n = 0, 1, 2, 3, and 4. He conjectured that Fn is prime for all n, but Euler factored F5, disproving Fermat’s conjecture.

Why did this post delay using standard terminology? The first reason was to unpack the math. The direct discussion can go by so quickly that it’s hard to appreciate what it is going on. The second reason was to illustrate how much easier things can become when we switch from discussing bit patterns to using math notation. Sometimes things go the other way, with bits being easier to understand than mathematical equations.

Certified sonnet primes

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.

Pratt Primality Certificates

The previous post implicitly asserted that J = 8675309 is a prime number. Suppose you wanted proof that this number is prime.

You could get some evidence that J is probably prime by demonstrating that

2J-1 = 1 mod J.

You could do this in Python by running the following [1].

    >>> J = 8675309
    >>> assert( pow(2, J-1, J) == 1 )

This shows J is a probable prime to base 2.

If you want more evidence, you could also show J is a probable prime to base 3.

    >>> assert( pow(3, J-1, J) == 1 )

But no matter how many bases you try, you won’t have proof that J is prime, only evidence. There are pseudoprimes, (rare) composite numbers that satisfy the necessary-but-not-quite-sufficient conditions of Fermat’s primality test.

Primality certificates

A primality certificate is a proof that a number is prime. To be practical, a certificate must be persuasive and efficient to compute.

We could show that J is not divisible by any integer less than √J. That would actually be practical because J is not that large.

    >>> for n in range(2, int(J**0.5)+1):
    ...     assert(J % n > 0)

But we’d like to use J to illustrate a method that scales to much larger numbers than J.

Pratt certificates

Pratt primality certificates are based on a theorem by Lucas [2] that says a number n is prime if there exists a number a such that two conditions hold:

an-1 = 1 mod n,

and

a(n-1)/p ≠ 1 mod n

for all primes p that divide n-1.

How do you find a? See this post.

Example

To find a Pratt certificate for J, we have to factor J-1. I assert that

J-1 = 8675308 = 4 × 2168827

and that 2168827 is prime. Here’s verification that Lucas’ theorem holds:

    >>> assert( pow(2, (J-1)//2, J) != 1 )
    >>> assert( pow(2, (J-1)//2168827, J) != 1 )

What’s that you say? You’re not convinced that 2168827 is prime? Well then, we’ll come up with a Pratt certificate for 2168827.

Pratt certificates are generally recursive. To prove that p is prime, we have to factor p-1, then prove all the claimed prime factors of p-1 are prime, etc. The recursion ends when it gets down to some set of accepted primes.

Now I assert that

    2168827 - 1 = 2168826 = 2 × 3 × 11 × 17 × 1933

and that all these numbers are prime. I’ll assume you’re OK with that, except you’re skeptical that 1933 is prime.

The following code is proof that 2168827 is prime, assuming 1933 is prime.

    >>> m = 2168827
    >>> for p in [2, 3, 11, 17, 1933]:
    ...     assert( pow(3, (m-1)//p, m) != 1 )

Finally, we’ll prove that 1933 is prime.

You can verify that

    1933 - 1 = 1932 = 2² × 3 × 7 × 23

and I assume you’re convinced that each of these factors is prime.

    >>> m = 1933
    >>> for p in [2, 3, 7, 23]:
    ...     assert( pow(5, (m-1)//p, m) != 1 )

Pratt certificates can be written in a compact form that verification software can read. Here I’ve made the process more chatty just to illustrate the parts.

Update: Here’s a visualization of the process above, drawing arrows from each prime p to the prime factors of p-1.

In this post we’ve agree to recognize 2, 3, 7, 11, 17, and 23 as primes. But you only have to assume 2 is prime. This would make the software implementation more elegant but would make the example tedious for a human consumption.

Efficiency

A primality certificate does not have to be efficient to produce, though of course that would be nice. It has to be efficient to verify. You could imagine that the prime number salesman has more compute power available than the prime number customer. In the example above, I used a computer to generate the Pratt certificate, but it wouldn’t be unreasonable to verify the certificate by hand.

The brute force certificate above, trying all divisors up to √p, obviously takes √p calculations to verify. A Pratt certificate, however, takes about

4 log2 p

calculations. So verifying a 10-digit prime requires on the order of 100 calculations rather than on the order of 100,000 calculations.

Atkin-Goldwasser-Kilian-Morain certificates

Producing Pratt certificates for very large numbers is difficult. Other certificate methods, like Atkin-Goldwasser-Kilian-Morain certificates, scale up better. Atkin-Goldwasser-Kilian-Morain certificates are more complicated to describe because they involve elliptic curves.

Just as Pratt took a characterization of primes by Lucas and turned it into a practical certification method, Atkin and Morain turned a characterization of primes by Goldwasser and Kilian, one involving elliptic curves, and turned it into an efficient certification method.

These certificates have the same recursive nature as Pratt certificates: proving that a number is prime requires proving that another (smaller) number is prime.

Update: More on elliptic curve primaility proofs.

***

[1] This is a more efficient calculation than it seems. It can be done quickly using fast exponentiation. Note that it is not necessary to compute 2J-1 per se; we can carry out every intermediate calculation mod J.

[2] Lucas was French, and so the final s in his name is silent.

Quadratic reciprocity algorithm

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

\left(\frac{a}{p}\right)

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

\left(\frac{p}{q}\right) \left(\frac{q}{p}\right) = (-1)^{\frac{p-1}{2} \frac{q-1}{2}}

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

\left(\frac{p}{q}\right) = (-1)^{(p-1)(q-1)/4} \left(\frac{q}{p}\right)

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

\left(\frac{m}{p}\right) \left(\frac{n}{p}\right) = \left(\frac{mn}{p}\right)

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

\left(\frac{2}{p}\right) = (-1)^{\frac{p^2 - 1}{8}}

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

\begin{align*} \left(\frac{112023}{8675309}\right) &= \left(\frac{3^5\cdot 461}{8675309}\right) \\ &= \left(\frac{3}{8675309}\right)^5 \left(\frac{461}{8675309}\right)\\ &= \left(\frac{3}{8675309}\right) \left(\frac{461}{8675309}\right) \end{align*}

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.

Pascal’s triangle mod row number

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.