Number theory

Fermat’s last theorem is so named because it was the last of his asserted theorems to be proved or disproved. But there are variations on another conjectures of Fermat that remain unresolved.

Fermat conjectured that numbers

are always prime. We now call these “Fermat numbers.” Fermat knew that the first five, F₀ through F₄, were all prime.

In some ways, this conjecture failed spectacularly. Euler showed in 1732 that the next number in the sequence, F₅, is not prime by factoring it into 641 × 6700417. So are all the Fermat numbers prime? No.

But that’s not the end of the story. Now we go back and refine Fermat’s conjecture. Instead of asking whether all F_n are prime, we could ask which F_n are prime.

The next several values, F₅ through F₃₂, are all known to be composite. So we might turn Fermat’s original conjecture around: are all F_n composite (for n > 4)? Nobody knows.

Well, let’s try weakening the conjecture. Is F_n composite for infinitely many values of n? Nobody knows. Is F_n prime for infinitely many values of n? Nobody knows that either, though at least one of these two statements must be true!

Here’s why I find all this interesting.

1. It shows how proof by example fails. Fermat knew that the first five numbers in his series were prime. It was reasonable to guess from this that the rest might be prime, though that turned out not to be the case.
2. It shows how what appears to be a dead end can be opened back up with a small refinement of the original question.
3. Like many questions in number theory, the revised question is easy to state but has defied proof for centuries.

In 1989, Star Trek: The Next Generation aired The Royale. In this episode, we learn that Captain Picard tries his hand at proving Fermat’s Last Theorem (FLT) in his spare time. The writers must have believed that FLT would still be unresolved in the 24th century. Four years after The Royale, Andrew Wiles announced his proof of FLT. There was a flaw in Wiles’ first proof, but this was patched two years later in 1995.

Richard Feynman also tried his hand at FLT. He wrote a paper (unpublished) in which he gave a pseudo-proof of FLT based on probability. Feynman said that the probability of FLT being false was “certainly less than 10^-200.” The argument was high creative, sketchy, but ultimately nonsensical. Paul Nahin concludes

Feynman’s probabilistic analysis of Fermat’s Last Theorem would have no mathematical interest at all but for the fact it was Feynman who cooked it up.

Source: Number-Crunching by Paul Nahin.

Here’s an elegant proof from Paul Erdős that there are infinitely many primes. It also does more, giving a lower bound on π(N), the number of primes less than N. [click to continue...]

Suppose you take factorials of a lot of numbers and look at the leading digit of each result. You could argue that there’s no apparent reason that any digit would be more common than any other, so you’d expect each of the digits 1 through 9 would come up 1/9 of the time. Sounds plausible, but it’s wrong.

The leading digits of factorials follow Benford’s law as described in the previous post. In fact, factorials follow Benford’s law even better than physical constants do. Here’s a graph of the leading digits of the factorials of 1 through 500.

In the remainder of this post, I’ll explain why Benford’s law should apply to factorials, make an aside on statistics, and point out an interesting feature of the Python code used to generate the chart above.

Why Benford’s law applies

Here’s a hand-waving explanation. One way to justify Benford’s law is to say that physical constants are uniformly distributed, but on a logarithmic scale. The same is true for factorials, and it’s easier to see why.

The leading digits of the logarithms depend on on their logarithms in base 10. The gamma function extends the factorial function and it is log-convex. The logarithm of the gamma function is fairly flat (see plot here), and so the leading digits of the log-gamma function applied to integers are uniformly distributed on a logarithmic scale.  (I’ve mixed logs base 10 and natural logs here, but that doesn’t matter. All logarithms are the same up to a multiplicative constant. So if a plot is nearly linear on a log10 scale, it’s nearly linear on a natural log scale.)

Update: Graham gives a link in the comments below to a paper proving that factorials satisfy Benford’s law exactly in the limit.

Uniform on what scale?

This example brings up an important principle in statistics. Some say that if you don’t have a reason to assume anything else, use a uniform distribution. For example, some say that a uniform prior is the ideal uninformative prior for Bayesian statistics. But you have to ask “Uniform on what scale?” It turns out that the leading digits of physical constants and factorials are indeed uniformly distributed, but on a logarithmic scale.

Python integers and floating point

I used nearly the same code to produce the chart above as I used in its counterpart in the previous post. However, one thing had to change: I couldn’t compute the leading digits of the factorials the same way. Python has extended precision integers, so I can compute 500! factorial without overflowing. Using floating point numbers, I could only go up to 170!. But when I used my previous code to find the leading digit, it first tried to apply log10 to an integer larger than the largest representable floating point number and failed. Converting numbers such as 500! to floating point numbers will overflow. (See Floating point numbers are a leaky abstraction.)

The solution was to find the leading digit using only integer operations.

def leading_digit_int(n):
    while n > 9:
        n = n/10
    return n

