Natural density of the set of simple numbers

Question

Let us call $n>1$ simple if every prime power $q$ with $q-1 \mid n-1$ is a prime number. (Please let me know if there is already an established name for these numbers.) The simple numbers $\leq 100$ are

2,3,5,6,11,12,14,18,20,21,23,24,26,30,35,38,39,42,44,45,47,48,51,54,56,59,60,62,66,68,69,72,74,75,77,80,83,84,86,87,90,93,95,96,98

This is A366343. It seems that on average about every second number is simple. This leads to the following question:

Question. Does the set of simple numbers have a natural density? If yes, is there an explicit formula? If not, can we at least compute lower asymptotic density and the upper asymptotic density? These are: $$\liminf_{n \to \infty} \frac{\# \text{ simple numbers } \leq n}{n}, \quad \limsup_{n \to \infty} \frac{\# \text{ simple numbers } \leq n}{n}$$

In numerical experiments with SageMath, the value jumps a bit between values around 0.46211. Here is my code:

def n_powers(n):
    """Computes the prime powers q with q-1 | n-1"""
    return [x+1 for x in divisors(n-1) if (x+1).is_prime_power()]

def is_simple(n):
    """Checks if n is a simple number"""
    return n > 1 and all(q.is_prime() for q in n_powers(n))

def simple_density(limit):
   """prints the proportions of simple numbers among the numbers 2 <= n <= limit"""
   n = 0
   simple_count = 0
   while (n < limit):
      n+=1
      if (is_simple(n)):
         simple_count += 1
      print(float(simple_count / n))

Then, simple_density(1000000) prints a lot of approximations.

Background. For simple numbers, there is a simple (!) proof of Jacobson's Theorem. See Equational proofs of Jacobson's Theorem.

Answer 1

We can get a quite good upper bound, and also a reasonably close lower bound, by determining and excluding sufficient values of $n$ which are not simple, i.e., where there's a prime power $q$ with $q - 1 \mid n - 1$ with $q$ not being prime. In particular, we'll use the smaller resulting congruence values to significantly limit what later values are required to be checked and used for increasing the exclusion. First, with the $3$ smallest such $q$, i.e., $4$, $8$ and $9$, we get

$$4 - 1 \mid n - 1 \;\;\to\;\; n \equiv 1 \pmod{3}$$

$$8 - 1 \mid n - 1 \;\;\to\;\; n \equiv 1 \pmod{7}$$

$$9 - 1 \mid n - 1 \;\;\to\;\; n \equiv 1 \pmod{8}$$

Note we can exclude checking any larger prime power where the resulting modulo, i.e., $q - 1$, is an integral multiple of any of the previous modulos, starting with the $3$ above. Thus, since the squares of all odd integers are congruent to $1$ modulo $8$, we only need to check odd powers of odd primes. Also, since $p^2 \equiv 1 \pmod{3}$ for primes $p \gt 3$, we only need to check odd powers of those primes which are congruent to $2$ modulo $3$. Further, for powers of $2$, to avoid multiples of $7$, we can exclude those where the exponent is a multiple of $3$. This already then becomes a quite sparse list, with it becoming very sparse rather quickly when more modulos are considered (e.g., we only need to check where the exponents themselves are primes). The only other values less than $100$, which add new congruences to check, are

$$27 - 1 \mid n - 1 \;\;\to\;\; n \equiv 1 \pmod{26}$$

$$32 - 1 \mid n - 1 \;\;\to\;\; n \equiv 1 \pmod{31}$$

The only two more after this (note that, for example, since $125 - 1 = 4(31)$, it's excluded since $31$ is already being used) less than $1$,$000$ are

$$128 - 1 \mid n - 1 \;\;\to\;\; n \equiv 1 \pmod{127}$$

$$243 - 1 \mid n - 1 \;\;\to\;\; n \equiv 1 \pmod{242}$$

To get the natural density resulting from combining the $q - 1$ values, we need to use the appropriate inclusion-exclusion adjustments, with Gustave C. Pekara's 1972 PhD thesis The Asymptotic Density of Certain Integer Sequences giving the specific formula, and proof of it, as Theorem $4.1$, on its page $52$ (page $56$ in the linked PDF file). In particular, for $m \ge 1$ values, with $q_i - 1$ for $1 \le i \le m$ being the checked values, the natural density $d_m$ is given by

$$\begin{equation}\begin{aligned} d_m = & \sum_{i = 1}^{m}\frac{1}{q_i - 1} - \sum_{i \lt j}\frac{1}{\operatorname{lcm}(q_i - 1, q_j - 1)} + \\ & \sum_{i \lt j \lt k}\frac{1}{\operatorname{lcm}(q_i - 1, q_j - 1, q_k - 1)} - \ldots + \frac{(-1)^{m - 1}}{\operatorname{lcm}(q_1 - 1, q_2 - 1, \ldots, q_m - 1)} \end{aligned}\end{equation}$$

For $m = 3$, the natural density of integers being excluded is

$$d_3 = \frac{1}{3} + \frac{1}{7} + \frac{1}{8} - \frac{1}{3\cdot 7} - \frac{1}{3\cdot 8} - \frac{1}{7\cdot 8} + \frac{1}{3\cdot 7\cdot 8} = 0.5$$

The adjustment for the $m$'th value being added, call it $a_m$ (so then $d_m = d_{m-1} + a_m$), is basically as stated in Computing the density of a set of multiples, i.e., it's

$$\begin{equation}\begin{aligned} a_m = & \frac{1}{q_m - 1} - \sum_{i \lt m}\frac{1}{\operatorname{lcm}(q_i - 1, q_m - 1)} + \\ & \sum_{i \lt j \lt m}\frac{1}{\operatorname{lcm}(q_i - 1, q_j - 1, q_m - 1)} - \ldots + \frac{(-1)^{m - 1}}{\operatorname{lcm}(q_1 - 1, q_2 - 1, \ldots, q_m - 1)} \end{aligned}\end{equation}$$

Thus, the adjustment $a_4$ to include the next value, i.e., $26$, using the $\operatorname{lcm}$ to get the denominator values (as shown above, and also stated in your comment) for this and all larger values, is

$$\begin{equation}\begin{aligned} a_4 & = \frac{1}{26} - \frac{1}{3\cdot 26} - \frac{1}{7\cdot 26} - \frac{1}{8\cdot 13} + \frac{1}{3\cdot 7\cdot 26} + \\ & \;\;\;\;\; \frac{1}{3\cdot 8\cdot 13} + \frac{1}{7\cdot 8\cdot 13} - \frac{1}{3\cdot 7\cdot 8\cdot 13} \\ & \approx 0.0165 \end{aligned}\end{equation}$$

Thus, a rough simple numbers natural density upper bound is

$$1 - (0.5 + 0.0165) = 0.4835$$

We can continue this for each new value to keep getting better upper bounds. Note that your answer shows that, using a SageMath code for all prime powers up to one or five million, the result is about $0.462118$. This fairly closely agrees with the value, stated in the comment, obtained by checking for all simple numbers up to $75$,$000$,$000$.

---

The $a_m$ adjustment expression given earlier is for the natural density of numbers $n - 1$ where $q_m - 1 \mid n - 1$, but $q_i - 1 \nmid n - 1 \;\forall\; 1 \le i \lt m$. Since we're excluding cases where there's an $1 \le i \lt m$ with $q_i - 1 \mid q_m - 1$ (since then $a_m = 0$), we therefore have

$$0 \lt a_m \lt \frac{1}{q_m - 1}$$

I did several online searches for an algebraic proof of the above expression. Although I didn't find anything, and wasn't able to come up with a proof myself (e.g., due to difficulties dealing with $\operatorname{lcm}$), I suspect this has been proven previously. Regardless, an approach that may work is to equally divide each summation term into the number of non-$m$ index values being used. For example, with the second summation, since the denominator is $\operatorname{lcm}(q_i - 1, q_j - 1, q_m - 1)$, then split each term into two halves. Next, combine all of the terms where each $q_i - 1$ (for $1 \le i \lt m$) is being used, and then find appropriate lower & upper bounds on the sums of those terms.

The fifth and higher terms, as explained, only involve prime powers $q_i$ which are to at least the third power. For powers of $2$, where the exponent is $\ge 5$, we only need to consider odd prime powers. Considering the larger value where these are exponents are just not multiples of $3$, then using the sums of those powers of $2$ forming geometric series, we get that this is a bit less than $\frac{64}{63}\left(\frac{1}{31} + \frac{1}{127}\right) \approx 0.04077$. Thus, using this, a bit more than $d_4$, and a value slightly larger than Apéry's constant (since all the checked higher prime exponents (e.g., $5$, $7$, etc.) of the primes are more than compensated by the unused summation terms), we get for all $m$ that

$$\begin{equation}\begin{aligned} d_m & \lt 0.0408 + 0.5166 + \sum_{i=3}^{\infty}\frac{1}{i^3} \\ & \lt 0.0408 + 0.5166 + \left(1.2021 - 1 - \frac{1}{2^3}\right) \\ & = 0.6345 \end{aligned}\end{equation}$$

Since the $d_m$ form a strictly increasing sequence with an upper bound, the montone convergence theorem states that it converges. The inequality above shows that this limit of $d_{\infty} = \lim_{m\to\infty}d_m$ is less than $1$, so the natural density of simple numbers also exists, and it's the positive $1 - d_{\infty}$.

The $1 - d_m$ values provide upper bounds for the simple numbers natural density, with each new $a_m$ adjustment making it more accurate. After calculating $1 - d_m$ up to some reasonably large $q_m - 1$, we can also then get a fairly good value for the natural density lower bound by determining an upper bound on the remaining $a_m$ adjustments left to add, such as by using the Apéry's constant terms, or perhaps by using something more refined instead.

Regarding analytically determining the simple numbers natural density, the challenges involved with using $\operatorname{lcm}$'s suggests to me that it's unlikely to be successful with my approach here. Instead, a rather different method would be required.

Answer 2

This is a comment to John's answer, but since comments cannot have codeblocks, let me write a CW answer.

John's approach is to let $q_1,q_2,q_3,\dotsc$ be a list of proper prime powers so that $q_i - 1 \not\mid q_n-1$ for all $i < n$. Then, $$1 - \sum_i \frac{1}{(q_i-1)} + \sum_{i < j} \frac{1}{\mathrm{lcm}(q_i-1,q_j-1)} - \sum_{i < j < k} \frac{1}{\mathrm{lcm}(q_i-1,q_j-1,q_k-1)} \pm$$ approximates the natural density. Here is the SageMath code to calculate this sum:

def relevant_powers(limit):
    """lists the proper prime powers q such that q-1 is not a multiple of x-1 for previous powers"""
    powers = []
    for q in srange(2,limit+1):
        if q.is_prime_power() and not q.is_prime() and all([(q-1)%(x-1) > 0 for x in powers]):
            powers.append(q)
    return powers

def approximate_density(limit):
    """approximates the natural density of simple numbers using inclusion-exclusion principle"""
    sum = 0
    powers = relevant_powers(limit)
    for k in range(0,len(powers)+1):
        for selection in Subsets(powers,k):
            sum += (-1)^k * 1/lcm([x-1 for x in selection])
            print(float(sum)) # print progress
    return float(sum)

Then, executing approximate_density(1000000) yields 0.46211885256102014 at some point and it does not change anymore. (Of course, it does, but SageMath's precision is not enough.)

However, executing approximate_density(5000000) yields 0.462118293064357 at some point and it does not change anymore. I am not sure what this means.