Number theory

Tighter bounds in the prime number theorem

Posted on by John

The most elementary form of the prime number theorem says that π(x), the number of prime numbers less than x, is asymptotically equal to x / log(x). That’s true, but a more accurate result says π(x) is asymptotically equal to li(x) where

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

Five years ago I wrote about a result that was new at the time, giving a bound on |π(x) − li(x)| for x > exp(2000). This morning I saw a result in a blog post by Terence Tao that says

\left| \pi(x) - \text{li}(x) \right| \leq 9.2211\, x\sqrt{\log(x)} \exp\left( -0.8476 \sqrt{\log(x)} \right)

for all x ≥ 2. The result comes from this paper.

The new bound has the same form as the bound from five years ago but with smaller constants.

Efficiently computing multiple modular inverses at once

Posted on by John

Suppose you have a large prime number M and you need to find the inverse of several numbers mod M.  Montgomery’s trick is a way to combine the computation of the inverses to take less time than computing the inverses individually. Peter Montgomery (1947–2020) came up with this trick in 1985.

We will illustrate Montgomery’s trick by inverting three numbers—a, b, and c—though the trick extends to any number of numbers. It is commonly used in cryptography.

Modular inverses are much slower to calculate than modular products, so doing fewer of the former and more of the latter is a good tradeoff. Montgomery’s method only calculates one modular inverse, regardless of how many numbers need to be inverted.

The idea is to directly invert the product of all the numbers and use multiplication to find the inverses of the individual numbers. In our case, we compute

x = ab
y = cy = abc
x−1 = cy−1
b−1 = ax−1
a−1 = bx−1

To show that this actually saves time, we’ll run some Python code to invert three random numbers modulo a very large prime, much larger than occurs in practice. The reason is to make the computation time longer and easier to demonstrate. In practice, Montgomery’s trick saves a little time off of a lot of calculations. Here we’ll save a lot of time off a handful of calculations.

import sys
import time
from secrets import randbelow

# extend the default maximum integer size
sys.set_int_max_str_digits(100000)

# the 32nd Mersenne prime
M = 2**756839 - 1

def simple(a, b, c, M):
    return [pow(x, -1, M) for x in [a, b, c]]

def montgomery(a, b, c, M):
    x = a*b % M
    y = x*c % M
    yinv = pow(y, -1, M)
    cinv = x*yinv % M
    xinv = c*yinv % M
    binv = a*xinv % M
    ainv = b*xinv % M
    return [ainv, binv, cinv]
    
a = randbelow(M)
b = randbelow(M)
c = randbelow(M)

start = time.perf_counter()
result = simple(a, b, c, M)
elapsed = time.perf_counter() - start
print(elapsed)

start = time.perf_counter()
result = montgomery(a, b, c, M)
elapsed = time.perf_counter() - start
print(elapsed)

When we ran this, the direct approach took 121.8 seconds, and Montgomery’s trick took 47.6 seconds.

Related posts

Primecoin primality test

Posted on by John

When I wrote about how Primecoin uses prime chains for proof of work, I left out a detail.

To mine a new Primecoin block, you have to find a prime chain of specified length that starts with a number that is a multiple of the block header hash. According to the Primecoin whitepaper

Another important property of proof-of-work for cryptocurrency is non-reusability. … To achieve this, the prime chain is linked to the block header hash by requiring that its origin be divisible by the block header hash.

So given a hash h, you have to find k such that kh is the origin of a prime chain, where “origin” is defined in footnote [2] here.

Strictly speaking, the primes in a Primecoin prime chain are probable primes. Someone verifying a Primecoin blockchain will be satisfied that the block is authentic if the necessary numbers are prime according to the test used by the Primecoin software. Whether the numbers are actually prime is irrelevant.

Probabilistic primality tests are much more efficient than deterministic tests, and most applications requiring primes, such as RSA encryption, actually use probable primes. If a number is rejected as a probable prime, it’s certainly not a prime. If it is accepted as a prime, it very like is prime. More on probable primes here.

If you write your own Primecoin mining software using a different primality test, that’s fine as long as you actually find primes. But if one time you find pseudoprimes, composite numbers that pass your primality test, then your block will be rejected unless your pseudoprimes also pass the Primecoin software’s primality test.

