I wonder if there is any efficient way to calculate Möbius function for a array of number 1:1000000



Here are a few papers that might be helpful.

Shallit and Shamir, Number-theoretic functions which are equivalent to number of divisors, Inform. Process. Lett. 20 (1985), no. 3, 151–153, MR0801982 (86k:11076). The review, by Hale Trotter, says $\mu(n)$ can be calculated from a single value $d(n^q)$ where $d$ is the divisor function and $q$ is a prime greater than $1+\log_2n$.

Lioen and van de Lune, Systematic computations on Mertens' conjecture and Dirichlet's divisor problem by vectorized sieving, in From Universal Morphisms to Megabytes: a Baayen Space Odyssey, 421–432, Math. Centrum, Centrum Wisk. Inform., Amsterdam, 1994, MR1490603 (98j:11125). The review, by Jonathan P. Sorenson, says the authors present sieving algorithms for computing $\mu(n)$.

Herman te Riele, Computational sieving applied to some classical number-theoretic problems, in Number Theory in Progress, Vol. 2 (Zakopane-Kościelisko, 1997), 1071–1080, de Gruyter, Berlin, 1999, MR1689561 (2000e:11119). The review, by Marc Deléglise, starts, "This paper is a survey about different applications of sieving in number theory. Finding all the primes belonging to a given interval and computing all the values of the Möbius function on a given interval are obvious examples."


For range 1:1000000 the stupidest possible algorithm (factor each integer, check if it's square free, raise -1 to the right power) will terminate in a blink of an eye. For much larger ranges, an obvious modifican of the sieve of erathsosphenes will be very fast ( just remember to zero every $p$th number at the $p$th step.

  • $\begingroup$ I want to avoid factoring it. I'm also considering Sieve of Eratosthenes $\endgroup$ – piyush_sao Jun 13 '12 at 16:53
  • 21
    $\begingroup$ Calling any algorithm the "stupidest possible" underestimates the power of human ingenuity.... $\endgroup$ – Greg Martin Jun 13 '12 at 17:34
  • 2
    $\begingroup$ @Greg: I double dog dare ya :) $\endgroup$ – Igor Rivin Jun 13 '12 at 17:39
  • 4
    $\begingroup$ For calclulating mu(n) when n>1, pick a prime p at random. If n/p is an integer, use (-1) times mu(n/p), taking care to check if p divides n/p. Otherwise go back and try again. There are stupider variations of course. Gerhard "Dare You Double-Dog Dare Me?" Paseman, 2012.06.13 $\endgroup$ – Gerhard Paseman Jun 13 '12 at 18:03
  • 4
    $\begingroup$ Is erathsosphenes what you get when you cross Eratosthenes with Aristophanes? $\endgroup$ – Gerry Myerson Jun 13 '12 at 23:14

You don't need the bigger machinery of a segmented sieve for such a small range. Here is a simple $O(n\log\log n)$ algorithm to calculate all $\mu(i)$ up to $n$ based on the Sieve of Eratosthenes. In fact, the sieve does fully factor all the square-free numbers; it just doesn't do it one at a time.

Depending on your computer, this approach is practical up to around $2^{30}$ at which point you need to start using higher-precision arithmetic and computing a range of $\mu(i)$ values in batches.

public static int[] GetMu(int max)
    var sqrt = (int)Math.Floor(Math.Sqrt(max));
    var mu = new int[max + 1];
    for (int i = 1; i <= max; i++)
        mu[i] = 1;
    for (int i = 2; i <= sqrt; i++)
        if (mu[i] == 1)
            for (int j = i; j <= max; j += i)
                mu[j] *= -i;
            for (int j = i * i; j <= max; j += i * i)
                mu[j] = 0;
    for (int i = 2; i <= max; i++)
        if (mu[i] == i)
            mu[i] = 1;
        else if (mu[i] == -i)
            mu[i] = -1;
        else if (mu[i] < 0)
            mu[i] = 1;
        else if (mu[i] > 0)
            mu[i] = -1;
    return mu;

