# Mertens' function in time $O(\sqrt x)$

This MathOverflow question seems to indicate that the state of the art in computing $$M(x)=\sum_{n\le x}\mu(n)$$ takes time $\Theta(n^{2/3}(\log\log n)^{1/3}),$ which matches my understanding. Recently I came across a paper [1] which gives, almost as an afterthought, an algorithm for computing $M(x)$ in $O(\sqrt x)$ time (in section 3.2).

The algorithm itself seems to be correct, being derived in a straightforward way from the identity $$M(x)=1-\sum_{d\ge2}M(x/d).$$ But the claim that is runs in time $O(\sqrt x),$ or even $O(x^{1/2+\varepsilon}),$ seems unusual enough for me to ask for verification.

First, this would be a major breakthrough in computing $M(x),$ enough that I would think it would merit more mention than a substep of another algorithm. At the least, I would expect a mention that previous algorithms were slower.

Second, the algorithm is very simple, so that if the time is as indicated the implied constant should be low and the algorithm should be practical. In any case it is not hard to program.

Third, on coding the algorithm I found it to be very slow. In fact, it was slower for those values tested than sum(n=1, 10^7, moebius(n)) in GP which involves factoring each number up to $10^7.$ The time needed to factor numbers of those sizes is about $\sqrt x/\log x$ on average, so that's a $\Theta(n^{3/2}/\log n)$ algorithm (admittedly, well-optimized) beating a $\Theta(n^{1/2})$ algorithm. The constants would have to be worse by a factor of $5\cdot10^6$ for that to happen, and there's nothing in the algorithm to suggest anything that bad.

Of course I may have miscoded it (though I obtained the correct answer) or there may be reasons why the constant factors would be so high for this algorithm. But in any case this seemed suspicious enough to bring up here that I might be enlightened in any case.

Of course even if the result claimed is not correct this is no mark against the author, as the paper is only a preprint and the claim is peripheral in any case.

[1]: Jakub Pawlewicz, Counting square-free numbers (2011), arXiv:1107.4890.

-
I am totally not an expert here, but at first sight, it looks as though the issue is the phrase in the preprint "having all values $M(x/d)$ for $d\ge 2$". I take this to mean that you're supposed to have a pre-computed list of those values. If you did not do this, but your code was recursive instead then this would end up being much slower because for numbers with many divisors, you would end up computing $M(n)$ many times over. –  Anthony Quas May 2 '12 at 6:54
@Anthony: I think you're right. Would you post this as an answer so I can mark this resolved? –  Charles May 2 '12 at 13:14
cool... I get to add a computational # theory badge! –  Anthony Quas May 2 '12 at 13:39

I am totally not an expert here, but at first sight, it looks as though the issue is the phrase in the preprint "having all values $M(x/d)$ for $d\ge 2$". I take this to mean that you're supposed to have a pre-computed list of those values. If you did not do this, but your code was recursive instead then this would end up being much slower because for numbers with many divisors, you would end up computing $M(n)$ many times over.

-
Just to be clear: the article does not claim to give an $O(\sqrt{x})$ time algorithm for computing $M(x)$, and in fact does not consider this problem at all. The algorithm implied by the above formula is (quasi-)linear. –  Dror Speiser May 2 '12 at 22:29
@Dror: Good, that matches my empirical observations. –  Charles May 5 '12 at 16:41

Well,

