Very large table of the Möbius function

Question

The Möbius function is important in analytic number theory for many reasons.

I'd like to pre-compute a big table of values of the Möbius function to test a few things (sum of Möbius function, i.e. Mertens function, and other similar things, etc.).

Because building a table of Möbius function involves quite long computation (it requires to find the number of prime divisors of each integer):

1. Is there an online resource that provides a table of $\mu(n)$ for n = 1 ... 10^10 ?

2. If not, with which programming language would you do this? Maybe would you have open-source code for this?

Answer 1

I converted this code into Python. For the following tests, I used Python 64-bit because it uses more than 1 GB of RAM.

mobius(10**6) takes 2.5 seconds with pure Python + numpy.

mobius(10**6) takes 0.7 seconds with Python + numpy + numba (just-in-time compiler).

mobius(10**9) takes 39 seconds with Python + numpy + numba.

---

import numpy as np
from numba import jit       # comment this line if you don't have numba installed

@jit                        # idem
def mobius(N):
    SQRT = int(np.sqrt(N))
    mu = np.ones(N+1, dtype=np.int32)
    mu[0] = 0

    for i in range(2, SQRT+1):
        if mu[i] == 1:
            for j in range(1, int(N / i)+1):
                mu[i*j] *= -i
            for j in range(1, int(N / i / i)+1):
                mu[i*i*j] = 0

    for i in range(2, N+1):
        if mu[i] == i:
            mu[i] = 1
        elif mu[i] == -i:
            mu[i] = -1
        elif mu[i] < 0:
            mu[i] = 1
        elif mu[i] > 0:
            mu[i] = -1

    return mu

print mobius(10**6)

# Test that shows that it works, see https://oeis.org/A008683
OEISmobius = [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,-1,-1,0,1,1,1,0,-1,1,1,0,-1,-1,-1,0,0,1,-1,0,0,0,1,0,-1,0,1,0,1,1,-1,0,-1,1,0,0,1,-1,-1,0,1,-1,-1,0,-1,1,0,0,1,-1]
assert (mobius(78) == OEISmobius).all()

---

Edit: I also translated this code in Python:

import numpy as np
from numba import jit

@jit
def mobius2(N):
    mu = -np.ones(N+1, dtype=np.int8)
    mu[0] = 0
    mu[1] = 1
    for n in range(2, N+1):
        for m in range(n+n, N, n):
            mu[m] -= mu[n]
    return mu

print mobius2(10**9)   

mobius2(10**9) takes ? seconds with Python + numpy + numba.

Answer 2

A method along the lines of the sieve of Eratosthenes will do what you want quite efficiently. See https://mathoverflow.net/questions/99473/calculating-möbius-function for more details, references and sample code. You will need 64-bit integer arithmetic (or arbitrary magnitude integer arithmetic as provided in Python) to get up to $10^{10}$, but that shouldn't be a problem on modern hardware. You will also need either enough memory to hold the whole table or do something a bit more clever than the sample code in the linked answer.

Answer 3

This may help: Gevorg Hmayakyan wrote a short paper which gives a recursive formula for the $\mu(n)$.

https://ghmath.files.wordpress.com/2010/06/mobius.pdf

It's not for the faint of heart, however. But maybe a faster algorithm could come out of it.

Answer 4

B. Goddard makes reference to an interesting paper by Gevorg Hmayakyan on calculating the mobius function. The following source code is in c++ with armadillo and implements one of his equations. The source code leaves a lot to be desired as it involves multiplication of large polynomials etc. I am yet to implement his final equation. At present the code can only handle numbers up to 1000.