Running GetMu(1000000) takes about 10 msec on my computer.

  • $\begingroup$ Doing it the second-most-stupid-way, as described in my answer with Mathematica takes 1sec. I did not try it using Gerhard's extra-stupid way, but it is safe to say that anything works in this range. $\endgroup$ – Igor Rivin Jun 14 '12 at 2:04
  • 1
    $\begingroup$ So, how many different spellings of Eratosthenes can we accumulate in this thread? $\endgroup$ – Gerry Myerson Jun 14 '12 at 5:21
  • $\begingroup$ @Gerry: Thanks. Fixed here ... and in all my other documents. $\endgroup$ – Rick Sladkey Jun 14 '12 at 17:42
  • $\begingroup$ In Rick Sladkey's answer, it seems to me that after running through the sieve, mu[i] should be 0 or +/- i. So why have the extra conditions: else if (mu[i] < 0) mu[i] = 1; else if (mu[i] > 0) mu[i] = -1; ?? $\endgroup$ – Steve Robbins Mar 23 '14 at 4:33
  • $\begingroup$ It's not true that all the prime factors of a number are less than its square root. But, if there is such a factor, there can only be one. $\endgroup$ – Rick Sladkey Mar 23 '14 at 14:55

The key phrase here is for an array. Then try the Sieve of Eratosthenes:

  • initialize an array of length N of 1's
  • for each prime $p < N$ iterate over the array:

    • multiply each $x \equiv 0 \, (\mod p\,)$ by $-1$
    • set each $x \equiv 0 \, (\mod p^2)$ to $\;\,0$

We stop at power $2$ since we are just looking for squarefree (OEIS: A005117).

This is certainly not the fastest but maybe easiest to implement. I downloaded a list of primes off the internet (or you can generate or find your own sieve).

import numpy as np

L = np.ones(100000).astype(int)

for p in P:
    L[::p]    *= -1
    L[::p**2] *=  0 

The output looks good (See also OEIS:A008683)

0,  1,  -1, -1,  0, -1,  1, -1,  0,  0,  1, -1,  0, -1,  1,  
1,  0,  -1,  0, -1,  0,  1,  1, -1,  0,  0,  1,  0,  0, -1

The Sieve

N = np.ones(100000).astype(int)
N[:2] = 0

P = []

p = 1
while p < np.sqrt(100000):

    p = np.argmax(N)
    N[::p] = 0
    P += [p]

P = np.hstack((P, np.where(N > 0)[0]))

This is a sample implementation. The results, they look correct:

547, 557, 563, 569, 571

  • 3
    $\begingroup$ I think you want to not just 'replace 1 by -1' but just 'multiply each entry $\equiv 0\pmod p$ by -1'; your code itself gets it right, but the description suggests that once a value has been set to -1 it can't be flipped back. $\endgroup$ – Steven Stadnicki Mar 19 '15 at 1:05
  • 1
    $\begingroup$ @StevenStadnicki just noticed I have reproduce Section 3 in Systematic computations on Mertens' conjecture and Dirichlet's divisor problem by vectorized sieving referred to in the checked answer. $\endgroup$ – john mangual Mar 19 '15 at 2:18
  • $\begingroup$ You may sieve and find primes at the same time as computing the mobius function $\endgroup$ – qwr Jan 4 at 3:30
  • $\begingroup$ Also I think the primes must be up to $N$, since considering only $p < \sqrt N$ obviously misses all the primes above $\sqrt N$. $\endgroup$ – qwr Jan 4 at 3:50

I also like thisone, $\mathcal{O}(N \ln N)$ additions being not so different to $\mathcal{O}(N \ln \ln N)$ multiplications

% N(1+ln N) additions/affectations. an array of 2 bits elements would be enough

% matlab 
function mu = computeMu(N)
    mu = zeros(1,N); mu(1) = 1; 
    % Will compute the Dirichlet inverse of any sequence starting with 1
    for n = 1:N
        mu(n+n:n:end) = mu(n+n:n:end) - mu(n);
  • $\begingroup$ Very fast and elegant. This doesn't compute the Dirichlet inverse, but rather the Mobius transform $f(x) = x \ast \mu$, right? $\endgroup$ – user76284 Jul 11 at 20:40
  • $\begingroup$ You also only need to go up to N/2+1. $\endgroup$ – user76284 Jul 13 at 0:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.