Vietnamese math puzzle

Alex Bellos shared the following math puzzle that was given to Vietnamese eight-year-olds. The challenge is to place the digits from 1 to 9 into the grid.

Vietnamese math puzzle

Preliminary remarks

  1. Although it wasn’t stated explicitly, I assume that each digit should be used only once.
  2. The colon signifies division.
  3. There are many correct answers. I assume that students were meant to use enlightened trial and error to find a solution, but were not expected to find all solutions.
  4. I am assuming the standard order of operations, but I am not sure if this was the original intention. Would an eight-year-old be expected to understand that multiplication and division have higher precedence than addition and subtraction?
  5. Frankly, this puzzle is kind of boring for grown-ups.

Python code

Searching for solutions by hand did not appeal to me, so I wrote a Python program to search for all solutions. I hope that this example is helpful for people who are learning Python. (Perhaps it will inspire someone to write a better program than this.)

from itertools import permutations
standard_order = True
tolerance = 1e-7
target = 66

if standard_order:
    expr = "? + 13 * ? / ? + ? + 12 * ? - ? - 11 + ? * ? / ? - 10"
else:
    expr = "(((? + 13) * ? / ? + ? + 12) * ? - ? - 11 + ?) * ? / ? - 10"
expr = expr.replace("?", "%d")

for inputs in permutations(range(1,10)):
    expression = expr % inputs
    result = eval(expression)
    if abs(result - target) < tolerance:
        print ("%s = %s" % (expression, target))

Discussion

The itertools package provides many useful tools for iteration, including the permutation function, which generates all permutations of a list. Our script iterates through all 362880 permutations of the set of integers from 1 to 9.

Each permutation is substituted into the expression and evaluated. The program uses the eval function, which allows Python to run Python code within itself. This can be very dangerous, so use this function with care.

The division operations can lead to round-off error, so we only check that the result is very close to the target value. The variable tolerance sets the required precision.

Output

There are 136 solutions.