This code works fine for numbers like 500! or even larger.

Related posts:

Benford’s law and SciPy
Physical constants and factorials
Anatomy of a floating point number

Suppose you have a right triangle and every side has rational length. What can you say about the area? For example, is it possible that such a triangle could have area 1?

It turns out the answer is no, and here’s a proof. Suppose you had a right triangle with sides of length a, b, and c with c being the length of the hypotenuse. And suppose a, b, and c are all rational numbers.

Start with the equation

(a² – b²)² = (a² + b²)² – 4a²b²

Suppose the area of the triangle is 1. Then ab/2 = 1 and so ab = 2. Use this and the Pythagorean theorem to rewrite the right side of the equation above. Now we have

(a² – b²)² = c⁴ – 16

But this result contradicts a theorem of Fermat: in rational numbers, the difference of two fourth powers cannot be a square.

So a rational right triangle cannot have area 1. Also, it cannot have area r² for any rational number r. (If it did, you could divide each side by r and have a rational triangle with area 1.)

Now things are about to get a lot deeper.

What values can the area of a rational right triangle take on? Euler called such numbers congruent. Determining easily testable criteria for congruent numbers is an open problem. However, if the Birch and Swinnerton-Dyer conjecture is correct, then an algorithm of Jerrold Tunnell resolves the congruent number problem. (Incidentally, the Clay Mathematics Institute has placed a $1,000,000 bounty on the Birch and Swinnerton-Dyer conjecture.)

What if instead of limiting the problem to rational right triangles we considered all triangles with rational sides? See this paper for such a generalization.

A Mersenne prime is a prime number that is one less than a power of 2. These primes are indexed by the corresponding power of two, i.e. M_p = 2^p – 1. It turns out p must be prime before 2^p – 1 can be prime.

Here are five things I find interesting about Mersenne primes.

1. Record size primes

The largest known prime number is a Mersenne prime, M_43,112,609, proved prime in 2008. And ever since M₅₂₁ was proven prime in 1952, the largest known prime has always been a Mersenne prime (with one exception in 1989). See a history of prime number records.

One reason for the prevalence of Mersenne primes in the record books is that there is a special algorithm for testing whether a number of the form 2^p – 1 is prime, the Lucas-Lehmer test.

2. Finiteness

There may only be a finite number of Mersenne primes. Only 47 are known so far. There are conjectures that predict there are an infinite number of Mersenne primes, but these have not been settled.

3. Connection with perfect numbers

Euclid proved that if M is a Mersenne prime, M(M+1)/2 is a perfect number. Two millennia later, Euler proved that if N is an even perfect number, N must be of the form M(M+1)/2 where M is a Mersenne prime. (More details here.)

Since we only know of 47 Mersenne primes at the moment, and we don’t know of any odd perfect numbers, there are only 47 known perfect numbers.

4. Connection with random number generation

The Mersenne twister is a popular, high-quality random number generator. The generator is so named because its period is a Mersenne prime, M_19,937.

5. History

Mersenne primes are named after the French monk Marin Mersenne (1588–1648) who compiled a list of Mersenne primes. Mersenne wasn’t the first to be aware of such primes. As mentioned above, Euclid connected these primes with perfect numbers.

Marin Mersenne is one of my academic ancestors. I studied under Ralph Showalter, who studied under Tsuan Ting, and so forth back to Frans van Shooten Jr., who studied under Marin Mersenne.

What I find fascinating about this is not my particular genealogy, but that adequate records exist to construct such genealogies. The Mathematics Genealogy Project has over 150,000 records, some reaching back to the Middle Ages.

My wife told me about someone on the radio yesterday discussing voluntary water rationing. People in odd-numbered houses are being asked to water their yards only on odd-numbered days. This person said “I suppose they’re talking about the last digit.”

When she told me about this, my first two thoughts were:

1. Yes, that’s what it means to be odd.
2. Nearly every house number in suburban Houston starts with 1, so going by first digit would be a bad idea.

Strictly speaking, it’s a theorem that odd numbers are those that end in odd digits. The definition of an odd number is one that leaves a remainder of 1 when divided by 2. And in base 10, a number is odd if and only if it ends in an odd digit.

But what if you were using a base other than 10? If the base is even, then a number is odd if and only if the last digit is odd, just like base 10. But what if you’re using an odd base, say base 7? Then the theorem doesn’t hold. For example the number 122 in base 7 is odd, and the number 33 is even. And it’s not just the opposite of the rule for base 10 because 43 is also odd in base 7.

In an odd base, a number is odd iff it has an odd number of odd digits.

(In case you haven’t seen “iff” before, it’s an abbreviation for “if and only if.”)

So, for example, in base 7, the number 642341 is even because it contains two odd digits. And the number 744017 in base 9 is odd because it has three odd digits.

