Bounding an inclusion-exclusion formula for $M(n)$

Question

Some months ago, I derived the following formula for the Merten's function $M(n)$ using the inclusion-exclusion principle:

$$M(n)=1-\pi\left(n\right)+\sum_{p_{i}\leq\frac{n}{p_i}}\left(\pi\left(\frac{n}{p_{i}}\right)-i\right)-\sum_{p_{i}<p_{j}\leq\frac{n}{p_ip_j}}\left(\pi\left(\frac{n}{p_{i}p_{j}}\right)-j\right)+\sum_{p_{i}<p_{j}<p_{k}\leq\frac{n}{p_ip_jp_k}}\left(\pi\left(\frac{n}{p_{i}p_{j}p_{k}}\right)-k\right)-\dots$$

I suspect that applying bounds for $\pi(n)$ could be useful to bound $M(n)$ in this formulation. For instance, maybe Dussart's bound $\frac {x}{\log (x)-1}<\pi(x)<\frac {x}{\log (x) -1.1}$ could be applicable. However, it still seems a very hard task, and I am not sure if it would be worthy to follow that path.

Question. Is the way proposed (or other similar) worthy to bound $M(n)$ in this formulation, or are there better methods to try?

Thanks!

EDIT

The formula can be divided into two great expressions as $$M(n)=\left(-\pi\left(n\right)+\sum_{p_{i}\leq\frac{n}{p_i}}\left(\pi\left(\frac{n}{p_{i}}\right)\right)-\sum_{p_{i}<p_{j}\leq\frac{n}{p_ip_j}}\left(\pi\left(\frac{n}{p_{i}p_{j}}\right)\right)+\sum_{p_{i}<p_{j}<p_{k}\leq\frac{n}{p_ip_jp_k}}\left(\pi\left(\frac{n}{p_{i}p_{j}p_{k}}\right)\right)-\dots\right)-\left(\sum_{p_{i}\leq\frac{n}{p_i}}i-\sum_{p_{i}<p_{j}\leq\frac{n}{p_ip_j}}j+\sum_{p_{i}<p_{j}<p_{k}\leq\frac{n}{p_ip_jp_k}}k-\dots\right)$$

Applying the Prime Number Theorem to the first of the two great expressions between brackets, yields the expression $$\left(-\frac{n}{\log\left(n\right)}+\sum_{p_{i}<\frac{n}{p_{i}}}\left(\frac{n}{p_{i}\left(\log\left(n\right)-\log\left(p_{i}\right)\right)}\right)-\sum_{p_{i}<p_{j}<\frac{n}{p_{i}p_{j}}}\left(\frac{n}{p_{i}p_{j}\left(\log\left(n\right)-\log\left(p_{i}\right)-\log\left(p_{j}\right)\right)}\right)+\sum_{p_{i}<p_{j}<p_{k}<\frac{n}{p_{i}p_{j}p_{k}}}\left(\frac{n}{p_{i}p_{j}p_{k}\left(\log\left(n\right)-\log\left(p_{i}\right)-\log\left(p_{j}\right)-\log\left(p_{k}\right)\right)}\right)-...\right) (1)$$ I have evaluated sepparately this last expression and the second great expression between brackets for the first thousands of values of $n$, and they oscillate as quasi-symmetric functions over the x-axis, as shown in the following graph:

I attach below the code used for getting the plot showed:

    @author: juanmoreno
