Triangular Numbers Revisited

The $n$th triangular number is

$$T(n) = 1 + 2 + 3 + \cdots + n = \frac{n(n+1)}{2}.$$

Triangular numbers satisfy many interesting properties, including a product rule for $m, n \ge 1$:

$$T(mn) = T(m) T(n) + T(m-1) T(n-1).$$

This product rule can be interpreted geometrically by subdividing a large triangle into smaller triangles. This picture illustrates the case $m = 5$ and $n = 4$.

20th-triangular-number

An interesting question (in my opinion) is whether there are other sequences that satisfy this identity. It turns out that there are exactly five such sequences:

  1. $T(n) = 0$ for all $n \ge 0$.
  2. $T(n) = \frac12$ for all $n \ge 0$.
  3. $T(0) = 0$ and $T(2n) = T(2n-1) = n$ for all $n \ge 1$.
  4. $T(3n) = T(3n+2) = 0$ and $T(3n+1) = 1$ for all $n \ge 0$.
  5. $T(n) = \frac12 n(n+1)$ for all $n \ge 0$.

I wrote about this problem previously at MathBlag, but now I am trying to prepare the result for publication. I would appreciate feedback on my e-print.

Primes with average digit 1

In this note, we will use Python to explore the following problem: How many primes are there whose average digit is 1? In other words, how many \(k\)-digit prime numbers exist whose sum of digits is \(k\)?

Here is a direct approach. We generate all prime numbers up to some limit and calculate the sum of digits for each prime. We use the primerange function from the SymPy package to generate the prime numbers.

from sympy import primerange

def digitsum(n):
    s = 0
    while n > 0:
        s += (n % 10)
        n //= 10
    return s

def primes_avg1(N):
    for n in primerange(2, N):
        if digitsum(n) == len(str(n)):
            yield n

Listed below are the prime numbers less than a million whose average digit is 1. Note that there are no examples with 3 or 6 digits. In fact, there is never an example with \(3k\) digits, because a number is divisible by 3 if and only if the sum of its digits is divisible by 3.

print list(primes_avg1(10**6))

[11, 1021, 1201, 2011, 3001, 
10103, 10211, 10301, 11003, 
12011, 12101, 13001, 20021, 
20201, 21011, 21101, 30011]

