Sums of consecutive reciprocals

The sum of the reciprocals of consecutive integers is never an integer. That is, for all positive integers m and n with n > m, the sum

\sum_{k=m}^n \frac{1}{k}

is never an integer. This was proved by József Kürschák in 1908.

This means that the harmonic numbers defined by

H_n = \sum{k=1}^n \frac{1}{k}

are never integers for n > 1. The harmonic series diverges, so the sequence of harmonic numbers goes off to infinity, but it does so carefully avoiding all integers along the way.

Kürschák’s theorem says that not only are the harmonic numbers never integers, the difference of two distinct harmonic numbers is never an integer. That is, HnHm is never an integer unless m = n.

Related posts

Average sum of digits

The smaller the base you write numbers in, the smaller their digits sums will be on average.

This need not be true of any given number, only for numbers on average. For example, let n = 2021. In base 10, the sum of the digits is 5, but in base 2 the sum of the digits is 8 because n is written as 11111100101 in binary. So for our particular example, a smaller base lead to a larger digit sum. In fact, all smaller bases than 10 lead to a larger digit sum for 2021. But on average, numbers have smaller digit sums when written in binary than when written in base 10.

Let’s explore this with a little Python code, then look at a theoretical result.

First, we need some code to express numbers in different bases, and here’s a simple function for that.

    # From
    def numberToBase(n, b):
        if n == 0:
            return [0]
        digits = []
        while n:
            digits.append(int(n % b))
            n //= b
        return digits[::-1]

Using the code above, our digit sum function is trivial:

    def digitSum(n, b):
        return sum(numberToBase(n, b))

Now for our simulation.

    import numpy as np

    upper_bound = 1000000000
    sample_size = 1000
    bases = np.arange(2, 17)

    sample = np.random.randint(0, upper_bound, size = sample_size)
    out = [sum(digitSum(n, b) for n in sample)/sample_size for b in bases]

A theorem in [1] says that

\frac{S(b, N)}{N} \sim \frac{(r-1)\log N}{2\log b}

where S(b, N) is the sum of the number of digits of all numbers from 1 to N when written in base b. This value is asymptotically equal to the expression on the right.

We can compute the average number of digits we’d expect in our simulation using the asymptotic result as an approximation.

    asymp = [0.5*(b-1)*np.log(upper_bound) / np.log(b) for b in bases]

The table below compares the simulation results to theory.

    | base | simulation | asymptotic |
    |    2 |     14.899 |     14.949 |
    |    3 |     18.660 |     18.863 |
    |    4 |     22.317 |     22.423 |
    |    5 |     25.398 |     25.752 |
    |    6 |     28.345 |     28.915 |
    |    7 |     30.996 |     31.949 |
    |    8 |     34.839 |     34.880 |
    |    9 |     36.422 |     37.726 |
    |   10 |     40.275 |     40.500 |
    |   11 |     41.480 |     43.211 |
    |   12 |     44.353 |     45.868 |
    |   13 |     47.454 |     48.476 |
    |   14 |     49.949 |     51.041 |
    |   15 |     51.772 |     53.567 |
    |   16 |     53.610 |     56.058 |

[1] L. E. Bush. An Asymptotic Formula for the Average Sum of the Digits of Integers. The American Mathematical Monthly , Mar., 1940, Vol. 47, pp. 154-156.

A connected topology for the integers

You can define a topology on the positive integers by choosing as an open basis sets the series of the form an + b where a and b are relatively prime positive integers. Solomon Golumb defines this topology in [1] and proves that it is connected.

But that paper doesn’t prove that proposed basis really is a basis. That is implicitly an exercise for the reader, and we carry out that exercise here. The proof will say that certain things can be computed, and we’ll go a step further an actually compute them with Python.


To prove that these sets form the basis of a topology, we need to establish two things.

  1. Every point is contained in some basis element.
  2. The intersection of basis elements is another basis element or empty.

The first is easy. For any positive number x, the sequence (x+1)n + x contains x, and x+1 is relatively prime to x. Any other number relatively prime to x would work just as well.

The second is more interesting. Consider two sequences: an + b and cm + d. A number x in the intersection satisfies the pair of equations

x = b mod a
x = d mod c.

If a and c are relatively prime, the Chinese Remainder Theorem says that the equation has a solution y, and the solutions are unique mod ac. So the intersection of the two series is acn + y.

If a and c are not relatively prime, let r be their greatest common divisor. Then the intersection of an + b and cm + d is another sequence of the same form if b and d are congruent mod r and empty otherwise.

Separating points

