C-code for the SumOfPrimes function

We saw the algorithm for the generalized Prime counting function in the previous post. Here is the C code for the same, with some slight optimizations.


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

#define M 1000000
#define Mod 1000000000

long long Lower[M+1], Higher[M+1];
char Used[M+1];

long long SumNat(long long n)
{
    long long p = n, q = n+1;
    if(p%2)
        q /= 2;
    else
        p /= 2;
    return (p % Mod) * (q % Mod) % Mod;
}

long long SumPrimePi(long long n)
{
    long long v = sqrt(n), p, temp, q, j, end, i, d, t;
    Higher[1] = SumNat(n) - 1;
    for(p = 2; p <= v; p++) {
        Lower[p] = SumNat(p) - 1;
        Higher[p] = SumNat(n/p) - 1;
    }
    for(p = 2; p <= v; p++) {
        if(Lower[p] == Lower[p-1])
            continue;
        temp = Lower[p-1];
        q = p * p;
        Higher[1] -= (Higher[p] - temp) % Mod * p % Mod;
        if(Higher[1] < 0)
            Higher[1] += Mod;
        j = 1 + (p & 1);
        end = (v <= n/q) ? v : n/q;
        for(i = p + j; i <= 1 + end; i += j) {
            if(Used[i])
                continue;
            d = i * p;
            if(d <= v)
                Higher[i] -= (Higher[d] - temp) % Mod * p % Mod;
            else {
                t = n/d;
                Higher[i] -= (Lower[t] - temp) % Mod * p % Mod;
            }
            if(Higher[i] < 0)
                Higher[i] += Mod;
        }
        if(q <= v)
            for(i = q; i <= end; i += p*j)
                Used[i] = 1;
        for(i = v; i >= q; i--) {
            t = i/p;
            Lower[i] -= (Lower[t] - temp) % Mod * p % Mod;
            if(Lower[i] < 0)
                Lower[i] += Mod;
        }
    }
    return Higher[1] % Mod;
}

int main() {
    printf("%lld", SumPrimePi(1000000000000LL));
    return 0;
}


The code computes the value of $\pi_1(10^{12})$ in about 6 seconds which is way, way better than Mathematica's computation time. If you are looking for the C-code for $\pi(x)$, you can find it in this post.