Although this does not answer your particular question of whether the paper is right, it seems rather straightforward to obtain the $O(x^{1/2+\epsilon})$ result for calculating the sum $\sum_{n \leq x} \mu(n)$ by using the same idea as in my answer to Fastest Algorithm to Compute the Sum of Primes?. Let as in that answer $\Phi(x)$ be a smooth test function such that $\Phi(x)=1$ for $x<0$ and $\Phi(x)=0$ for $x>1$. The difference here is that we consider a sum $$\sum_{n\leq x} \mu(n) =\sum_{n=1}^\infty \mu(n)\Phi \left( \frac {n-x} {\sqrt x} \right ) - \sum_{x< n < x + \sqrt x } \mu(n) \Phi \left(\frac {n-x} {\sqrt x} \right)$$ over the Möbius function instead of a sum over primes. The first sum can be calculated by the integral $$\sum_{n=1}^\infty \mu(n)\Phi \left( \frac {n-x} {\sqrt x} \right )= \frac 1 {2 \pi i} \int_{c-\infty i}^{c+\infty i} \frac 1 {\zeta(s)} \int_0^\infty \Phi \left(\frac{y-x}{\sqrt x} \right)f(y) y^{s-1} dyds, (c>1)$$ which can be calculated by the Odlyzko-Schönhage algorithm in time $O(x^{1/2+\epsilon})$ in the same way as in that answer. The remaining sum will be over an interval of length $\sqrt x$. Sieving techniques can determine the Möbius function for all these numbers fast and the time it will take to calculate that sum is of the order of $O(\sqrt x \log x)$.

Update: I looked at Lagarias-Odlyzkos 1987 paper "Computing $\pi(x)$ an analytic method". On page 8 at the bottom of the page I quote "The basic ideas underlying analytic $\pi(x)$ algorithms can be used for computing other arithmetical functions, such as $$M(x) =\sum_{n\leq x} \mu(n).$$ Thus they were certainly aware of the $O(x^{1/2+\epsilon})$ time complexity for calculating $M(x)$. What I describe above is indeed a variant of the Lagarias-Odlyzko method.