1 + 13 * 2 / 6 + 4 + 12 * 7 - 8 - 11 + 3 * 5 / 9 - 10 = 66
1 + 13 * 2 / 6 + 4 + 12 * 7 - 8 - 11 + 5 * 3 / 9 - 10 = 66
1 + 13 * 3 / 2 + 4 + 12 * 5 - 8 - 11 + 7 * 9 / 6 - 10 = 66
1 + 13 * 3 / 2 + 4 + 12 * 5 - 8 - 11 + 9 * 7 / 6 - 10 = 66
1 + 13 * 3 / 2 + 9 + 12 * 5 - 6 - 11 + 4 * 7 / 8 - 10 = 66
1 + 13 * 3 / 2 + 9 + 12 * 5 - 6 - 11 + 7 * 4 / 8 - 10 = 66
1 + 13 * 3 / 4 + 7 + 12 * 6 - 5 - 11 + 2 * 9 / 8 - 10 = 66
1 + 13 * 3 / 4 + 7 + 12 * 6 - 5 - 11 + 9 * 2 / 8 - 10 = 66
1 + 13 * 3 / 6 + 2 + 12 * 7 - 9 - 11 + 4 * 5 / 8 - 10 = 66
1 + 13 * 3 / 6 + 2 + 12 * 7 - 9 - 11 + 5 * 4 / 8 - 10 = 66
1 + 13 * 3 / 9 + 4 + 12 * 7 - 8 - 11 + 2 * 5 / 6 - 10 = 66
1 + 13 * 3 / 9 + 4 + 12 * 7 - 8 - 11 + 5 * 2 / 6 - 10 = 66
1 + 13 * 4 / 8 + 2 + 12 * 7 - 9 - 11 + 3 * 5 / 6 - 10 = 66
1 + 13 * 4 / 8 + 2 + 12 * 7 - 9 - 11 + 5 * 3 / 6 - 10 = 66
1 + 13 * 5 / 2 + 3 + 12 * 4 - 8 - 11 + 7 * 9 / 6 - 10 = 66
1 + 13 * 5 / 2 + 3 + 12 * 4 - 8 - 11 + 9 * 7 / 6 - 10 = 66
1 + 13 * 5 / 2 + 8 + 12 * 4 - 7 - 11 + 3 * 9 / 6 - 10 = 66
1 + 13 * 5 / 2 + 8 + 12 * 4 - 7 - 11 + 9 * 3 / 6 - 10 = 66
1 + 13 * 5 / 3 + 9 + 12 * 4 - 2 - 11 + 7 * 8 / 6 - 10 = 66
1 + 13 * 5 / 3 + 9 + 12 * 4 - 2 - 11 + 8 * 7 / 6 - 10 = 66
1 + 13 * 8 / 3 + 7 + 12 * 4 - 5 - 11 + 2 * 6 / 9 - 10 = 66
1 + 13 * 8 / 3 + 7 + 12 * 4 - 5 - 11 + 6 * 2 / 9 - 10 = 66
1 + 13 * 9 / 6 + 4 + 12 * 5 - 8 - 11 + 3 * 7 / 2 - 10 = 66
1 + 13 * 9 / 6 + 4 + 12 * 5 - 8 - 11 + 7 * 3 / 2 - 10 = 66
1 + 13 * 9 / 6 + 7 + 12 * 5 - 2 - 11 + 3 * 4 / 8 - 10 = 66
1 + 13 * 9 / 6 + 7 + 12 * 5 - 2 - 11 + 4 * 3 / 8 - 10 = 66
2 + 13 * 1 / 4 + 3 + 12 * 7 - 9 - 11 + 5 * 6 / 8 - 10 = 66
2 + 13 * 1 / 4 + 3 + 12 * 7 - 9 - 11 + 6 * 5 / 8 - 10 = 66
2 + 13 * 3 / 6 + 1 + 12 * 7 - 9 - 11 + 4 * 5 / 8 - 10 = 66
2 + 13 * 3 / 6 + 1 + 12 * 7 - 9 - 11 + 5 * 4 / 8 - 10 = 66
2 + 13 * 4 / 8 + 1 + 12 * 7 - 9 - 11 + 3 * 5 / 6 - 10 = 66
2 + 13 * 4 / 8 + 1 + 12 * 7 - 9 - 11 + 5 * 3 / 6 - 10 = 66
2 + 13 * 6 / 9 + 8 + 12 * 5 - 1 - 11 + 4 * 7 / 3 - 10 = 66
2 + 13 * 6 / 9 + 8 + 12 * 5 - 1 - 11 + 7 * 4 / 3 - 10 = 66
2 + 13 * 8 / 6 + 9 + 12 * 4 - 1 - 11 + 5 * 7 / 3 - 10 = 66
2 + 13 * 8 / 6 + 9 + 12 * 4 - 1 - 11 + 7 * 5 / 3 - 10 = 66
2 + 13 * 9 / 6 + 3 + 12 * 5 - 1 - 11 + 4 * 7 / 8 - 10 = 66
2 + 13 * 9 / 6 + 3 + 12 * 5 - 1 - 11 + 7 * 4 / 8 - 10 = 66
3 + 13 * 1 / 4 + 2 + 12 * 7 - 9 - 11 + 5 * 6 / 8 - 10 = 66
3 + 13 * 1 / 4 + 2 + 12 * 7 - 9 - 11 + 6 * 5 / 8 - 10 = 66
3 + 13 * 2 / 1 + 5 + 12 * 4 - 7 - 11 + 8 * 9 / 6 - 10 = 66
3 + 13 * 2 / 1 + 5 + 12 * 4 - 7 - 11 + 9 * 8 / 6 - 10 = 66
3 + 13 * 2 / 4 + 8 + 12 * 5 - 1 - 11 + 7 * 9 / 6 - 10 = 66
3 + 13 * 2 / 4 + 8 + 12 * 5 - 1 - 11 + 9 * 7 / 6 - 10 = 66
3 + 13 * 2 / 8 + 6 + 12 * 5 - 1 - 11 + 7 * 9 / 4 - 10 = 66
3 + 13 * 2 / 8 + 6 + 12 * 5 - 1 - 11 + 9 * 7 / 4 - 10 = 66
3 + 13 * 5 / 2 + 1 + 12 * 4 - 8 - 11 + 7 * 9 / 6 - 10 = 66
3 + 13 * 5 / 2 + 1 + 12 * 4 - 8 - 11 + 9 * 7 / 6 - 10 = 66
3 + 13 * 6 / 4 + 9 + 12 * 5 - 8 - 11 + 1 * 7 / 2 - 10 = 66
3 + 13 * 6 / 4 + 9 + 12 * 5 - 8 - 11 + 7 * 1 / 2 - 10 = 66
3 + 13 * 9 / 2 + 8 + 12 * 1 - 5 - 11 + 6 * 7 / 4 - 10 = 66
3 + 13 * 9 / 2 + 8 + 12 * 1 - 5 - 11 + 7 * 6 / 4 - 10 = 66
3 + 13 * 9 / 6 + 2 + 12 * 5 - 1 - 11 + 4 * 7 / 8 - 10 = 66
3 + 13 * 9 / 6 + 2 + 12 * 5 - 1 - 11 + 7 * 4 / 8 - 10 = 66
4 + 13 * 2 / 6 + 1 + 12 * 7 - 8 - 11 + 3 * 5 / 9 - 10 = 66
4 + 13 * 2 / 6 + 1 + 12 * 7 - 8 - 11 + 5 * 3 / 9 - 10 = 66
4 + 13 * 3 / 2 + 1 + 12 * 5 - 8 - 11 + 7 * 9 / 6 - 10 = 66
4 + 13 * 3 / 2 + 1 + 12 * 5 - 8 - 11 + 9 * 7 / 6 - 10 = 66
4 + 13 * 3 / 9 + 1 + 12 * 7 - 8 - 11 + 2 * 5 / 6 - 10 = 66
4 + 13 * 3 / 9 + 1 + 12 * 7 - 8 - 11 + 5 * 2 / 6 - 10 = 66
4 + 13 * 9 / 6 + 1 + 12 * 5 - 8 - 11 + 3 * 7 / 2 - 10 = 66
4 + 13 * 9 / 6 + 1 + 12 * 5 - 8 - 11 + 7 * 3 / 2 - 10 = 66
5 + 13 * 1 / 2 + 9 + 12 * 6 - 7 - 11 + 3 * 4 / 8 - 10 = 66
5 + 13 * 1 / 2 + 9 + 12 * 6 - 7 - 11 + 4 * 3 / 8 - 10 = 66
5 + 13 * 2 / 1 + 3 + 12 * 4 - 7 - 11 + 8 * 9 / 6 - 10 = 66
5 + 13 * 2 / 1 + 3 + 12 * 4 - 7 - 11 + 9 * 8 / 6 - 10 = 66
5 + 13 * 3 / 1 + 7 + 12 * 2 - 6 - 11 + 8 * 9 / 4 - 10 = 66
5 + 13 * 3 / 1 + 7 + 12 * 2 - 6 - 11 + 9 * 8 / 4 - 10 = 66
5 + 13 * 4 / 1 + 9 + 12 * 2 - 7 - 11 + 3 * 8 / 6 - 10 = 66
5 + 13 * 4 / 1 + 9 + 12 * 2 - 7 - 11 + 8 * 3 / 6 - 10 = 66
5 + 13 * 4 / 8 + 9 + 12 * 6 - 7 - 11 + 1 * 3 / 2 - 10 = 66
5 + 13 * 4 / 8 + 9 + 12 * 6 - 7 - 11 + 3 * 1 / 2 - 10 = 66
5 + 13 * 7 / 2 + 8 + 12 * 3 - 9 - 11 + 1 * 6 / 4 - 10 = 66
5 + 13 * 7 / 2 + 8 + 12 * 3 - 9 - 11 + 6 * 1 / 4 - 10 = 66
5 + 13 * 9 / 3 + 6 + 12 * 2 - 1 - 11 + 7 * 8 / 4 - 10 = 66
5 + 13 * 9 / 3 + 6 + 12 * 2 - 1 - 11 + 8 * 7 / 4 - 10 = 66
6 + 13 * 2 / 8 + 3 + 12 * 5 - 1 - 11 + 7 * 9 / 4 - 10 = 66
6 + 13 * 2 / 8 + 3 + 12 * 5 - 1 - 11 + 9 * 7 / 4 - 10 = 66
6 + 13 * 3 / 1 + 9 + 12 * 2 - 5 - 11 + 7 * 8 / 4 - 10 = 66
6 + 13 * 3 / 1 + 9 + 12 * 2 - 5 - 11 + 8 * 7 / 4 - 10 = 66
6 + 13 * 9 / 3 + 5 + 12 * 2 - 1 - 11 + 7 * 8 / 4 - 10 = 66
6 + 13 * 9 / 3 + 5 + 12 * 2 - 1 - 11 + 8 * 7 / 4 - 10 = 66
7 + 13 * 1 / 4 + 9 + 12 * 6 - 5 - 11 + 2 * 3 / 8 - 10 = 66
7 + 13 * 1 / 4 + 9 + 12 * 6 - 5 - 11 + 3 * 2 / 8 - 10 = 66
7 + 13 * 2 / 8 + 9 + 12 * 6 - 5 - 11 + 1 * 3 / 4 - 10 = 66
7 + 13 * 2 / 8 + 9 + 12 * 6 - 5 - 11 + 3 * 1 / 4 - 10 = 66
7 + 13 * 3 / 1 + 5 + 12 * 2 - 6 - 11 + 8 * 9 / 4 - 10 = 66
7 + 13 * 3 / 1 + 5 + 12 * 2 - 6 - 11 + 9 * 8 / 4 - 10 = 66
7 + 13 * 3 / 2 + 8 + 12 * 5 - 9 - 11 + 1 * 6 / 4 - 10 = 66
7 + 13 * 3 / 2 + 8 + 12 * 5 - 9 - 11 + 6 * 1 / 4 - 10 = 66
7 + 13 * 3 / 4 + 1 + 12 * 6 - 5 - 11 + 2 * 9 / 8 - 10 = 66
7 + 13 * 3 / 4 + 1 + 12 * 6 - 5 - 11 + 9 * 2 / 8 - 10 = 66
7 + 13 * 5 / 2 + 8 + 12 * 4 - 9 - 11 + 1 * 3 / 6 - 10 = 66
7 + 13 * 5 / 2 + 8 + 12 * 4 - 9 - 11 + 3 * 1 / 6 - 10 = 66
7 + 13 * 6 / 4 + 8 + 12 * 5 - 9 - 11 + 1 * 3 / 2 - 10 = 66
7 + 13 * 6 / 4 + 8 + 12 * 5 - 9 - 11 + 3 * 1 / 2 - 10 = 66
7 + 13 * 8 / 3 + 1 + 12 * 4 - 5 - 11 + 2 * 6 / 9 - 10 = 66
7 + 13 * 8 / 3 + 1 + 12 * 4 - 5 - 11 + 6 * 2 / 9 - 10 = 66
7 + 13 * 9 / 6 + 1 + 12 * 5 - 2 - 11 + 3 * 4 / 8 - 10 = 66
7 + 13 * 9 / 6 + 1 + 12 * 5 - 2 - 11 + 4 * 3 / 8 - 10 = 66
8 + 13 * 2 / 4 + 3 + 12 * 5 - 1 - 11 + 7 * 9 / 6 - 10 = 66
8 + 13 * 2 / 4 + 3 + 12 * 5 - 1 - 11 + 9 * 7 / 6 - 10 = 66
8 + 13 * 3 / 2 + 7 + 12 * 5 - 9 - 11 + 1 * 6 / 4 - 10 = 66
8 + 13 * 3 / 2 + 7 + 12 * 5 - 9 - 11 + 6 * 1 / 4 - 10 = 66
8 + 13 * 5 / 2 + 1 + 12 * 4 - 7 - 11 + 3 * 9 / 6 - 10 = 66
8 + 13 * 5 / 2 + 1 + 12 * 4 - 7 - 11 + 9 * 3 / 6 - 10 = 66
8 + 13 * 5 / 2 + 7 + 12 * 4 - 9 - 11 + 1 * 3 / 6 - 10 = 66
8 + 13 * 5 / 2 + 7 + 12 * 4 - 9 - 11 + 3 * 1 / 6 - 10 = 66
8 + 13 * 6 / 4 + 7 + 12 * 5 - 9 - 11 + 1 * 3 / 2 - 10 = 66
8 + 13 * 6 / 4 + 7 + 12 * 5 - 9 - 11 + 3 * 1 / 2 - 10 = 66
8 + 13 * 6 / 9 + 2 + 12 * 5 - 1 - 11 + 4 * 7 / 3 - 10 = 66
8 + 13 * 6 / 9 + 2 + 12 * 5 - 1 - 11 + 7 * 4 / 3 - 10 = 66
8 + 13 * 7 / 2 + 5 + 12 * 3 - 9 - 11 + 1 * 6 / 4 - 10 = 66
8 + 13 * 7 / 2 + 5 + 12 * 3 - 9 - 11 + 6 * 1 / 4 - 10 = 66
8 + 13 * 9 / 2 + 3 + 12 * 1 - 5 - 11 + 6 * 7 / 4 - 10 = 66
8 + 13 * 9 / 2 + 3 + 12 * 1 - 5 - 11 + 7 * 6 / 4 - 10 = 66
9 + 13 * 1 / 2 + 5 + 12 * 6 - 7 - 11 + 3 * 4 / 8 - 10 = 66
9 + 13 * 1 / 2 + 5 + 12 * 6 - 7 - 11 + 4 * 3 / 8 - 10 = 66
9 + 13 * 1 / 4 + 7 + 12 * 6 - 5 - 11 + 2 * 3 / 8 - 10 = 66
9 + 13 * 1 / 4 + 7 + 12 * 6 - 5 - 11 + 3 * 2 / 8 - 10 = 66
9 + 13 * 2 / 8 + 7 + 12 * 6 - 5 - 11 + 1 * 3 / 4 - 10 = 66
9 + 13 * 2 / 8 + 7 + 12 * 6 - 5 - 11 + 3 * 1 / 4 - 10 = 66
9 + 13 * 3 / 1 + 6 + 12 * 2 - 5 - 11 + 7 * 8 / 4 - 10 = 66
9 + 13 * 3 / 1 + 6 + 12 * 2 - 5 - 11 + 8 * 7 / 4 - 10 = 66
9 + 13 * 3 / 2 + 1 + 12 * 5 - 6 - 11 + 4 * 7 / 8 - 10 = 66
9 + 13 * 3 / 2 + 1 + 12 * 5 - 6 - 11 + 7 * 4 / 8 - 10 = 66
9 + 13 * 4 / 1 + 5 + 12 * 2 - 7 - 11 + 3 * 8 / 6 - 10 = 66
9 + 13 * 4 / 1 + 5 + 12 * 2 - 7 - 11 + 8 * 3 / 6 - 10 = 66
9 + 13 * 4 / 8 + 5 + 12 * 6 - 7 - 11 + 1 * 3 / 2 - 10 = 66
9 + 13 * 4 / 8 + 5 + 12 * 6 - 7 - 11 + 3 * 1 / 2 - 10 = 66
9 + 13 * 5 / 3 + 1 + 12 * 4 - 2 - 11 + 7 * 8 / 6 - 10 = 66
9 + 13 * 5 / 3 + 1 + 12 * 4 - 2 - 11 + 8 * 7 / 6 - 10 = 66
9 + 13 * 6 / 4 + 3 + 12 * 5 - 8 - 11 + 1 * 7 / 2 - 10 = 66
9 + 13 * 6 / 4 + 3 + 12 * 5 - 8 - 11 + 7 * 1 / 2 - 10 = 66
9 + 13 * 8 / 6 + 2 + 12 * 4 - 1 - 11 + 5 * 7 / 3 - 10 = 66
9 + 13 * 8 / 6 + 2 + 12 * 4 - 1 - 11 + 7 * 5 / 3 - 10 = 66