A topological space is connected if it is not the union of two disjoint open sets. An easy way to make a space connected is to not have any disjoint open sets. But Golumb’s topology has plenty of disjoint open sets.

In fact, his topology is Hausdorff, meaning that any two points can be separated by disjoint open sets. The proof is simple. Take two points s and t. Let p be any prime bigger than both s and t, and consider the two sets pn + s and pn + t.

Python calculation

Our proof split into two cases: when the strides of the two series are relatively prime and when they are not.

Relatively prime strides

To get concrete, suppose we have the sequences

2, 9, 16, 23, …


3, 14, 25, 35, …

or in other words 7n + 2 and 11n + 3. According to the theory above, these sequences intersect because 7 and 11 are relatively prime, and the intersection should be of the form 77n + y, and we should be able to compute y from the Chinese Remainder Theorem. We will use the function crt from SymPy. (See another example of using this function here.)

    >>> from sympy.ntheory.modular import crt
    >>> crt([7,11], [2, 3], symmetric=False)
    >>> (58, 77)

This reports that y = 58.

Now let’s verify that the intersection of our two series looks like 77n + 58.

    >>> A = set(2+7*n for n in range(100))
    >>> B = set(3+11*n for n in range(100))
    >>> sorted(A.intersection(B))
    [58, 135, 212, 289, 366, 443, 520, 597, 674]

Strides with a common factor: empty intersection

Now for the case where a and c have greatest common divisor r > 1. Consider, for example,

1, 16, 31, 46, …


2, 14, 26, 38, …

The sequences 15n + 1 and 12m + 2 do not intersect. The greatest common divisor of 15 and 12 is 3. Numbers in the first series are congruent to 1 mod 3 and numbers in the second series are congruent to 2 mod 3, so no number can be in both. If we didn’t realize this, we could call

    crt([15,12], [1, 2], symmetric=False)

and see that it returns None.

Strides with a common factor: finding the intersection

But now let’s look at

1, 16, 31, 46, …


13, 25, 37, 49, …

Now both sequences are congruent to 1 mod 3, so it’s at least feasible that the two series may intersect. And in fact they do.

    >>> crt([15,12], [1, 13], symmetric=False)
    (1, 60)

This says that the set of solutions are all congruent to 1 mod 60. Now 1 itself is not a solution, but 61 is, and so is 60n + 1 for all positive integers n. We can illustrate this as before by computing the intersection of the two series.

    >>> A = set(1+15*n for n in range(30))
    >>> B = set(13+12*n for n in range(30))
    >>> sorted(A.intersection(B))
    [61, 121, 181, 241, 301, 361]

Related posts

[1] Solomon W. Golomb. A Connected Topology for the Integers. The American Mathematical Monthly , Oct., 1959, pp. 663-665

Factors of consecutive products

Pick a positive integer k and take the product of k consecutive integers greater than k. Then the result is divisible by a prime number greater than k. This theorem was first proved 128 years ago [1].

For example, suppose we pick k = 5. If we take the product


then it’s obviously divisible by a prime larger than 5 because 23 is one of the factors. But we could also take, for example,


Although the sequence {32, 33, 34, 35, 36} doesn’t contain any primes, it contains a couple numbers with prime factors larger than 5, namely 33 = 3*11 and 34 = 2*17. Sylvester’s theorem guarantees that there will always be one number out of 5 consecutive integers that is divisible by a prime greater than 5, provided the first number is itself greater than 5.

For one more example, we can pick k = 8 and look at


One of the numbers between 140 and 147 (inclusive) must be divisible by a prime greater than 8, and in fact five of them are, starting with 141 = 3*47.

In both of our examples there were multiple numbers that had a prime factor large enough to satisfy Sylvester’s theorem. Are there examples where only one number has a factor large enough? Yes, but not many.

For example, {8, 9} is a pair of consecutive integers, only one of which has a prime factor bigger than 2. And {8, 9, 10} is an example of 3 consecutive integers, only one of which has a prime factor larger than 3. I wrote a program to search for examples with k = 4 and didn’t find any.

In any case, a theorem dating back to 1897 [2] proves that there could only be finitely many examples for k > 2. Lehmer extends the results [2] in his paper [3].


[1] J. J. Sylvester, On arithmetical series, Messenger Math., 21 (1892) 1-19, and 87-120; Collected Mathematical Papers, 4 (1912) 687-731.

[2] G. St0rmer, Quelques theoremes sur l’equation de Pell x² – Dy² = ± 1 et leurs applications,
Skr. Norske Vid. Akad. Oslo, I no. 2 (1897).

