18
$\begingroup$

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

http://en.wikipedia.org/wiki/M%C3%B6bius_function

$\endgroup$

5 Answers 5

15
$\begingroup$

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."

$\endgroup$
12
$\begingroup$

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.

$\endgroup$
5
  • $\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, 2012 at 2:04
  • 1
    $\begingroup$ So, how many different spellings of Eratosthenes can we accumulate in this thread? $\endgroup$ Jun 14, 2012 at 5:21
  • $\begingroup$ @Gerry: Thanks. Fixed here ... and in all my other documents. $\endgroup$ Jun 14, 2012 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$ Mar 23, 2014 at 4:33
  • 1
    $\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$ Mar 23, 2014 at 14:55
9
$\begingroup$

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.

$\endgroup$
7
  • $\begingroup$ I want to avoid factoring it. I'm also considering Sieve of Eratosthenes $\endgroup$ Jun 13, 2012 at 16:53
  • 21
    $\begingroup$ Calling any algorithm the "stupidest possible" underestimates the power of human ingenuity.... $\endgroup$ Jun 13, 2012 at 17:34
  • 2
    $\begingroup$ @Greg: I double dog dare ya :) $\endgroup$
    – Igor Rivin
    Jun 13, 2012 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$ Jun 13, 2012 at 18:03
  • 4
    $\begingroup$ Is erathsosphenes what you get when you cross Eratosthenes with Aristophanes? $\endgroup$ Jun 13, 2012 at 23:14
7
$\begingroup$

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]))
P[100:105]

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

547, 557, 563, 569, 571

$\endgroup$
4
  • 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$ Mar 19, 2015 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$ Mar 19, 2015 at 2:18
  • $\begingroup$ You may sieve and find primes at the same time as computing the mobius function $\endgroup$
    – qwr
    Jan 4, 2019 at 3:30
  • 1
    $\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, 2019 at 3:50
4
$\begingroup$

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);
    end
end
$\endgroup$
2
  • $\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, 2019 at 20:40
  • $\begingroup$ You also only need to go up to N/2+1. $\endgroup$
    – user76284
    Jul 13, 2019 at 0:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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