void numbertheory::Mobius_Hmayakyan(uword n)
 {
 uword i,j;       //counters
 uword x, y;      // size variables
 ivec Mz(1);      //current Hmayakyan polynomial
 Mz(0)=-1;        // set Hmayakyan polynomial for i=4 Mz=-1
 ivec Mz_1;       //previous Hmayakyan polynomial
 ivec lambda1;     // polynomial ..λ_i-1 = (1+z+z^2+...+z^i-1)
 uvec lambda2;    // polynomial ..λ_i-2 = (1+z+z^2+...+z^i-2)
 uvec phi;      // poynomial phi_3 = λ_1*λ_2...*λ_i-3
 uvec phi_past = {1,1};     // poynomial phi_4 = λ_1*λ_2...*λ_i-4,
 sword mu=0;      //set mobius function for mu(4)=0
 if( n >=1)
    cout << " mu(1) =  1 " <<endl;
 if( n >=2)
    cout << " mu(2) = -1" <<endl;
 if( n >=3)
    cout << " mu(3) = -1" <<endl;
 if(n >=4 )
    cout << " mu(4) =  0" <<endl;

 i=5;
 while(i <= n)
 {
// copy vector Mz to vector Mz_1 and update size of vector Mz and set to
//zero
    Mz_1.zeros(Mz.size());
    Mz_1= Mz;
    Mz.zeros(1+ ((i-2)*(i-3)/2));

 //Generate Mahonian polynomial phi(z)
   lambda2.ones(i-2);
   phi = conv(phi_past, lambda2);

//copy vector phi to vector phi_past
  phi_past.zeros(phi.size());
  phi_past = phi;

 // multiply phi by z by shifting upwards by one
  ivec B;
  B.zeros(phi.size()+1);
  for(j=0;j< phi.size(); j++)
     B(j+1)=phi(j);


 lambda1.ones(i-1);           // Generate polynomial (1+z+z^2+....z^i-2)
 ivec A = conv(Mz_1,lambda1);   //Multiply previous Hmayakyan polynomials
                               //Mz_1 and lambda by convolution

 //resize vectors to same size and pad with zeros at high end
// in order to perform subtraction
  x = A.size();
  y = B.size();
  if(x >y)   B.insert_rows(y, x-y);  
  if(x <y )  A.insert_rows(x, y-x);

// Update the Hmayakyan polynomial according to the following formulae
 // M_i(z) = M_i-1(z).λ_i−2(z) − μ(i-1). z.ψ_i−3(z)
 Mz = A - mu*B;

//retrieve mobius function from current Hmayakyan polynimial
 mu = Mz((i-2)*(i-3)/2);

 cout << " mu(" << i << ") = " << mu << endl;

 i++;
 }

 }

Answer 5

The second approach of Gevorg Hmayakyan may be of interest [see his comments] as it suggests that the mobius function can be determined from the sign of the coefficients of the pentagonal numbers. Gevorg stated that if $$ ∏^{n−1}_{i=1}(1−x^i)=∑^{n(n−1)/2}_{k=1}a_{n,k}x^k$$ then $$\mu(n)=a_{n,1}+a_{n,n+1}+a_{n,2n+1}+a_{n,3n+1}+...$$ I notice that this function is equivalent to a row of Euler Truncated Product Triangle (see the paper by Alex Mennen at http://www.alex.mennen.org/mahoniantri.pdf). Using Mennens notation we can then state $$\mu(n)=\sum_{k=0}^{\frac{n(n-1)}{2}}P(n-2,kn+1)$$ where $$P(n-2,c)=p(c)+\sum_{k_{1}=0}^{c-n}p(k_{1})+\sum_{k_{2}=0}^{2k_{1}-c}p(k_{2})+\sum_{k_{3}=0}^{2k_{2}-k_{1}}p(k_{3})+\sum_{k_{4}=0}^{2k_{3}-k_{2}}p(k_{4})...$$ and where in accordance to Eulers Pentagonal theorem $p(k)=(-1)^{a}$ if there is some integer $a$ such that $a(3a-1)/2=k$ otherwise it is $0$.

I have generated recursive code based on this function but it highly unsatisfactory. It gives reliable results but the recursion increases quadratically with multiple calls to the same steps and crashes after about n equals 30. One thing of note is that any $\sum_{k=0}^n{p(k)}$ is either $+1,0,-1$. But the function $P(n-2,c)$ is not that simple.

Perhaps someone can take it one step further.