Primecoin’s primality test

So what is Primecoin’s (probable) primality test? We quote the Primecoin whitepaper:

The classical Fermat test of base 2 is used together with Euler-Lagrange-Lifchitz test to verify probable primality for the prime chains. Note we do not require strict primality proof during verification, as it would unnecessarily burden the efficiency of verification. Composite number that passes Fermat test is commonly known as pseudoprime. Since it is known by the works of Erdös and Pomerance that pseudoprimes of base 2 are much more rare than primes, it suffices to only verify probable primality.

Fermat’s test

The Fermat test with base b says to test whether a candidate number n satisfies

bn − 1 = 1 mod n.

If n is prime, the equation above will hold. The converse is usually true: if the equation above hold, n is likely to be prime. There are exceptions, such as b = 2 and n = 561 = 3 × 11 × 17.

Euler-Lagrange-Lifchitz test

But what is the Euler-Lagrange-Lifchitz test? The whitepaper links here, a page by Henri Lifchitz that gives results generalizing those of Leonard Euler and Joseph-Louis Lagrange.

Testing 2p + 1

Suppose p is an odd prime. Then the Euler-Lagrange test says that if p = 3 mod 4, then q = 2p + 1 is also prime if and only if

2p = 1 mod q.

Lifchitz gives three variations on this result. First, if p = 1 mod 4, then q = 2p + 1 is also prime if and only if

2p = −1 mod q.

Together these two theorems give a way to test with certainty whether 2p + 1, i.e. the next term in a Cunningham chain of the first kind, is prime. However, the test is still probabilistic if we don’t know for certain that p is prime.

Testing 2p − 1

The last two results of Lifchitz are as follows. If p and q = 2p − 1 are prime, then

2p − 1 = 1 mod q

if p = 1 mod 4 and

2p − 1 = −1 mod q

if p = 3 mod 4.

The converses of these two statements give probable prime tests for q if we know p is prime.

So we can verify a (probable) prime chain by using Fermat’s test to verify that the first element is (probably) prime and use the Euler-Lagrange-Lifchitz test that the rest of the numbers in the chain are (probably) prime.

Related posts

Bi-twin prime chains

Posted on by John

I mentioned bi-twin prime chains in the previous post, but didn’t say much about them so as not to interrupt the flow of the article.

A pair of prime numbers are called twins if they differ by 2. For example, 17 and 19 are twin primes.

A bi-twin chain is a sequence of twin primes in which the average of each twin pair is double that of the preceding pair. For example,

5, 7, 11, 13

is a bi-twin chain because (5, 7) and (11, 13) are twin pairs of primes, and the average of the latter pair is twice the average of the former pair.

There is some variation in how to describe how long a bi-twin chain is. We will define the length of a bi-twin chain to be the number prime pairs it contains, and so the chain above has length 2. The BiTwin Record page describes a bi-twin chain of length k as a chain with k − 1 links.

It is widely believed that there are infinitely many twin primes, but this has not been proven. So we don’t know for certain whether there are infinitely many bi-twin prime chains of length 1, much less chains of longer length.

The bi-twin chain of length three with smallest beginning number is

211049, 211051, 422099, 422101, 844199, 844201

Bi-twin records are expressed in terms of primorial numbers. I first mentioned primorial numbers a few days ago and now they come up again. Just as n factorial is the product of the positive integers up to n, n primorial is the product of the primes up to n. The primorial of n is written n#. [1]

The longest known bi-twin chains start with

112511682470782472978099, 112511682470782472978101

and has length 9.

We can verify the chains mentioned above with the following Python code.

def bitwin_length(average):
    n = average
    c = 0
    while isprime(n-1) and isprime(n+1):
        c += 1
        n *= 2
    return c

for n in [6, 211050, 112511682470782472978100]:
    print(bitwin_length(n))

[1] There are two conventions for defining n#. The definition used here, and on the BiTwin Record page, is the product of primes up to n. Another convention is to define n# to be the product of the first n primes.

The priomorial function in Sympy takes two arguments and will compute either definition, depending on whether the second argument is True or False. The default second argument is True, in which case primorial(n) returns the product of the first n primes.

Prime chains

Posted on by John