[3] D. H. Lehmer, The Prime Factors of Consecutive Integers, The American Mathematical Monthly , Feb., 1965, Vol. 72, No. 2, Part 2, pp. 19-20.

Ruling out gaps in the list of Mersenne primes

The largest known primes are Mersenne primes, in part because the Lucas-Lehmer algorithm makes it more efficient to test Mersenne numbers for primality than to test arbitrary numbers. These numbers have the form 2n – 1.

Trend in Mersenne primes

The Great Internet Mersenne Prime Search (GIMPS) announced today that they have tested whether 2n – 1 is prime for all values of n less than 100,000,000 [1]. More specifically, they’ve tested each exponent at least once, though they’d like to test each exponent twice before giving up on it.

If we assume all of their initial calculations are correct, we can say, for example, that 282,589,933 – 1 is the 51st Mersenne prime. Before we had to say that it was the 51st known Mersenne prime because there could be a smaller, undiscovered Mersenne prime, in which case 282,589,933 – 1 might be the 52nd Mersenne prime. Mersenne primes have not always been discovered in order. For example, the 29th Mersenne prime was discovered after the 30th and 31st such numbers.

So for now, again assuming all the GIMPS calculations have been correct, we know all Mersenne primes less than 2100,000,000, and there are 51 of them.

Related posts

[1] There’s no need to test every value of n. We know that if 2n – 1 is prime, then n is prime. So we only need to check prime values of n.

Density of square-free numbers, cube-free numbers, etc.

An integer n is said to be square-free if no perfect square (other than 1) divides n. Similarly, n is said to be cube-free if no perfect cube (other than 1) divides n.

In general, n is said to be k-free if it has no divisor d > 1 where d is a kth power [1].

As x gets larger, the proportion of numbers less than x that are k-free approaches 1/ζ(k). That is

\lim_{x\to\infty} \frac{Q_k(x)}{x} = \frac{1}{\zeta(k)}

where Qk(x) is the number of k-free positive integers less than or equal to x and ζ is the Riemann zeta function.

The values of ζ at even integers are known in closed form. So, for example, the proportion of square-free numbers less than x approaches 6/π² since ζ(2) = π²/6. For odd integers the values of ζ must be calculated numerically.

Once you have a theorem that says one thing converges to another, it’s natural to ask next is how fast it converges. How well does the limit work as an approximation of the finite terms? Are there any bounds on the errors?

In [2] the authors give new estimates the remainder term

R_k(x) = Q_k(x) - \frac{x}{\zeta(k)}

Python code

Let’s play with this using Python. The code below is meant to be clear rather than efficient.

    from sympy import factorint

    def kfree(n, k):
        exponents = factorint(n).values()
        return max(exponents) < k

    def Q(k, x):
        q = 1 # because 1 is not k-free
        for n in range(2, x+1):
            if kfree(n, k):
                q += 1
        return q

    print(Q(4, 1000000))

This says Q4(1,000,000) = 923,939.

So Q4(106)/106 = 0.923939. Now ζ(4) = π4/90, so 1/ζ(4) = 0.9239384….

We have R4(106) = 0.597 and so in this case R is quite small.

(If we didn’t know ζ in closed form, we could calculate it in python with scipy.special.zeta.)

More posts involving Riemann zeta

[1] It would be clearer if we said kth power-free, but people shorten this to k-free.

In my opinion the term “k-free” is misleading. The term ought to mean a number not divisible by k, or maybe a number that has no digit k when written in some base. But the terminology is what it is.

[2] Michael Mossinghoff, Tomás Oliviera e Sliva, and Timothy Trudgian. The distribution of k-free numbers. Mathematics of Computation. Posted November 24, 2020.

Digitally delicate primes

A digitally delicate prime is a prime number such that if any of its digits are changed by any amount, the result is not prime. The smallest such number is 294001. This means that variations on this number such as 394001 or 291001 or 294081 are all composite.

A paper [1] posted last week gives several new results concerning digitally delicate primes. Here I’ll just mention a few things that are known.

First, there are infinitely many digitally delicate primes. They are uncommon, but they have positive density. That is, as x goes to infinity, the ratio of the number of digitally delicate primes less than x to the total number of primes less than x converges to a small positive number.

This result has been extended to bases other than base 10. In [1] the authors extend the result to widely digitally delicate primes, meaning that the result holds even if you consider the “digits” of a number to include leading zeros.

The following Python code tests whether a number is a digitally delicate prime. It then takes the first six digitally delicate primes and verifies that they are indeed digitally delicate.

from sympy import isprime
from math import floor, log10

