calculating Möbius function

Question

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

Answer 1

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

Answer 2

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.

Answer 3

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.

Answer 4

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

Answer 5

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

Answer 6

Here's one way to do it. My indexing might not be the best (I modified it so that the first element of the list would be n=1).

import math
def f(a,b):
  if a>b/2:
    return -1
  if abs(a)<=b/2:
    return 0
  if -a>b/2:
    return 1 
def mobius(n):
  mobiuslist=[]
  counter=-1/3
  m=0
  while m<n:
    mobiuslist.append(f(counter, -math.log(1-3**(-m-1))/(m+1)))
    counter-=mobiuslist[m]*math.log(1-3**(-m-1))/(m+1)
    m+=1
  print(mobiuslist)
mobius(100)

In python, this gets the correct value up until n=34 or so due to rounding errors. This runs with O(n) operations (However, I don't know the time complexity of exponentiation or logarithms, so it may well be more than that!). This works because of something related to Lambert Series. Namely, if one replaces the q^n/(1-q^n) in the definition with ln(1-q^n), and works for a little bit, they can get the above. On the other hand, one can copy this method for regular Lambert series to get something similar without the logs. I believe that that would look something like this:

import math
def f(a,b):
  if a>b/2:
    return -1
  if abs(a)<=b/2:
    return 0
  if -a>b/2:
    return 1 
def mobius(n):
  mobiuslist=[]
  counter=-1/3
  m=0
  while m<n:
    mobiuslist.append(f(counter, 1/(3**(m+1)-1)))
    counter+=mobiuslist[m]/(3**(m+1)-1)
    m+=1
  print(mobiuslist)
mobius(100)

Oddly enough, this fails at n=35 (I'm guessing that this is because of rounding errors). I hope that this helps!