Russian peasant multiplication and loop invariants

Russian peasant multiplication is a method for multiplying numbers based on repeated doubling and halving. Nobody is quite sure why it is called “Russian peasant multiplication,” but the name has stuck.

Here is the procedure:

  1. Write each number at the head of a column.
  2. Double the number in the first column.
  3. Halve the number in the second column. (If odd, ignore the remainder.)
  4. Repeat (2) and (3) until the number in the second column is 1.
  5. Cross out the rows with an even number in the second column.
  6. Add up the remaining numbers in the first column.

Example: Multiply 43 by 18.

First number Second number
 43   18 
\(86\) \(9\)
 172   4 
 344   2 
\(688\) \(1\)

\(\implies 43 \times 18 = 86 + 688 = 774.\)

Here is a Python implementation of the algorithm.

def multiply(m, n):
    p = 0
    while n > 0:
        if n % 2: # n is odd
            p = p + m
        m = m * 2
        n = n // 2
    return p

Correctness

Let’s see what happens to the values of \(m\), \(n\), and \(p\) in one iteration. If \(n\) is even, then \(m\) is doubled and \(n\) is halved, so the product \(mn\) is unchanged. If \(n\) is odd, then \(m\) is doubled and \(n\) is replaced with \((n-1)/2\), so the value of the product \(mn\) is decreased.
\[(2m) \cdot (n-1)/2 = mn – m\]
To compensate for this decrease, the value of \(p\) is increased by \(m\) when \(n\) is odd. Consequently, the value of \(Q = mn + p\) does not change – it is a loop invariant. At the beginning of the loop we have \(Q = mn\) (since \(p=0\)) and at the end we have \(Q = p\) (since \(n=0\)). Therefore, the product \(mn\) is equal to the final value of \(p\).