"""

import math
import matplotlib.pyplot as plt

def sieve_of_eratosthenes(n):
is_prime = [True] * (n + 1)
is_prime[0] = is_prime[1] = False
primes = []

for p in range(2, int(math.sqrt(n)) + 1):
    if is_prime[p]:
        for i in range(p * p, n + 1, p):
            is_prime[i] = False

for i in range(2, n + 1):
    if is_prime[i]:
        primes.append(i)

return primes

def calculate_expression_1(n):
return (-n / math.log(n))

def calculate_expression_2(n, primes):
result = 0
result_2 = 0
for pi in primes:
    if pi < n / pi:
        result += (n / (pi * (math.log(n) - math.log(pi))))
        result_2 -= primes.index(pi)
        
return result, result_2

def calculate_expression_3(n, primes):
result = 0
result_3 = 0
for i in range(len(primes)):
    for j in range(i + 1, len(primes)):
        pij = primes[i] * primes[j]
        if primes[j] < n / pij:
            result -= (n / (pij * (math.log(n) - 
math.log(primes[i]) - math.log(primes[j]))))
            result_3 += j
            
return result, result_3

def calculate_expression_4(n, primes):
result = 0
result_4 = 0
for i in range(len(primes)):
    for j in range(i + 1, len(primes)):
        pij = primes[i] * primes[j]
        if pij < n / pij:
            for k in range(j + 1, len(primes)):
                pijk = pij * primes[k]
                if primes[k] <= n / pijk:
                    result += (n / (pijk * (math.log(n) - 
 math.log(primes[i]) - math.log(primes[j]) - math.log(primes[k]))))
                    result_4 -= k
                
return result, result_4

def calculate_expression_5(n, primes):
result = 0
result_5 = 0
for i in range(len(primes)):
    for j in range(i + 1, len(primes)):
        pij = primes[i] * primes[j]
        if pij < n / pij:
            for k in range(j + 1, len(primes)):
                pijk = pij * primes[k]
                if pijk < n / pijk:
                    for l in range(k + 1, len(primes)):
                        pijkl = pijk * primes[l]
                        if primes[l] <= n / pijkl:
                            result -= (n / (pijkl * (math.log(n) - math.log(primes[i]) - math.log(primes[j]) - math.log(primes[k])-math.log(primes[l]))))
                            result_5 += l
                    
return result, result_5

# Create lists to store results and square roots of n
results_1 = []
results_2 = []
sqrt_n = []
x_axis = []

for n in range(2, 10001):
    primes = sieve_of_eratosthenes(n)
    result1 = calculate_expression_1(n)
    result2 = calculate_expression_2(n, primes)[0]
    result3 = calculate_expression_3(n, primes)[0]
    result4 = calculate_expression_4(n, primes)[0]
    result5 = calculate_expression_5(n, primes)[0]
    final_result_1 =1 + result1 + result2 + result3 + result4 + result5

    result21 = calculate_expression_2(n, primes)[1]
    result31 = calculate_expression_3(n, primes)[1]
    result41 = calculate_expression_4(n, primes)[1]
    result51 = calculate_expression_5(n, primes)[1]
    final_result_2 = result21 + result31 + result41 + result51


results_1.append(final_result_1)
results_2.append(final_result_2)
sqrt_n.append(math.sqrt(n))
x_axis.append(0)

 # Create a plot to compare results and sqrt(n)
 plt.figure(figsize=(10, 6))
 plt.plot(range(2, 10001), results_1, label='Expression (1)', 
 linewidth=2)
 plt.plot(range(2, 10001), results_2, label='Expression (2)', 
 linewidth=2)
 plt.plot(range(2, 10001), sqrt_n, label='square root of n', 
 linestyle='--', linewidth=2)
 plt.plot(range(2, 10001), x_axis, label='x-axis', linestyle='--', 
 linewidth=2)
 plt.xlabel('n')
 plt.ylabel('Value')
 plt.legend()
 plt.title('Comparison of Final Result and x-axis')
 plt.grid(True)
 plt.show()

Answer 1

This isn't really a conclusive answer, but is too long for a comment.

---

I believe your formula is equivalent to

$$M(x)=-\sum\limits_{k\le x\land Gpf(k)\leq \frac{x}{k}} \mu(k)\left(\pi\left(\frac{x}{k}\right)-\left(\left\{\begin{array}{cc} 1 & k=1 \\ \pi(Gpf(k)) & k\ge 2 \\ \end{array}\right.\right)\right)\tag{1}$$

where $Gpf(k)$ returns the greatest prime factor of $k\ge 2$ or $1$ if $k=1$ (see OEIS entry A006530). I've validated formula (1) above for all integers in the range $0\le x\le 27,000$.

---

I suspect your formula isn't particularly useful for bounding $M(x)$ from $\pi(x)$ as the formula involves the Moebius function $\mu(k)$ and is more complicated than the Mertens function definition

$$M(x)=\sum\limits_{n\le x} \mu(n)\tag{2}$$

or perhaps even the formulas

$$M(x)=\sum\limits_{n\le x} \left(\sum_{d|n} \mu(d)\, \mu\left(\frac{n}{d}\right)\right) \left\lfloor\frac{x}{n}\right\rfloor\tag{3}$$

and

$$M(x)=\sum\limits_{n\le x} \mu(n)\, n\, \Phi\left(\frac{x}{n}\right)\tag{4}$$

where

$$\Phi(x)=\sum\limits_{n\le x} \phi(n)=\sum\limits_{n\le x} \frac{\mu(n)}{n} \left\lfloor\frac{x}{n}\right\rfloor=\sum\limits_{n\le x} n\, M\left(\frac{x}{n}\right)\tag{5}$$

is the totient_summatory_function.

---

You can find quite a bit of information on bounds for $M(x)$ using other methods by searching the Internet.

---

There are several prime-counting functions related to the Mertens function $M(x)$ including the fundamental prime-counting function

$$\pi(x)=\sum\limits_{p\le x} 1=\sum\limits_{n=1}^{log_2(x)} \mu(n)\, K(x^{1/n})=\sum\limits_{n\le x} \nu(n)\, M\left(\frac{x}{n}\right)\tag{6},$$

the prime-power counting function

$$K(x)=\sum\limits_{p^k\le x} 1=\sum\limits_{n=1}^{log_2(x)} \pi(x^{1/n})=-\sum\limits_{n\le x} \mu(n)\, \nu(n)\left\lfloor \frac{x}{n}\right\rfloor=\sum\limits_{n\le x} \Omega(n)\, M\left(\frac{x}{n}\right)\tag{7},$$

and the first and second Chebyshev functions

$$\vartheta(x)=\sum\limits_{p\le x} \log(p)=\sum\limits_{n=1}^{log_2(x)} \mu(n)\, \psi(x^{1/n})=\sum\limits_{n\le x} \log(rad(n))\, M\left(\frac{x}{n}\right)\tag{8}$$

and

$$\psi(x)=\sum\limits_{n\le x} \Lambda(n)=\sum\limits_{n=1}^{log_2(x)} \vartheta(x^{1/n})=-\sum\limits_{n\le x} \mu(n)\, \log(n)\left\lfloor \frac{x}{n}\right\rfloor$$ $$=\sum\limits_{n\le x} \mu(n)\, T\left(\frac{x}{n}\right)=\sum\limits_{n\le x} \log(n)\, M\left(\frac{x}{n}\right)\tag{9}$$

where $\nu(n)$ is the number of distinct primes dividing $n$ (see OEIS entry A001221), $\Omega(n)$ is the number of prime divisors of $n$ counting multiplicities (see OEIS entry A001222), $rad(n)$ is the radical of n also known as the square-free kernel or largest square-free divisor, $\Lambda(n)$ is the von Mangoldt function, and

$$T(x)=\sum\limits_{n\le x} \log(n)=\sum\limits_{n\le x} \psi\left(\frac{x}{n}\right)\tag{10}.$$

---

Note $\lfloor x\rfloor$ and $T(x)$ have very precise asymptotics, and yet none of the related formulas above seem to be useful in proving the Riemann hypothesis.

---

I'm curious whether formulas (7) to (9) above can be inverted in a manner analogous to your inversion of formula (6) above. Noting that $\pi(x)=\nu\left(e^{\theta(x)}\right)=\nu\left(e^{\psi(x)}\right)$, I'm most interested in the inversion of the prime-power counting function $K(x)$ defined in formula (7) above. The $\pi(x)$ function in formula (1) above can be replaced with a sum of $K(x)$ functions forming a nested sum, but I suspect there's a simpler non-nested sum inversion of formula (7) above more analogous to formula (1) above.

---

I believe the Mertens function $M(x)$ can also be evaluated in terms of the prime-counting function $\pi(x)$ and prime-power counting function $K(x)$ as

$$M(x)=\sum\limits_{n=1}^x a(n)\, \left(\pi\left(\frac{x}{n}\right)+1\right)\tag{11}$$

and

$$M(x)=\sum\limits_{n=1}^x b(n)\, \left(K\left(\frac{x}{n}\right)+1\right)\tag{12}$$

where $a(n)$ and $b(n)$ are the Dirichlet inverses of $\nu(n)+1$ and $\Omega(n)+1$ respectively.

---

The following code illustrates a Python implementation of the formula for $M(x)$. The code below is a revision of the Python code posted in the question above modified to calculate $M(x)$ instead of the PNT version of first great expression defined in formula (1) of the question above and also includes a few optimizations. The Python code below is evaluated with $maxn=81$ but has been tested for $1\le n\le 10000$.

#    @author: juanmoreno
# """

