Summing the totient function $\sum_{k=1}^n \varphi(k)$

Question

I explored some computer science/number theory challenges sites for fun, and they presented the following problem, exactly as follows:

Let $$P(n) = \sum_{k=1}^n \varphi(k)$$

Find $P(10^{16})$

I searched for quite a while about this and tried different approaches:

1. Using the formula for $$\varphi(n)= n \cdot \prod_{i=1}^k \frac{p_i-1}{p_i}$$ I tried to calculate each $\varphi(n)$ in range, but this becomes very inefficient for large $n$. I could get as far as $10^7$ with this approach. Beyond this it just gets too slow.

2. I tried a different one, more direct. Wikipedia and Wolfram Alpha suggest similiar formulas for directly calculating $P(n)$: $$P(n) = \sum_{k=1}^n \varphi(k)= \frac12 \cdot \biggl (1+ \sum_{k=1}^n\mu (k)\cdot \lfloor {\frac nk} \rfloor^2\biggl)$$ This formula seemed a lot more promising. I tried it and managed to get alot further than $10^7$ but still far from the target. With pre-calculating a sieve of the Moebius function, I could get to a bit less than $10^9$. My memory was insufficient, and couldn't compute anymore values in a sieve. And even if I could, it still takes a long time and is very far from $10^{16}$.

Here is part of the code that I used for my 2nd approach written in Java:

public static BigInteger PhiSummatoryFunction (long limit)
{
    BigInteger sum = BigInteger.ZERO;
    int [] m = MoebiusSieve(limit);
    for (int i=1;i<m.length;i++)
        sum=sum.add(BigInteger.valueOf((long) (m[i]*Math.floor(limit/i)*Math.floor(limit/i))));
    return sum.add(BigInteger.ONE).divide(BigInteger.ONE.add(BigInteger.ONE));
}

Where MoebiusSieve is a function that computes the Moebius function values up to a certain limit in a sieve, using an eratosthenes-like method.

---

3. After understanding and implementing the recursive method suggested in a link provided in the comments: $$P(n) = \frac {n(n+1)}{2} - \sum_{i=2}^\sqrt n P(\lfloor \frac ni \rfloor) - \sum_{j=1}^\sqrt n P(j) \cdot (\lfloor \frac nj \rfloor - \lfloor \frac n{j+1} \rfloor)$$

I can compute values up to $P(10^{11})$, and with maximum memory allocation, pre-computing as many $\varphi(n)$ as possible and consequently all $P(n)$ that I can for memoization, I can compute $P(10^{12})$ in just over 20 minutes. A major improvement but still a little far from $P(10^{16})$. It's ok if the computation takes a bit longer, but I fear $P(10^{16})$ would take exponentially longer time, judging by the "jump" in computation time between $P(10^{11})$ and $P(10^{12})$. My memory allows me to "save" up to $350,000,000 \space φ(n)$ values OR up to $700,000,000 \space μ(k)$ values. Perhaps there is a way to perform the summation using μ(k) values rather than φ(n)?.

All my computations suggest and show that my recursion is the prominent time consumer. This is obvious, but I am sure it takes longer than it should, as pointed out by qwr. So I am posting below the recursion code, with some documentation. It seems to me that this is the right way to do this computation, but my implementation is not optimal.

    public static BigInteger phiR (long limit, long [] s) // limit is 10^t, s is the sieve of precomputed values of `P(n)`. Can store maximum 350,000,000 values
{                                                                                                                                               
    if (limit<s.length)                                 
        return BigInteger.valueOf(s[(int) limit]);
    BigInteger sum = BigInteger.valueOf(limit).multiply(BigInteger.valueOf(limit).add(BigInteger.ONE)).divide(BigInteger.valueOf(2)); // this corresponds to the n'th triangular number
    BigInteger midsum1=BigInteger.ZERO;
    BigInteger midsum2=BigInteger.ZERO;
    for (long m=2;m*m<=limit;m++) // computing the first sum, first for changing floor(limit/m) values
        midsum1=midsum1.add(phiR((long) Math.floor(limit/m),s));
    for (long d=1;d*d<=limit;d++) // computing the second sum
        if ((double)d!=Math.floor(limit/d))
            midsum2=midsum2.add(BigInteger.valueOf((long) (Math.floor(limit/d)-Math.floor(limit/(d+1)))).multiply(phiR(d,s)));
    sum=sum.subtract(midsum1).subtract(midsum2);
    return sum;
}

I was suggested to use dictinaries by qwr, in addition to the array, for big values of $n$, but I don't know anything about it. Can another improvement be made to make the time frame a day or so?

---

Answer 1

The recursive solution you showed takes $O(n^{2/3})$ time (see my answer). To the best of my knowledge this is as fast as you can get for exact answers using recursive approaches. My implementation in PyPy computes $P(10^{12})$ in 22 seconds (extrapolating $P(10^{16})$ will take 1-2 days with enough memory), so you'll need to work on the efficiency of your implementation.

My algorithm is almost a straightforward translation from my linked answer. Currently I have a cutoff set at $n^{2/3}$. For all $n$ less than this cutoff, the totients/totient sums are pre-calculated using a Sieve of Eratosthenes approach. For values of $n$ larger than the cutoff, compute recursively and store in a dictionary. For $P(10^{16})$ you will probably run out of memory trying to store all totient values and instead have to load and store from disk.

You can scour the references in OEIS A002088 for resources which I often find useful.

On OEIS, only terms up to $P(10^{18})$ are known (A064018). The highest values were computed by Hiroaki Yamanouchi (cursory research indicates he may be the legendary Min_25 on Project Euler, SPOJ, CodeChef, etc.!) on OEIS; you may try contacting him to see how he did it.

Answer 2

Asymptotically you may be able to do better using the approach described in Lagarias and Odlyzko, Computing $\pi(x)$: an analytic method (1987) J. Algorithms vol.8 pp.173-191. In section 2 they discuss conditions which suffice to adapt the approach to functions other than $\pi(x)$, and $\Phi(x)$ would seem to meet them. However, it should be noted that although asymptotically their algorithm is $O(x^{1/2 + \varepsilon})$, the hidden constant may be rather higher than the algorithm you're currently using.