A decimal variation

The Russian Peasant method is associated with binary (base two) arithmetic, since it involves multiplying and dividing by two. Is there a decimal (base ten) version? Consider the following:

def multiply_decimal(m, n):
    p = 0
    while n > 0:
        r = n % 10
        p = p + m * r
        m = m * 10
        n = n // 10
    return p

Let’s analyze this algorithm. Each time through the loop, the values of \(m\), \(n\), and \(p\) are replaced with new values \(m’\), \(n’\), and \(p’\). The following relations hold between the new values and the old values.
\[
\begin{align*}
m’ &= 10m\\
n &= 10n’ + r\\
p’ &= p + mr
\end{align*}
\]

As before, the quantity \(Q = mn+p\) is a loop invariant, since
\(
mn + p = m(10n’ + r) + p = (10m)n’ + (p + mr) = m’n’ + p’.
\)
Therefore, the function returns the product \(mn\).

Example: Multiply 123 by 456.

\(m\) \(n\) \(mr\)
\(123\) \(456\) \(123 \times 6 = 738\)
\(1230\) \(45\) \(1230 \times 5 = 6150\)
\(12300\) \(4\) \(12300 \times 4 = 49200\)

\(\implies 123 \times 456 = 738 + 6150 + 49200 = 56088.\)