import math
import matplotlib.pyplot as plt

def sieve_of_eratosthenes(n):
    is_prime = [True] * (n + 1)
    is_prime[0] = is_prime[1] = False
    primes = []
    for p in range(2, int(math.sqrt(n)) + 1):
        if is_prime[p]:
            for i in range(p * p, n + 1, p):
                is_prime[i] = False

    for i in range(2, n + 1):
        if is_prime[i]:
            primes.append(i)

    return primes

def primePi(x, primes):
    if x < 2:
        n = 0
    elif x >= primes[len(primes) - 1]:
        n = len(primes)
    else:
        n = math.floor(x / math.log(x))
        while primes[n] < x:
            n += 1
        while primes[n] > x:
            n -= 1
        n += 1

    return n

def calculate_expression_1(n, primes):
    return (-primePi(n, primes)+1)

def calculate_expression_2(n, primes):
    result = 0
    i = 0
    pi = primes[i]
    while pi <= n / pi:
        result += primePi(n / pi, primes) - (i + 1)
        i += 1
        pi = primes[i]

    return result

def calculate_expression_3(n, primes):
    result = 0
    i = 0
    pi = primes[i]
    while pi**2 <= n / pi:
        j = i + 1
        pij = pi * primes[j]
        while primes[j] <= n / pij:
            result -= primePi(n / pij, primes) - (j+1)
            j += 1
            pij = pi * primes[j]
        i += 1
        pi = primes[i]

    return result

