C-Code for Merten's Function

The C-Code for Merten's function proceeds much in the same way as that of Totient Summatory function and uses the formula derived in Merten's function. This code evaluates the value of $M(10^{11})$ in approx. $8$ seconds and $M(10^{12})$ in approx. $80$ seconds. Here is the code.

#include<stdio.h>
#include<math.h>

#define M 60000000
#define Mod 1000000007

int MertCache[M+1], Mert[M+1];
char Mu[M+1], flag[M+1];

void precalc()
{
    int i, j;
    Mu[1] = 1;
    Mert[1] = 1;
    MertCache[1] = 56789;
    for(i = 2; i <= M; i++) {
        if(!flag[i])
            for(j = i; j <= M; j += i) {
                if( !flag[j] )
                    Mu[j] = 1;
                if(j/i%i == 0)
                    Mu[j] = 0;
                else
                    Mu[j] = -Mu[j];
                flag[j] = 1;
        }
        Mert[i] = Mert[i-1] + Mu[i];
        MertCache[i] = 56789;
    }
}

long long Merten(long long n, long long k)
{
    long long m = n/k;
    if(m <= M)
        return Mert[m];
    else if(MertCache[k] != 56789)
        return MertCache[k];
    long long p = sqrt(m), q = m/(p+1), i, res = 0;
    for(i = 2; i <= p; i++)
        res += Merten(n, k*i);
    for(i = 1; i <= q; i++)
        res = res + (m/i) * Mu[i];
    res = 1 - res;
    res += p * Merten(n, k*(p+1));
    return MertCache[k] = res;
}

int main()
{
    precalc();
    printf("%lld", Merten(100000000000LL,1));
    return 0;
}

Note: $\text{Merten}(n,k) = M\left(\left\lfloor\dfrac{n}{k}\right\rfloor\right)$ in the code.

This code doesn't include the Phi function at all. But if you plan to use both the Totient summatory function and the Merten's function in the same code, I think code may slow down. I'vnt tried it yet but I'll also post the code that uses both the functions.

UPDATE (13 May 2017): As in Totient Summatory function, I made an iterative version of the Merten's function. You can check the performance of the code in this Online IDE.

from timeit import default_timer

start = default_timer()

n = 10 ** 6
s = int(n ** 0.5)


def sieve():
    Lows, Flag = [0] + [1] * s, [0] * (1 + s)
    Lows[1] = 1
    for i in range(2, 1 + s):
        if Flag[i] == 0:
            for j in range(i, 1 + s, i):
                Flag[j] = 1
                if (j // i % i) == 0:
                    Lows[j] = 0
                else:
                    Lows[j] *= -1
        Lows[i] += Lows[i - 1]
    return Lows


Lows = sieve()
Highs = [0] * (1 + s)

for k in range(s, 0, -1):
    tempn = n // k
    u = int(tempn ** 0.5)
    v = tempn // (1 + u)
    res = 1
    for i in range(1, 1 + v):
        res -= (Lows[i] - Lows[i - 1]) * (tempn // i)
    for i in range(2, 1 + u):
        res -= Highs[k * i] if k * i <= s else Lows[tempn // i]
    i = 1 + u
    res += Lows[tempn // i] * u
    Highs[k] = res

print(Highs[1])
print(default_timer() - start)

Note: $\text{Highs}[k] = M\left(\left\lfloor\dfrac{n}{k}\right\rfloor\right)$ in the code. Also, the value $-962982$ is purely arbitrary just to make sure that value has not been touched yet.


Until then
Yours Aye
Me