Note that this is essentially the standard algorithm for multiplying numbers.

  123
× 456
—————
  738
 6150
49200
—————
56088

Fast exponentiation

We have seen that we multiply numbers via repeated doubling. In a similar way, we can raise a number to a power via repeated squaring. Here is a Python implementation:

def power(m, n):
    p = 1
    while n > 0:
       if n % 2: # n is odd
            p = p * m
        m = m * m
        n = n // 2
    return p

The loop invariant \(Q = m^n p\) can be used to prove that the function correctly returns the value of \(m^n\).

Emoji and Zipf’s Law

Emojitracker.com is a very interesting site that tracks how often each emoji is used on Twitter. The statistics are updated in real time, and it’s really quite amazing to watch.

Screenshot from 2015-05-06 23:40:09

(In case the reader is unaware, an emoji is a small icon that is often used to convey ideas or emotions in electronic communications. It should not be confused with an emoticon, which uses ordinary text characters for a similar purpose.)

Is there is a formula that describes these frequencies? In 1935, the American linguist George Zipf proposed that, in any corpus of words, the frequency of a word is inversely proportional to its rank in the frequency table. That is, the most common word should be \(N\) times more frequent than the \(N\)th word. More generally, we may suppose that the frequencies follow an inverse power law of the form \(y = Cx^{-r}\), where \(y\) is the frequency and \(x\) is the rank.