This method is inefficient for large numbers because prime numbers are relatively common, even when the number of digits is large. According to the prime number theorem, the number of prime numbers less than \(N\) is approximately \(N/\ln(N)\). The exact number of primes having \(k\) digits or fewer has been computed up to \(k = 26\). (See http://en.wikipedia.org/wiki/Prime-counting_function)

But the number of \(k\)-digit numbers with average digit 1 increases more slowly as a function of \(k\) than the number of \(k\)-digit prime numbers. We can estimate this number as follows. Any positive integer can be represented by a pattern of stars and bars, where the bars separate the stars into groups. For example, the pattern shown below represents the number 4102 because there are 4 stars in the first group, 1 star in the second group, 0 stars in the third group, and 2 stars in the fourth group.

★ ★ ★ ★ | ★ |   | ★ ★

If a number has \(k\) digits that sum to \(k\), then it can be represented by a certain arrangement of \(k\) stars and \((k-1)\) bars. The total number of permutations of \(k\) stars and \((k-1)\) bars is

\[\binom{2k-1}{k-1} = \frac{(2k-1)!}{(k-1)!\,k!}\]

which is less than \(2^{2k-1}\), since the sum of the numbers in row \((2k-1)\) of Pascal’s triangle is \(2^{2k-1}\). Therefore, the number of \(k\) digits with sum of digits \(k\) is less than \(2^{2k-1}\). Although this quantity grows exponentially with \(k\), it grows much more slowly than the number of \(k\)-digit primes.

The exact count of \(k\)-digit numbers with digit sum \(k\) is given by the following summation, which can be obtained by means of the inclusion-exclusion principle.

\[\sum_{j\ge0} (-1)^j \binom{k}{j} \binom{2k-10j-2}{k-1} .\]

Wait, where are you going with this?

Since \(k\)-digit primes are common, while \(k\)-digit numbers with digit sum \(k\) are rare, we realize that filtering a list of primes was a mistake. It is much better to generate a list of numbers having average digit 1, and test each of these numbers for primality. Fortunately, this list can be generated using an efficient recursive algorithm.

def gen(k, s):
    """k-digit numbers
       with digit sum s"""
    if k == 1:
        if s < 10:
            yield s
    else:
        for a in range(min(s, 10)):
            for m in gen(k-1, s-a):
                yield 10*m + a

As a further refinement, we can use the fact that the last digit of a prime greater than 5 is always 1, 3, 7, or 9. Here is our Python code for generating \(k\)-digit prime numbers with sum of digits \(s\).

from sympy.ntheory import primetest

def primes(k, s):
    for a in 1, 3, 7, 9:
        if s <= a:
            break
        for m in gen(k-1, s-a):
            n = 10*m + a
            if primetest.isprime(n):
                yield n

We can use this to function to count the number of \(k\)-digit primes with digit sum \(k\). Since we are using generators instead of lists, our code does not use a lot of memory. However, it is time consuming to test a large number for primality, so a probabilistic test is used.

for k in range(2, 18):
    if k % 3 > 0:
        print k, sum(1 
          for p in primes(k, k)) 

The output is shown below. The code uses the Miller-Rabin primality test, which is guaranteed to be correct for \(n < 10^{16}\) but has a small chance of error \((<10^{-28})\) for large inputs. Consequently, the last value in the table should be viewed as highly probable, but not certain.

Digits Count
2 1
4 4
5 12
7 95
8 212
10 2395
11 10657
13 126068
14 375941
16 4943357
17 20513691

The square root of 11

Here are four proofs of the irrationality of the square root of 11. Can you think of another proof?

Unique factorization

Suppose that the square root of 11 is rational. Then we can express the square root of 11 as a fraction of integers.

\[\sqrt{11} = \frac{a}{b}\]

Therefore, \(a^2 = 11b^2\).

Now, 11 occurs an even number of times in the prime factorization of \(a^2\), but it occurs an odd number of times in the prime factorization of \(11b^2\). This is a contradiction, since prime factorizations are unique.

See James Tanton’s essay for a more detailed explanation.

Euclid’s lemma

Euclid’s lemma states that if a prime number divides the product of two numbers, then it divides one or both of the factors.

Suppose that the square root of 11 is rational. Then we can express the square root of 11 as a fraction of integers \(a/b\). We may assume that \(a\) and \(b\) have no common factor other than 1.

As before, we write \(a^2 = 11b^2\). Since 11 is a prime number and it divides \(a^2\), Euclid’s lemma implies that 11 also divides \(a\).

Let \(a = 11c\). Then \((11c)^2 = 11b^2\), hence \(11c^2 = b^2\). This time, Euclid’s lemma implies that 11 divides \(b\). This contradicts the assumption that \(a\) and \(b\) have no common factor.

Rational root theorem

The rational root theorem states that if \(x\) is a rational root of a polynomial with integer coefficients, then \(x = p/q\), where \(p\) divides the constant term of the polynomial and \(q\) divides the leading coefficient.

The square root of 11 is a root of the polynomial \(x^2 – 11\). According to the rational root theorem, the only possible rational roots are \(\pm 1\) or \(\pm 11\). Since none of these are equal to the square root of 11, we conclude that the square root of 11 is irrational.

Infinite descent

Suppose that the square root of 11 is rational. Then there exists a positive integer \(n\) so that \(m := n\sqrt{11}\) is an integer. Let \(n\) be the smallest positive integer with this property.

Then
\[(m – 3n)\sqrt{11} = (n\sqrt{11}-3n)\sqrt{11} = 11n-3n\sqrt{11} = 11n-3m\]
which is also an integer.

On the other hand, \(0 < m - 3n < n\), since \(3 < \frac{m}{n} < 4\). This is a contradiction, because we assumed that \(n\) is the smallest positive integer so that \(n\sqrt{11}\) is an integer.

BONUS: Continued fraction

The square root of 11 has an infinite continued fraction representation.
\[
\sqrt{11} = [3; \overline{3, 6}]
= 3 + \cfrac{1}{3
+ \cfrac{1}{6
+ \cfrac{1}{3
+ \cfrac{1}{6 + \cdots} } } }
\]

But the continued fraction representation of a rational number is finite. Therefore, the square root of 11 is irrational.

Mahler’s Quinary Conundrum

The following question was posed by the German mathematician Kurt Mahler (1903 – 1988) in a letter that he wrote to Alf van der Poorten.

I am interested in the problem of whether there are squares of integers which, to the base g = 5, have only digits 0 or 1. I could not find a single example although I went quite far on my calculator.

Now some examples come immediately to mind. If \( n = 5^k\) then \( n^2 = 5^{2k}\) which is expressed in base 5 as \( 1000\ldots0\) with \( 2k\) zeros. But these solutions are trivial. We are more interested in solutions that are not divisible by 5, since they generate all other solutions. If \( n^2\) has only digits 0 or 1, then the same is true for \( (5^k \cdot n)^2\) and conversely. (The word “digits” always refers to base 5 digits in this note.)

I wrote a Python script to search for examples. The first step is to write a function to check if a number has only digits 0 or 1 in a given base.

def check(n, base):
    while n > 0:
        if n % base > 1:
            return False
        n = n // base
    return True

We could use this function to search for examples directly, but this is inefficient.

def naive_search(N):
    for n in range(N):
        if n % 5 > 0 and check(n*n, 5):
            yield n

It is inefficient because we are checking many values that could have been excluded from the outset. For example, the last digit of n must be 1 or 4, otherwise the last digit of \( n^2\) would be 4. We can find similar conditions on the last two digits, the last three digits, and so on.

  • The last two digits of n must be 01, 14, 31, or 44.
  • The last three digits of n must be 001, 031, 114, 144, 301, 331, 414, or 444.
  • The last four digits of n must be 0001, 0301, 0331, 0414, 1444, 2031, 2114, 2144, 2301, 2331, 2414, 3001, 4031, 4114, 4144, or 4444.

In fact, the number of possible \( k\)-digit endings is \( 2^k\). Each \( k\)-digit ending gives birth to two \( (k+1)\)-digit endings by adding digits to the left. I will leave it to the reader to prove this fact.

Here is a Python function that generates the possible \( k\)-digit endings. It works recursively by generating and extending the possible \( (k-1)\)-digit endings.

def gen(k):
    if k == 1:
        yield 1
        yield 4
    else:
        P = 5 ** (k-1)
        for n in gen(k-1):
            x = (n*n) // P
            delta = 2 * (n % 5) % 5
            d = x % 5
            m = n
            for a in range(5):
                if d < 2:
                    yield m
                d = (d + delta) % 5
                m += P        

The final step is to feed these candidates into our digit checker. By this means, I was able to find all positive integers n less than \( 5^{32} \approx 2.3 \times 10^{22}\) and not divisible by 5 such that \( n^2\) has only digits 0 or 1 when written in base 5.

results = sorted([n for n in gen(32) 
                    if check(n*n, 5)])

Results: 1, 972799, 3051273374, 6132750376, 839228035909, 3818357814376, 2384643515634376, 1490173338867234376, 931329727148437734376.

If we write these numbers in base 5, a pattern becomes apparent. Here is a Python script that converts the numbers to base 5.

def toBase5(n):
    if n < 5: 
        return n
    return 10 * toBase5(n//5) + (n % 5)

for n in results:
    print (toBase5(n))

# Base 10 Base 5
1 1 1
2 972799 222112144
3 3051273374 22222111221444
4 6132750376 100024441003001
5 839228035909 102222214334122114
6 3818357814376 1000024444100030001
7 2384643515634376 10000024444410000300001
8 1490173338867234376 100000024444441000003000001
9 931329727148437734376 1000000024444444100000030000001

One sees a similarity between the fourth, sixth, seventh, eighth, and ninth terms of this sequence. The runs of 0s and 4s increase by one in each step. This can be described by a polynomial function. If \( P(x) = 25x^4 + 15x^3 – 4x^2 + 3x + 1\) then these terms are \( P(5^3), P(5^4), P(5^5), P(5^6)\), and \( P(5^7)\).

We can verify that

$$
\begin{align*}
[P(x)]^2 & = 625x^8+750x^7+25x^6+30x^5 \\
& + 156x^4+6x^3+x^2+6x+1.
\end{align*}
$$

Note that the coefficients are positive integers less than \( 5^5\) and they can be written in base 5 using only 0s and 1s. Therefore, \( (P(5^k))^2\) has only digits 0 or 1 for all \( k \ge 3\).

This gives a striking, albeit partial, answer to Mahler’s question. There exist infinitely many squares that are not divisible by 5, but which have only the digits 0 or 1 when written in base 5.

Credits: Thanks to Gary Davis (@republicofmath) for bringing this problem to my attention. The problem is mentioned in the article The Legacy of Kurt Mahler, which appears in the May 2015 edition of the Notices of the American Mathematical Society.