The title of this post has a double meaning. We will look at chains in the sense of number theory and in the sense of cryptocurrency, i.e. Cunningham chains and blockchains, that involve prime numbers.

Cunningham chains

A chain of primes is a sequence of prime numbers in which each is almost double its predecessor. That is, the next number after p is 2p ± 1.

In a Cunningham chain of the first kind, the successor of p is 2p + 1. For example,

41, 83, 167.

In a Cunningham chain of the second kind, the successor of p is 2p − 1. For example,

19, 37, 73.

Two questions come up immediately. First, are there infinitely many Cunningham chains? Second, how long can a Cunningham chain be? What is known and what is conjectured are at opposite ends of the spectrum. It is unknown whether there are infinitely many Cunningham chains of length 2, but it is conjectured that there are infinitely many Cunningham chains of all lengths.

According to this page, the longest known Cunningham chains of the first kind has length 17, and the longest known Cunningham chain of the second kind has length 19. We can verify these results with the following Python code.

from sympy import isprime

def chain_length(start, kind):
    p = start
    c = 0
    while isprime(p):
        c += 1
        p = 2*p + kind
    return c

print(chain_length(2759832934171386593519, 1))
print(chain_length(79910197721667870187016101, -1))

Bi-twin chains

A number n is the basis of a bi-twin chain of length k if n − 1 is the start of a Cunningham chain of the first kind of length k and n + 1 is the start of a Cunningham chain of the second kind of length k.

I say more about bi-twin prime chains in the next post.

Primecoin

Primecoin was one of the first cryptocurrencies, coming out four years after Bitcoin. Primecoin is still going, though its market cap is six orders magnitude smaller than that of Bitcoin.

What’s interesting about Primecoin is that it uses finding prime chains as its proof of work task [1]. To mint a new Primecoin block, you must find a prime chain of the required length whose origin a multiple of the hash of the block header [2].

Primecoin allows any of the three kinds of prime chains mentioned above: Cunningham chains of the first or second kind, or a bi-twin prime chain. Primecoin adjusts its mining difficulty over time by varying the length of Cunningham or bi-twin chain needed to mint a new block.

There are also cryptocurrencies based on finding prime clusters and prime gaps.

Related posts

[1] Strictly speaking, Primecoin requires finding probable prime chains, as explained here.

[2] The origin of a prime chain is n if the first item in a Cunningham chain of the first kind is n + 1, or if the first item in a Cunningham chain of the first kind or a bi-twin chain is n − 1.

Largest known compositorial prime

Posted on by John

I ran across a blog post here that said a new record has been set for the largest compositorial prime. [1]

OK, so what is a compositorial prime? It is a prime number of the form

n! / n# + 1

where n# denotes n primorial, the product of the prime numbers no greater than n.

The newly discovered prime is

N = 751882!/751882# + 1

It was described in the article cited above as

751882!/751879# + 1,

but the two numbers are the same because there are no primes greater than 751879 and less than 751882, i.e.

751879# = 751882#.

About how large is N? We can calculate the log of the numerator easily enough:

>>> import scipy.special
>>> scipy.special.loggamma(751883)
9421340.780760147

However, the denominator is harder to compute. According to OIES we have

n# = exp((1 + o(1)) n)

which would give us the estimate