-
That's neat. I'm looking for an exact value for M(x), though. –  Charles May 2 '12 at 14:48
This method gives an exact value of $M(x)$. Since $M(x)$ is an integer it is sufficient to do the calculations sufficiently precise so that the error is $<1/2$. The time complexity is $O(x^{1/2+\epsilon})$, not to be confused with an error term in $M(x)$. –  Johan Andersson May 2 '12 at 14:57
Charles, Maybe my wording was unfortunate in my answer. I changed estimate to calculate. –  Johan Andersson May 2 '12 at 15:01
It is interesting to note that the article is actually quite interesting in regards to this answer. In the same way as you wrote, we can count the number of squarefree integers up to $x$ in time $x^{1/2+\epsilon}$. The article, on the other hand, does something much more elementary, but tricky, and gets the better time of $x^{2/5+\epsilon}$. Even though this is a very specific problem, it makes me wonder what I find it would naturally make anyone wonder. –  Dror Speiser May 2 '12 at 22:24
@Dror, I have written an answer to your comment as another answer (does not fit as a comment) where I show that by combining the elementary and analytic method even better results can be obtained. –  Johan Andersson May 10 '12 at 18:48
This is really a response to Dror Speiser's comment (but it is too long to give as a comment), and gives an analytic $O(n^{1/3+\epsilon})$ time and $O(n^{1/6+\epsilon})$ space algorithm to count the number of square free numbers less than $x$ (thus improving on the complexity of the algorithm given in the Jakub Pawlewicz paper cited in the question). If we do exactly in the same way as in my other answer for $$\sum_{n\leq x} |\mu(n)|$$ instead of $$\sum_{n\leq x} \mu(n)$$ by using the fact that $$\frac{\zeta(s)}{\zeta(2s)} = \sum_{n=1}^\infty |\mu(n)| n^{-s}$$ it is true that the time complexity will be $O(x^{1/2+\epsilon})$. There is however a way to combine the analytic method (which is essentially the same as in Lagarias-Odlyzkos 1987 paper "Computing $\pi(x)$ an analytic method" http://www.dtc.umn.edu/~odlyzko/doc/arch/analytic.pi.of.x.pdf) with the elementary identity $$S(n)=\sum_{d=1}^{\lfloor \sqrt n\rfloor} \mu(d) \lfloor \frac n {d^2} \rfloor$$ that is Pawleviec's starting point in obtaining his algorithm. Let $\Phi$ be defined as in my other answer, i.e. let $\Phi(x)$ be a smooth test function such that $\Phi(x)=1$ for $x<0$ and $\Phi(x)=0$ for $x>1$. Consider the identity $$\sum_{n\leq x} |\mu(n)| =\sum_{n=1}^\infty |\mu(n)| \Phi \left( \frac {n-x} {x^{2/3}} \right ) - \sum_{x< n < x + x^{2/3} } |\mu(n)| \Phi \left(\frac {n-x} {x^{2/3}} \right)$$
The first sum can be estimated by a complex integral of length $x^{1/3+\epsilon}$ and by the Odlyzko-Schönhage algorithm it can be calculated in $O(x^{1/3+\epsilon})$ time (and $O(x^{1/6+\epsilon})$ space). The remaining sum will have length $O(x^{2/3})$ so might at first glance not be so simple to treat. I claim however that in fact it can be calculated in $O(x^{1/3+\epsilon})$ time. This is where using the elementary identity is handy. Let $\Psi(x)=0$ for $x<0$ and $\Psi(x)=\Phi(x)$ for $x \geq 0$ (this will have a discontinuity at $x=0$) By a similar elementary identity we obtain
$$\sum_{x< n < x + x^{2/3} } |\mu(n)| \Phi \left(\frac {n-x} {x^{2/3}} \right) = \sum_{d=1}^{\lfloor \sqrt {x+x^{2/3}}\rfloor} \mu(d) \sum_{k=1}^\infty \Psi \left(\frac {d^2k-x} {x^{2/3}} \right).$$ For any $d$ the inner sum can be calculated fast (certainly in $O(x^\epsilon)$ time), by either noticing that the sum is empty, just contains one element or the Poisson summation formula. The trick now is to see that there will only be $O(x^{1/3})$ integers $d$ where the inner sum (in $k$) is non empty, namely $d$ must either be $O(x^{1/3})$ or belong to an interval $(\sqrt{ x/k},\sqrt{(x+x^{2/3})/k})$ for $k \leq x^{1/3}$ and there are only $O(x^{1/3})$ such $d$. Now if we can calculate the values of $\mu(d)$ fast (simple sieving seems more difficult in this case) we obtain the desired algorithm. Determining $\mu(d)$ basically comes down to factoring the number. However it is sufficient to determine the factors less than $d^{1/3}$ since if we know that all prime factors are greater than $d^{1/3}$, then either it has one or two prime factors. We can determine if it is a prime fast by the Agrawal-Kayal-Saxena algorithm, then of course $\mu(d)=-1$. If it has two prime factors either $\mu(d)=0$, but then it has to be a square, and it can be easily checked, or $\mu(d)=1$. Now factoring the $O(x^{1/3})$ numbers of size $O(x^{1/2})$ into all prime factors less than $x^{1/6}$ will take $O(x^{1/3+\epsilon})$ time and $O(x^{1/6+\epsilon})$ space by Bernstein's algorithm http://cr.yp.to/papers/sf.pdf. I might possibly write up this answer as a paper proper. If anyone has seen these arguments anywhere else or have references I would be interested.
This follows from the identity $\sum_d \mu(d) \sum_k f(d^2k)=\sum_n |\mu(n)| f(n)$ which follows from the fact that $\sum_{d^2|n} \mu(d)=|\mu(n)|$. –  Johan Andersson May 10 '12 at 21:07
Is it possible to use the Lagarias-Odlyzkos method to evaluate $\sum_{d} \mu(d) \Phi(\frac{d^2k-x}{x^{2/3}})$ in $x^{1/3-\delta}$ time (note the $d^2$ here)? If so, it seems you might be able to split the very last sum up and use this for large values of $d$ (where there aren't too many $k$) to get an improvement. –  Mark Lewko May 11 '12 at 3:40