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

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}

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}

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} 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()

• This formula seems familiar, Did you previously ask another question about it? Commented Aug 31, 2023 at 14:35
• @StevenClark yes, about its correctness at MO, and indeed you validated it ;) I am intending to do something with it, if possible! Commented Aug 31, 2023 at 15:08
• Using formula (1) in my answer below (which I've verified for $0\le x\le 27,000$), I believe formula (1) in your question above should be equivalent to $$f(x)=-\sum\limits_{k\le x\land Gpf(k)\leq \frac{x}{k}} \mu(k)\frac{x/k}{\log(x/k)},$$ but this formula oscillates around the $x$-axis instead of $\sqrt{x}$ so something seems to be wrong somewhere. Commented Sep 12, 2023 at 18:14
• @StevenClark I have reviewed the Python program I have used for the calculations and I find no error... have you checked yours? If you want I can post my code as an edit for you to check Commented Sep 13, 2023 at 16:28
• I believe in calculate_expression_3, "if pij < n / pij:" should be "if primes[j] <= n / pij:". There's a similar problem in calculate_expression_4 and calculate_expression_5. Commented Sep 14, 2023 at 16:08

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.

• thanks! Even not being conclusive, your comments are really useful. +1! Commented Aug 31, 2023 at 20:35
• I have just realized your improvement on the code, great job! Commented Sep 18, 2023 at 13:07