Power laws can be identified using a log-log plot, which has logarithmic scales on both axes. When we take the logarithm of both sides, we get

\[\log y = \log C – r \log x\]

which appears on a log-log plot as a straight line with slope \((-r)\) .

Does the frequency of emojis follow a power law? To answer this question, we need to find a way to scrape the data from the webpage. My solution was crude but effective. I selected all of the text on the page, and used copy-and-paste to get it into a text file. I wrote a short Python script to extract the numbers into a list.

import re
with open("emoji.txt") as f:
    text = '\n'.join(f.readlines())
y = [int(m) for m in re.findall(r"\d+", text)]
x = range(1, len(y)+1)

Next, I used matplotlib to plot the data on a log-log scale, and scipy.stats to compute the line of best fit for the transformed data.

from math import log, exp
from scipy import stats
import matplotlib.pyplot as plt

logx = map(log, x)
logy = map(log, y)
slope, intercept, r_value, p_value, std_er = stats.linregress(logx[:40], logy[:40])
y_pred = [exp(slope*t + intercept) for t in logx]
plt.loglog(x, y)
plt.loglog(x, y_pred)
plt.xlabel('Rank')
plt.ylabel('Frequency')
plt.title('Frequency of Emoji characters on Twitter')
plt.grid(True)
plt.text(100, 3e7, r'$y = (%.3f \times 10^8)\ x^{%.3f}$' % (exp(intercept)/10**8, slope)) 
plt.text(100, 1e8, r'$(r^2 = \ %.4f)$' % r_value**2)
plt.savefig('emoji_zipf.png')
plt.show()