def num_digits(n):
    return floor(log10(n))+1

def digit_in_pos(n, pos):
    return (n // 10**pos) % 10

def set_digit(n, pos, new):
    "change digit in given position to new value"
    old = digit_in_pos(n, pos)
    return n + (new - old)*10**pos

def is_digitally_delicate_prime(n):
    if not isprime(n):
        return False
    nd = num_digits(n)
    for pos in range(nd):
        for newdigit in range(1, 10):
            if digit_in_pos(n, pos) == newdigit:
            m = set_digit(n, pos, newdigit)
            if isprime(m):
                print("m = ", m)
                return False
    return True

for p in [294001, 505447, 584141, 604171, 971767, 1062599]:

[1] Michael Filaseta and Jeremiah Southwick. Primes that become composite after changing an arbitrary digit. Mathematics of Computation. Article electronically published on November 24, 2020.

Refinements to the prime number theorem

Let π(x) be the number of primes less than x. The simplest form of the prime number theorem says that π(x) is asymptotically equal to x/log(x), where log means natural logarithm. That is,

\lim_{x\to\infty} \frac{\pi(x)}{x/\log(x)} = 1

This means that in the limit as x goes to infinity, the relative error in approximating π(x) with x/log(x) goes to 0. However, there is room for improvement. The relative approximation error goes to 0 faster if we replace x/log(x) with li(x) where

\text{li}(x) = \int_0^x \frac{dt}{\log t}

The prime number theorem says that for large x, the error in approximating π(x) by li(x) is small relative to π(x) itself. It would appear that li(x) is not only an approximation for π(x), but it is also an upper bound. That is, it seems that li(x) > π(x). However, that’s not true for all x.

Littlewood proved in 1914 that there is some x for which π(x) > li(x). We still don’t know a specific number x for which this holds, though we know such numbers exist. The smallest such x is the definition of Skewes’ number. The number of digits in Skewes’ number is known to be between 20 and 317, and is believed to be close to the latter.

Littlewood not only proved that li(x) – π(x) is sometimes negative, he proved that it changes sign infinitely often. So naturally there is interest in estimating li(x) – π(x) for very large values of x.

A new result was published a few days ago [1] refining previous bounds to prove that

|\pi(x) - \text{li}(x)| \leq 235 x (\log x)^{0.52} \exp(-0.8\sqrt{\log x})

for all x > exp(2000).

When x = exp(2000), the right side is roughly 10857 and π(x) is roughly 10865, and so the relative error is roughly 10-8. That is, the li(x) approximation to π(x) is accurate to 8 significant figures, and the accuracy increases as x gets larger.


[1] Platt and Trudgian. The error term in the prime number theorem. Mathematics of Computation. November 16, 2020.

The smallest number with a given number of divisors

Suppose you want to find the smallest number with 5 divisors. After thinking about it a little you might come up with 16, because

16 = 24

and the divisors of 16 are 2k where k = 0, 1, 2, 3, or 4.

This approach generalizes: For any prime q, the smallest number with q divisors is 2q-1.

Now suppose you want to find the smallest number with 6 divisors. One candidate would be 32 = 25, but you could do better. Instead of just looking at numbers divisible by the smallest prime, you could consider numbers that are divisible by the two smallest primes. And in fact

12 = 22 3

is the smallest number with 6 divisors.

This approach also generalizes. If h is the product of 2 primes, say h = pq where pq, then the smallest number with h divisors is

2p-1 3q-1.

The divisors come from letting the exponent on 2 range from 0 to p-1 and letting the exponent on 3 range from 0 to q-1.

For example, the smallest number with 35 divisors is

5184 = 27-1 35-1.

Note that we did not require p and q to be different. We said pq, and not p > q. And so, for example, the smallest number with 25 divisors is

1296 = 25-1 35-1.

Now, suppose we want to find the smallest number with 1001 divisors. The number 1001 factors as 7*11*13, which has some interesting consequences. It turns out that the smallest number with 1001 divisors is

213-1 311-1 57-1.

Does this solution generalize? Usually, but not always.

Let h = pqr where p, q, and r are primes with pqr. Then the smallest number with h divisors is

2p-1 3q-1 5r-1

with one exception. The smallest number with 8 divisors would be 30 = 2*3*5 if the theorem always held, but in fact the smallest number with 8 divisors is 24.

In [1] M. E. Gorst examines the exceptions to the general pattern. We’ve looked at the smallest number with h divisors when h is the product of 1, or 2, or 3 (not necessarily distinct) primes. Gorst considers values of h equal to the product of up to 6 primes.

We’ve said that the pattern above holds for all h the product of 1 or 2 primes, and for all but one value of h the product of 3 primes. There are two exceptions for h the product of 4 primes. That is, if h = pqrs where pqrs are primes, then the smallest number with h divisors is

2p-1 3q-1 5r-1 7s-1

with two exceptions. The smallest number with 24 divisors is 23 × 3 × 5, and the smallest number with 3 × 23 divisors is 23 × 32 × 5.

When h is the product of 5 or 6 primes, there are infinitely many exceptions, but they have a particular form given in [1].

The result discussed here came up recently in something I was working on, but I don’t remember now what. If memory serves, which it may not, I wanted to assume something like what is presented here but wasn’t sure it was true.


[1] M. E. Grost. The Smallest Number with a Given Number of Divisors. The American Mathematical Monthly, September 1968, pp. 725-729.

Test for divisibility by 13

There are simple rules for telling whether a number is divisible by 2, 3, 4, 5, and 6.

  • A number is divisible by 2 if its last digit is divisible by 2.
  • A number is divisible by 3 if the sum of its digits is divisible by 3.
  • A number is divisible by 4 if the number formed by its last two digits is divisible by 4.
  • A number is divisible by 5 if its last digit is divisible by 5.
  • A number is divisible by 6 if it is divisible by 2 and by 3.

There is a rule for divisibility by 7, but it’s a little wonky. Let’s keep going.

  • A number is divisible by 8 if the number formed by its last three digits is divisible by 8.
  • A number is divisible by 9 if the sum of its digits is divisible by 9.
  • A number is divisible by 10 if its last digit is 0.

There’s a rule for divisibility by 11. It’s a little complicated, though not as complicated as the rule for 7. I describe the rule for 11 in the penultimate paragraph here.

A number is divisible by 12 if it’s divisible by 3 and 4. (It matters here that 3 and 4 are relatively prime. It’s not true, for example, that a number is divisible by 12 if it’s divisible by 2 and 6.)

But what do you do when you get to 13?

Testing divisibility by 7, 11, and 13

We’re going to kill three birds with one stone by presenting a rule for testing divisibility by 13 that also gives new rules for testing divisibility by 7 and 11. So if you’re trying to factor a number by hand, this will give a way to test three primes at once.

To test divisibility by 7, 11, and 13, write your number with digits grouped into threes as usual. For example,


Then think of each group as a separate number — e.g. 11, 37, and 989 — and take the alternating sum, starting with a + sign on the last term.

989 – 37 + 11

The original number is divisible by 7 (or 11 or 13) if this alternating sum is divisible by 7 (or 11 or 13 respectively).

The alternating sum in our example is 963, which is clearly 9*107, and not divisible by 7, 11, or 13. Therefore 11,037,989 is not divisible by 7, 11, or 13.

Here’s another example. Let’s start with


The alternating sum is

518 – 498 + 894 – 4 = 910

The sum takes a bit of work, but less work than dividing a 10-digit number by 7, 11, and 13.

The sum 910 factors into 7*13*10, and so it is divisible by 7 and by 13, but not by 11. That tells us 4,894,498,518 is divisible by 7 and 13 but not by 11.

Why this works

The heart of the method is that 7*11*13 = 1001. If I subtract a multiple of 1001 from a number, I don’t change its divisibility by 7, 11, or 13. More than that, I don’t change its remainder by 7, 11, or 13.

The steps in the method amount to adding or subtracting multiples of 1001 and dividing by 1000. The former doesn’t change the remainder by 7, 11, or 13, but the latter multiplies the remainder by -1, hence the alternating sum. (1000 is congruent to -1 mod 7, mod 11, and mod 13.) See more formal argument in footnote [1].

So not only can we test for divisibility by 7, 11, and 13 with this method, we can also find the remainders by 7, 11, and 13. The original number and the alternating sum are congruent mod 1001, so they are congruent mod 7, mod 11, and mod 13.

In our first example, n = 11,037,989 and the alternating sum was m = 963. The remainder when m is divided by 7 is 4, so the remainder when n is divided by 7 is also 4. That is, m is congruent to 4 mod 7, and so n is congruent to 4 mod 7. Similarly, m is congruent to 6 mod 11, and so n is congruent to 6 mod 11. And finally m is congruent to 1 mod 13, so n is congruent to 1 mod 13.

Related posts

[1] The key calculation is
10^3 \equiv -1 \bmod 1001 \implies \sum_{i=0}^n a_i10^{3i} \equiv \sum_{i=0}^n (-1)^ia_i \bmod 1001