Why does this rule work? Suppose, for example, you have a 4-digit number number pqrs in base b where b is odd. Then pqrs represents

p b³ + q b² + r b + s

All the powers of b are odd, so a number like p times a power of b is odd iff p is odd. So every odd digit in the number contributes an odd number to the sum that expands what the number means. Even digits contribute even terms. A sum is odd iff it has an odd number of odd terms, so a number in an odd base is odd iff it has an odd number of odd digits.

Similar posts:

Divisibility by 7
Divisibility rules in hex
Casting out z’s

Michael Lugo pointed out that the telephone number 867-5309 is prime and may be the largest prime number to appear in the title of a popular song. (The song 867-5309/Jenny peaked at #4 on Billboard in 1982.)

David Radcliffe added “The phone number of Jenny’s twin sister is 8675311″ because 8675309 and 8675311 are twin primes.

How likely is it for a telephone number to be prime? The Prime Number Theorem says that for large values of x, the probability that a number less than x is prime is approximately 1/log(x). Since 1/log(10⁷) = 0.062, about 6% of phone numbers are prime.

We could try to be more accurate. We could look at the probability that a seven-digit number is prime rather than simply a number less than 10⁷ (i.e. excluding numbers with less than seven digits). Or we could use the exact number of primes in a certain range (say using Mathematica’s PrimePi function) rather than using the Prime Number Theorem approximation. But these refinements would still give us estimates of about 6%. Note that not all seven-digit numbers have been assigned as phone numbers, so an exact calculation still gives only an approximate answer.

What about phone numbers with area codes? The Prime Number Theorem suggests about 4.3% of 10-digit numbers are prime. But the US has on the order of 300 area codes, so most 10-digit numbers are not telephone numbers. Also, area codes were not selected randomly. They were selected with a preference for smaller numbers, which means our estimate of 4.3% may be a little low. (We’d expect more prime numbers to start with small area codes.) But I imagine the area code bias has little effect.

Related posts:

Words that are primes base 36
Sonnet primes
Limerick primes

“Casting out nines” is a trick for determining the remainder when a number is divided by nine. Just add the digits of the number together. For example, what’s the remainder when 3896 is divided by 9? The same as when 3+8+9+6 = 26 is divided by 9. We can apply the same trick again: 2+6 = 8, so the remainder is 8.

Casting out nines works because 9 is one less than 10, i.e. one less than the base we use to represent numbers. The analogous trick would work casting out (b-1)’s in base b. So you could cast out 7’s in base 8, or F’s in base 16, or Z’s in base 36.

Why can you cast out (b-1)’s in base b? First, a number written is base b is a polynomial in b. If the representation of a number x is a_na_n-1 … a₁a₀ then

x = a_nbⁿ + a_n-1b^n-1 + … + a₁b + a₀.

Since

b^m – 1 = (b – 1)(b^m-1 + b^m-2 + … + 1)

it follows that b^m leaves a remainder of 1 when divided by b – 1. So a_mb^m leaves the same remainder as a_m when divided by b – 1. If follows that

a_nbⁿ + a_n-1b^n-1 + … + a₁b + a₀

has the same remainder when divided by b – 1 as

a_n + a_n-1 + … + a₁ + a₀

does when it is divided by b – 1.

Related posts:

Words that are prime base 36
Divisibility rules in hex
Divisibility by 7

This morning on Twitter, Alexander Bogomolny posted a link to his article that gives examples of words that are prime numbers when interpreted as numbers in base 36. Some examples are “Brooklyn”, “paleontologist”, and “deodorant.” (Numbers in base 36 are written using 0, 1, 2, …, 9, A, B, C, …, Z as “digits.” )

Tim Hopper replied with a snippet of Mathematica code that lists all words with up to four letters that correspond to base 36 primes.

Rest[ Flatten[ Union[
    DictionaryLookup /@ IntegerString[
        Table[Prime[n], {n, 1, 300000}], 36]]]]

That made me wonder whether you could estimate how many such words there are without doing an exhaustive search.

The Prime Number Theorem says that the probability of a number less than N being prime is approximately 1/log(N). If we knew how many English words there were of a certain length, then we could guess that 1/log(N) of that those words would be prime when interpreted as base 36 numbers. This assumes that forming an English word and being prime have independent probabilities, which may be approximately true.

How well would our guess have worked on Tim’s example? He prints out all the words corresponding to the first 300,000 primes. The last of these primes is 4,256,233. The exact probability that a number less than that upper limit is prime is then

300,000 / 4,256,233 ≈ 0.07.

There are about 4200 English words with four or fewer letters. (I found this out by running

grep -ciE '^[a-z]{1,4}$'

on the words file on a Linux box. See similar tricks here.) If we estimate that 7% of these are prime, we’d expect 294 words from Tim’s program. His program produces 275 words, so our prediction is pretty good.

If we didn’t know the exact probability of a number in our range being prime, we could have estimated the probability at

1/log(4,256,233) ≈ 0.0655

using the Prime Number Theorem. Using this approximation we’d estimate 4200*0.0655 = 275.1 words; our estimate would be exactly correct! There’s good reason to believe our estimate would be reasonably close, but we got lucky to get this close.

Related posts:

Limerick primes
Sonnet primes

Here’s an interesting account of the largest known primes over time. Thanks to @mathematicsprof for pointing this out.

Ever since 1952, the largest known prime has been a Mersenne prime, with one exception in 1989. One reason is that it is simple to test whether Mersenne numbers are prime using the Lucas-Lehmer test. The algorithm is described in seven lines of pseudo-code here.

Here are a couple connections with Mersenne and his primes I’ve written about before. First, Mersenne is one of my mathematical ancestors. Second, Mersenne primes are intimately connected with even perfect numbers, a connection that has been known since Euclid.

Related posts:

Algorithm used for world record pi calculations
Probability that a number is prime

The previous post showed how to list all limerick primes. This post shows how to list all sonnet primes. 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.

Here’s the Mathematica code.

number[s_] := 10100000000000 s[[1]] + 1010000000000 s[[2]] +
    1010000000 s[[3]] + 101000000 s[[4]] +
    101000 s[[5]] + 10100 s[[6]] + 11 s[[7]]

test[n_] := n > 10^13 && PrimeQ[n]

candidates = Permutations[Table[i, {i, 0, 9}], {7}];

list = Select[Map[number, candidates], test];

The function Permutations[list, {n}] creates a list of all permutations of list of length n. In this case we create all permutations the digits 0 through 9 that have length 7. These are the digits a through g.

The function number turns the permutation into a number. This function is applied to each candidate. We then select all 14-digit prime numbers from the list of candidates using test.

If we ask for Length[list] we see there are 16,942 sonnet primes.

As mentioned before, the smallest of these primes is 10102323454577.
The largest is 98987676505033.

Related post:

Limerick primes

The other day, Futility Closet posted this observation:

10102323454577 is the smallest 14-digit prime number that follows the rhyme scheme of a Shakespearean sonnet (ababcdcdefefgg).

I posted this on AlgebraFact and got a lot of responses. One was from Matt Parker who replied that 11551 was the smallest prime with a limerick rhyme scheme.

So how many limerick primes are there? Well, there aren’t many candidates. A limerick prime has to have the form AABBA where A is an odd digit and B is any digit other than A. So for each of five choices of A, there are nine possible B’s. Here’s a Mathematica program to do a brute force search for limerick primes.

For[ j = 0, j < 5, j++,
    For[ k = 0, k < 10, k++,
        x = (2 j + 1)*11001 + 110 k;
        If[ PrimeQ[x], Print[x] ]]]

It turns out there are eight limerick primes:

• 11551
• 33113
• 33223
• 33773
• 77447
• 77557
• 99119
• 99559

See the next post for Mathematica code to list all sonnet primes.

Related posts:

Divisibility by 7
Odd perfect numbers
Twin prime conjecture and the Pentium division bug

Twin primes are pairs of primes that differ by 2. For example, 3 and 5 are twin primes, as are 17 and 19. Importantly, so are 824633702441 and 824633702443. More on that in a minute.

No one knows whether there is a largest pair of twin primes. The twin prime conjecture says that there are infinitely many pairs of twin primes, but the conjecture has not been proven.

Now suppose we take the reciprocals of the twin primes and add them up.

If there were only finitely many twin primes, the sum would have finitely many terms and hence a finite sum. But the sum might converge even though it has infinitely many terms. On the other hand, if we could show that the sum diverges, we’d have a proof of the twin prime conjecture. Viggo Brun showed that the sum does converge. Its sum, known as Brun’s constant, is a little more than 1.9.

In 1994, Thomas Nicely was studying Brun’s constant when he found that his computer incorrectly computed 1/824633702441 beyond the eighth significant figure. Nicely had discovered the infamous Pentium division bug.

Intel responded by saying the division errors were inconsequential. Intel was absolutely correct, but the public couldn’t understand that. They only knew that the chips were “wrong.”

The error was estimated to occur once in every 9 billion divisions. (I doubt any large program has ever been written that is as bug-free as the buggy Pentium chips.) And when an error did occur, the result was not entirely wrong, only less accurate than usual. The public only understood that sometimes the answers were “wrong.” Most people do not understand that floating point arithmetic is nearly always “wrong” in the sense of being less than perfectly accurate.

At first Intel said it would only replace the chips for people who could show they were effected by the bug, i.e. almost nobody. Eventually Intel gave in to pressure and replaced the chips. The episode cost Intel half a billion dollars.

Related posts:

Odd perfect numbers
Probability that a number is prime