From the graph shown below, it appears that the numbers do not follow a power law. However, the 40 most common emojis do appear to follow a power law. The graph shows the frequencies for all of the emojis, with a line of best fit for the 40 most common emojis.

emoji_zipf

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.

Intersecting secants and coordinate geometry

The intersecting secants theorem states that when two secant lines intersect each other, the products of their segments are equal. In the picture below,

\[ PA \cdot PB = PC \cdot PD. \]

Secant lines AB and CD intersecting at a point P outside the circle.

The standard proof is as follows. Angles PAD and PCB are congruent because they are inscribed angles which cut off the same arc BD. The triangles PAD and PCB are similar because two pairs of corresponding angles are congruent. Therefore, PA/PD = PC/PB, hence \(PA \cdot PB = PC \cdot PD\).

We could also prove this result using coordinate geometry. Assume that P is at the origin, and suppose that the circle has equation \(x^2 + y^2 + ax + by + c = 0 \).

Any line through the origin can be parametrized as \(x = r\cos \theta\) and \(y = r\sin \theta\) where \(r\) is the distance from the origin and \(\theta\) is fixed.

Substituting into the equation of the circle yields
\[r^2 + ar \cos \theta + br \sin \theta + c = 0.\]

This is a quadratic equation in \(r\), and the product of the roots is \(c\). Therefore \(PA \cdot PB = PC \cdot PD = c\).