log(751882#) ≈ 751882.

So

log10(N) = log(N) / log(10) ≈ (9421340 − 751882) / log(10) ≈ 3765097.

According to this page,

log10(N) = 3765620.3395779

and so our approximation above was good to four figures.

So N has between 3 and 4 million digits, making it much smaller than the largest known prime, which has roughly 41 million digits. Overall, N is the 110th largest known prime.

 

[1] I misread the post at first and thought it said there was a new prime record (skipping over the “compositorial” part) and was surprised because the number is not a Mersenne number. For a long time now the largest known prime has been a Mersenne prime because there is a special algorithm for testing whether Mersenne numbers are prime, one that is much more efficient than testing numbers in general.

 

log2(3) and log2(5)

Posted on by John

AlmostSure on X pointed out that

log2 3 ≈ 19/12,

an approximation that’s pretty good relative to the size of the denominator. To get an approximation that’s as accurate or better requires a larger denominator for log2 5.

log2 5 ≈ 65/28

This above observations are correct, but are they indicative of a more general pattern? Is log2 3 easier to approximate than log2 5 using rational numbers? There are theoretical ways to quantify this—irrationality measures—but they’re hard to compute.

If you look at the series of approximations for both numbers, based on continued fraction convergents, the nth convergent for log2 5 is more accurate than the nth convergent for log2 3, at least for the first 16 terms. After that I ran out of floating point precision and wasn’t sufficiently interested to resort to extended precision.

Admittedly this is a non-standard way to evaluate approximation error. Typically you look at the approximation error relative to the size of the denominator, not relative to the index of the convergents.

Here’s a more conventional comparison, plotting the log of approximation error against the log of the denominators.

Continued fraction posts

Multiples with no large digits

Posted on by John

Here’s a curious theorem I stumbled across recently [1]. Take an integer N which is not a multiple of 10. Then there is some multiple of N which only contains the digits 1, 2, 3, 4, and 5.

For example, my business phone number 8324228646 has a couple 8s and a couple 6s. But

6312 × 8324228646 = 52542531213552

which contains only digits 1 through 5.

For a general base b, let p be the smallest prime factor of b. Then for every integer N that is not a multiple of b, there is some multiple of N whose base b representation contains only the digits 1, 2, 3, …, b/p.

This means that for every number N that is not a multiple of 16, there is some k such that the hex representation of kN contains only the digits 1 through 8. For example, if we take the magic number at the beginning of every Java class file, 0xCAFEBABE, we find

1341 × CAFEBABEhex = 42758583546hex.

In the examples above, we’re looking for multiple containing only half the possible digits. If the largest prime dividing the base is larger than 2 then we can find a multiples with digits in a smaller range. For example, in base 35 we can find a multiple containing only the digits 1 through 7.

[1] Gregory Galperin and Michael Reid. Multiples Without Large Digits. The American Mathematical Monthly, Vol. 126, No. 10 (December 2019), pp. 950-951.

Zero knowledge proof of compositeness

Posted on by John

A zero knowledge proof (ZKP) answers a question without revealing anything more than answer. For example, a digital signature proves your possession of a private key without revealing that key.

Here’s another example, one that’s more concrete than a digital signature. Suppose you have a deck of 52 cards, 13 of each of spades, hearts, diamonds, and clubs. If I draw a spade from the deck, I can prove that I drew a spade without showing which card I drew. If I show you that all the hearts, diamonds, and clubs are still in the deck, then you know that the missing card must be a spade.

Composite numbers

You can think of Fermat’s primality test as a zero knowledge proof. For example, I can convince you that the following number is composite without telling you what its factors are.

n = 244948974278317817239218684105179099697841253232749877148554952030873515325678914498692765804485233435199358326742674280590888061039570247306980857239550402418179621896817000856571932268313970451989041

Fermat’s little theorem says that if n is a prime and b is not a multiple of n, then

bn−1 = 1 (mod n).

A number b such that bn−1 ≠ 1 (mod n) is a proof that n is not prime, i.e. n is composite. So, for example, b = 2 is a proof that n above is composite. This can be verified very quickly using Python:

    >>> pow(2, n-1, n)
    10282 ... 4299

I tried the smallest possible base [1] and it worked. In general you may have to try a few bases. And for a few rare numbers (Carmichael numbers) you won’t be able to find a base. But if you do find a base b such that bn−1 is not congruent to 1 mod n, you know with certainty that n is composite.

Prime numbers

The converse of Fermat’s little theorem is false. It can be used to prove a number is not prime, but it cannot prove that a number is prime. But it can be used to show that a number is probably prime. (There’s some subtlety as to what it means for a number to probably be prime. See here.)

Fermat’s little theorem can give you a zero knowledge proof that a number is composite. Can it give you a zero knowledge proof that a number is prime? There are a couple oddities in this question.

First, what would it mean to have a zero knowledge proof that a number is prime? What knowledge are you keeping secret? When you prove that a number is composite, the prime factors are secret (or even unknown), but what’s the secret when you say a number is prime? Strictly speaking a ZKP doesn’t have to keep anything secret, but in practice it always does.

Second, what about the probability of error? Zero knowledge proofs do not have to be infallible. A ZKP can have some negligible probability of error, and usually do.

It’s not part of the definition, but practical ZKPs must be easier to verify than the direct approach to what they prove. So you could have something like a primality certificate that takes far less computation to verify than the computation needed to determine from scratch that a number is prime.

Proving other things

You could think of non-constructive proofs as ZKPs. For example, you could think of the intermediate value theorem as a ZKP: it proves that a function has a zero in an interval without giving you any information about where that zero may be located.

What makes ZKPs interesting in application is that they can prove things of more general interest than mathematical statements [2]. For example, cryptocurrencies can provide ZKPs that accounting constraints hold without revealing the inputs or outputs of the transaction. You could prove that nobody tried to spend a negative amount and that the sum of the inputs equals the sum of the outputs.

Related posts

[1] You could try b = 1, but then bn−1 is always 1. This example shows that the existence of a base for which bn−1 = 1 (mod n) doesn’t prove anything.

[2] You might object that accounting rules are mathematical statements, and of course they are. But they’re of little interest to mathematicians and of great interest to the parties in a transaction.

Closest consecutive reciprocal sum to an integer

Posted on by John

József Kürschák proved in 1908 that the function

f(m, n) = \sum_{k = m}^n \frac{1}{k}

is never an integer for 0 < m < n. In particular, the harmonic numbers

H_n = f(1, n)

are never integers for n > 1.

The function f(m, n) can get arbitrarily close to any integer value by taking m and n large enough, but it can never exactly equal an integer.

For this post, I’d like to look at how close  f(m, n) comes to an integer value when 0 < m < n ≤ N for some large N, say N = 100,000.

Computation strategy

The most naive way to approach this would be to compute f(m, n) for all m and n and keep track of which values were closest to integers. This would be wasteful since it would recompute the same terms over and over. Instead, we could take advantage of the fact that

f(m, n + 1) = f(m, n) + \frac{1}{n+1}

Instead of working with f(m, n), it will be convenient to work with just its fractional part

g(m, n) = f(m, n) - \lfloor f(m, n) \rfloor

because it won’t hurt to throw away the integer parts as we go. The values of m and n minimizing g(m, n) will be the values for which f(m, n) comes closest to an integer from above. The values of m and n maximizing g(m, n) will be the values for which f(m, n) comes closest to an integer from below.

We could calculate a matrix with all values of g(m, n), but this would take O(N²) memory. Instead, for each n we will calculate g(m, n), save the maximum and minimum values, then overwrite that memory with g(m, n + 1). This approach will only take O(N) memory.

Floating point error

When we compute f(m, n) for large values of n, can we rely on floating point arithmetic?

If N = 100,000, f(m, n) < 16 = 24. A floating point fraction has 53 bits, so we’d expect each addition to be correct to within an error of 2−49 and so we’d expect our total error to be less than 2−49 N.

Python code

The following code computes the values of g(m, n) closest to 0 and 1.

import numpy as np

N = 100_000
f_m = np.zeros(N+1) # working memory

# best values of m for each n
min_fm = np.zeros(N+1)
max_fm = np.zeros(N+1)

n = 2
f_m[1] = 1.5

for n in range(3, N+1):
    f_m[n-1] = 1/(n-1)
    f_m[1:n] += 1/n
    f_m[1:n] -= np.floor(f_m[1:n])
    min_fm[n] = np.min(f_m[1:n])
    max_fm[n] = np.max(f_m[1:n])    

print(min(min_fm[3:]))
print(max(max_fm))

This reports a minimum value of 5.2841953035454026e-11 and a maximum value of 0.9999999996613634. The minimum value is closer to 0 than our (pessimistic) error estimate, though the maximum value is further from 1 than our error estimate.

Modifying the code a bit shows that the minimum occurs at (27134, 73756), and this the input to g that is within our error estimate. So we can be confident that it is the minimum, though we can’t be confident of its value. So next we turn to Mathematica to find the exact value of g(27133, 73756) as a rational number, a fraction with 32024 digits in the numerator and denominator, and convert it to a floating point number. The result agrees with our estimate in magnitude and to four significant figures.

So in summary

\sum_{k = 27134}^{73756} \frac{1}{k} \approx 1

with an error on the order of 10−11, and this is the closest value of f(m, n) to an integer, for 0 < m < n ≤ 100,000.

Related posts