def calculate_expression_4(n, primes):
    result = 0

    i = 0
    pi = primes[i]
    while pi**3 <= n / pi:
        j = i + 1
        pij = pi * primes[j]
        while primes[j]**2 <= n / pij:
            k = j + 1
            pijk = pij * primes[k]
            while primes[k] <= n / pijk:
                result += primePi(n / pijk, primes) - (k+1)
                k += 1
                pijk = pij * primes[k]
            j += 1
            pij = pi * primes[j]
        i += 1
        pi = primes[i]

    return result

def calculate_expression_5(n, primes):
    result = 0

    i = 0
    pi = primes[i]
    while pi**4 <= n / pi:
        j = i + 1
        pij = pi * primes[j]
        while primes[j]**3 <= n / pij:
            k = j + 1
            pijk = pij * primes[k]
            while primes[k]**2 <= n / pijk:
                l = k + 1
                pijkl = pijk * primes[l]
                while primes[l] <= n / pijkl:
                    result -= primePi(n / pijkl, primes) - (l + 1)
                    l += 1
                    pijkl = pijk * primes[l]
                k += 1
                pijk = pij * primes[k]
            j += 1
            pij = pi * primes[j]
        i += 1
        pi = primes[i]

    return result

# Create list to store results
results = []

maxn = 81
primes = sieve_of_eratosthenes(maxn)
for n in range(1, maxn + 1):
    result1 = calculate_expression_1(n, primes)
    result2 = calculate_expression_2(n, primes)
    result3 = calculate_expression_3(n, primes)
    result4 = calculate_expression_4(n, primes)
    result5 = calculate_expression_5(n, primes)
    final_result = result1 + result2 + result3 + result4 + result5
    results.append(final_result)

print(results)

if True:
    # Create a plot to illustrate results
    plt.figure(figsize=(10, 6))
    plt.plot(range(1, maxn + 1), results, 'bo', label='Final Result', linewidth=2)
    plt.xlabel('n')
    plt.ylabel('Value')
    plt.legend()
    plt.title('Illustration of Final Result')
    plt.grid(True)
    plt.show()

The code above produces the results

$$[1, 0, -1, -1, -2, -1, -2, -2, -2, -1, -2, -2, -3, -2, -1, -1, -2, -2, -3, -3, -2, -1, -2, -2, -2, -1, -1, -1, -2, -3, -4, -4, -3, -2, -1, -1, -2, -1, 0, 0, -1, -2, -3, -3, -3, -2, -3, -3, -3, -3, -2, -2, -3, -3, -2, -2, -1, 0, -1, -1, -2, -1, -1, -1, 0, -1, -2, -2, -1, -2, -3, -3, -4, -3, -3, -3, -2, -3, -4, -4, -4]$$

which are consistent with the values listed for $M(x)$ in OEIS entry A002321.