How to calculate these totient summation sums efficiently?

Question

I am trying to find good ways to tackle sums of the form

$\sum_{k=1}^{N}k^j\varphi(k)$

$j$ can be anything but I am largely concerned about cases 0, 1, and 2.

$\varphi(k)$ is the Euler totient function.

Can this be done without needing to calculate $k^j\varphi(k)$ manually for every single step of $k$? Is there any optimization opportunity? Any identities that apply here that might help?

Answer 1

$\phi$ is multiplicative with $\operatorname{Id} = \phi * \mathbf 1 = \mathbf 1 * \phi$, i.e.

$$n = \sum_{d|n} \phi(d) = \sum_{d|n} \phi \left( \frac n d \right)$$

so the Dirichlet hyperbola method suggests how we can rewrite the sum of $\phi$ in terms of easier functions $\mathbf 1$ and $\operatorname{Id}$.

Starting with the sum of $\operatorname{Id}(m)$ as the triangle numbers $T(n)$:

$$T(n) = \sum_{m=1}^n m = \sum_{m=1}^n \sum_{d|m} \phi \left( \frac m d \right) $$

The trick is now to exchange order of summations. Writing out this sum in rows for $n=6$: \begin{matrix} & d=1 & d=2 & d=3 & d=4 & d=5 & d=6 \\ m=1: & \phi(1) \\ m=2: & \phi(2) & \phi(1) \\ m=3: & \phi(3) & & \phi(1) \\ m=4: & \phi(4) & \phi(2) & & \phi(1) \\ m=5: & \phi(5) & & & & \phi(1) \\ m=6: & \phi(6) & \phi(3) & \phi(2) & & & \phi(1) \end{matrix}

The first column, $d=1$, has totients $\phi(m)$. The second column, $d=2$, has totients $\phi(m/2)$ for $m$ divisible by 2. The third column, $d=3$, has totients $\phi(m/3)$ for $m$ divisible by 3, etc. Formally, by summing over the columns,

$$ \newcommand{\floor}[1]{\left\lfloor #1 \right\rfloor} \sum_{m=1}^n \sum_{d|m} \phi \left( \frac m d \right) = \sum_{m=1}^n \sum_{d=1}^n [d | m] \phi \left( \frac m d \right) = \sum_{d=1}^n \sum_{m=1}^n [d | m] \phi \left( \frac m d \right) = \sum_{d=1}^n \sum_{m=1}^{\floor{n/d} } \phi(m) $$

Let $\Phi(n) = \sum_{m=1}^n \phi(m)$. This is the starting point of all the recursive algorithms: $$ T(n) = \sum_{m=1}^n \Phi \left(\floor{\frac n m}\right) $$

By a generalization to the Möbius inversion formula,

$$\Phi(n) = \sum_{m=1}^n \mu(m) \ T \left( \floor{\frac n m} \right)$$

Since $\mu(m)$ for range $[1,n]$ can be computed in linear time using a linear sieve, this gives a linear time and space algorithm. A small variation on the Sieve of Eratosthenes computes $\mu$ and primes in $O(n \log \log n)$ time.

However, the hyperbola shape of $\floor{n/m}$ suggests a sub-linear algorithm. Rearrange to solve for $\Phi(n)$: $$\Phi(n) = T(n) - \sum_{m=2}^n \Phi \left(\floor{\frac n m} \right)$$

The observation is that for large $m$ (consider $m \ge \sqrt n$), $\lfloor n/m \rfloor$ is constant for many values. We can calculate precisely how many times each $\Phi( k )$ value occurs.

$$ \floor{\frac n m} = \begin{cases} 1 & \text{if} \floor{n/2} < m \le n \\ 2 & \text{if} \floor{n/3} < m \le \floor{n/2} \\ \vdots \end{cases} $$

Using this observation, we only recurse $O(\sqrt n)$ times:

$$\Phi(n) = T(n) - \sum_{m=2}^{\floor{\sqrt n}} \Phi \left(\floor{ \frac n m } \right) - \sum_{k=1}^{\floor{n / \floor{\sqrt n}} -1} \left( \floor{\frac n k} - \floor{ \frac n {k+1} } \right) \Phi(k) $$

$k$ ranges from $1$ up to (but not including) the point where $\Phi(k)$ would overlap with $\Phi(\floor{n/m})$, which is at most one off from $\floor{\sqrt n}$.

The top-level call $\Phi(N)$ requires recursively computing and memoizing $\Phi(\floor{N/m})$ for $2 \le m \le \sqrt N$, and $\Phi(k)$ for $1 \le k \le \sqrt N$. For $\Phi(\floor{N/m})$ term, all recursive calls don't use any values not needed for $\Phi(N)$, as $\floor{\floor{N/m}/m'} = \floor{N/mm'}$. The same applies to recursive calls to $\Phi(k)$, as those contribute directly to $\Phi(N)$. So in total, to evaluate the about $2 \sqrt N$ terms $\Phi(\floor{N/m})$, each requires about $O(\sqrt{N/m})$ operations, assuming constant time memoized storage and lookup. The total computation is bounded by something like $\sum_{m=2}^\sqrt N O(\sqrt{N/m}) = O(N^{3/4})$. See this answer by Erick Wong for more details.

See this blog post by adamant for a more detailed explanation and details on using pre-computed prefix sums for an $O(n^{2/3})$ algorithm. The idea is that you can compute $\Phi$ and prefix sums up to $n^{2/3}$ in linear time with a linear sieve. Then we only need $\Phi(\floor{n/m})$ for $m$ up to $\sqrt[3]{n}$. In total the sums are computed in $\sum_{m=1}^{\sqrt[3]{n}} \sqrt{n/m} = O(n^{2/3})$. This is a little better than $O(n^{3/4})$ but at this point constant factors matter.

Answer 2

For the case $j=0$, you can define some auxiliary summations to formulate an algorithm that runs in $O(n^{3/4})$ time:

$$F(N) = \lvert \{ a,b : 0 < a < b \le N \} \rvert$$

$$R(N) = \lvert \{ a,b : 0 < a < b \le N, \gcd(a,b) = 1 \} \rvert$$

You can see that we are looking for $R(N) + 1$. Also, $F(N)$ is $\dfrac{N(N-1)}{2}$.

Now observe something nice:

R$\left( \Big\lfloor\dfrac{N}{m}\Big\rfloor \right)$ = $\lvert \{ a,b : 0 < a < b \le N, \gcd(a,b) = m \} \rvert$

Why? This is because you can multiply every coprime pair of $(a,b)$ by $m$.

This fact lets you write $F$ in terms of $R$:

F(N) = $\displaystyle\sum_{m=1}^N{ R\left(\Big\lfloor\dfrac{N}{m}\Big\rfloor\right) } $

Since we are looking for $R(N)$, we solve for the first term of the right summation.

$R(N) = F(N) - \displaystyle\sum_{m=2}^N{ R\left(\Big\lfloor\dfrac{N}{m}\Big\rfloor\right) } $

Note this interesting property of the floor function here: $\Big\lfloor\dfrac{N}{m}\Big\rfloor$ will stay constant for a range of $m$. This lets us calculate the summation in chunks. example:

$\Big\lfloor\dfrac{1000}{m}\Big\rfloor$ is constant for $m$ in the range of [501,1000].

Here's a program I wrote in C++ that caches R to trade O(log n) memory for a large speedup

Answer 3

You should be able to calculate all of the values $\phi(1),\dots,\phi(N)$ simultaneously in time $O(N\log\log N)$, assuming you have sufficient memory.

To do so, set up a Sieve of Eratosthenes-type calculation, but instead of only recording whether every integer is prime or not, keep track of each step in the sieve that "crosses off" a given integer. The result is that you will have stored the list of all primes dividing $n$, for all $1\le n\le N$. (You can modify the sieve to get the complete factorization, but it's not important for this problem.)

Once you have this list in storage, you can calculate all the $\phi(n)$ by the formula $\phi(n)=n\prod_{p\mid n}(1-1/p)$. The total number of multiplications is $\sum_{n\le N} \omega(n)$, where $\omega(n)$ is the number of distinct prime factors of $n$; this sum is known to be $O(N\log\log N)$. (I'm sloppily counting a multiplication of two rational numbers as $1$ step.)

A similar setup will allow you to compute the values of any multiplicative function over an interval in time not much longer than the length of the interval.

Answer 4

Andy already described pretty well how to calculate $R_0$ in sublinear time. It is important to memoize results of already calculated $R_j$ in order to avoid duplicate computations. Here is a general formula for any $j\geq0$:

$$R_j(n) = \sum\limits_{k=1}^n k^j \varphi(k) $$

$$S_j(n) = \sum\limits_{k=1}^n k^j$$ $$R_j(n) = S_{j+1}(n) - \sum\limits_{k=2}^n k^j R_j\left ( \left \lfloor \frac{n}{k} \right \rfloor \right )$$

Notice that $S_j$ is just a polynomial of degree $j+1$: $$S_1(n) = \frac{1}{2}n(n+1)$$ $$S_2(n) = \frac{1}{6}n(n+1)(2n+1)$$ $$S_3(n) = \frac{1}{4}n^2(n+1)^2$$

Note that the inner sum can be calculated in $O(\sqrt{n})$ steps (without taking into account calculation of values of $R_j$ which have to be calculated anyway) if we take $q=\lfloor \sqrt{n}\rfloor $:

$$\sum\limits_{k=2}^n k^j R_j\left ( \left \lfloor \frac{n}{k} \right \rfloor \right )$$ $$= \sum\limits_{k=2}^{\lfloor n/q \rfloor} k^j R_j\left ( \left \lfloor \frac{n}{k} \right \rfloor \right ) + \sum\limits_{m=1}^{q-1} \sum\limits_{k=\lfloor \frac{n}{m+1} \rfloor + 1}^{\lfloor \frac{n}{m} \rfloor} k^j R_j(m)$$ $$= \sum\limits_{k=2}^{\lfloor n/q \rfloor} (S_j(k) - S_j(k-1)) R_j\left ( \left \lfloor \frac{n}{k} \right \rfloor \right ) + \sum\limits_{m=1}^{q-1} \left ( S_j \left ( \left \lfloor \frac{n}{m} \right \rfloor \right ) - S_j \left ( \left \lfloor \frac{n}{m+1} \right \rfloor \right ) \right ) R_j(m)$$

The total complexity is $O(n^{3/4})$, but this can be easily improved to $O(n^{2/3})$ if we preprocess values smaller than $O(n^{2/3})$ with a sieve.

Answer 5

In the reply by plamenko to Andy above, we must be careful that R(n) is not the same as $R_0(n)$. I don't really understand why Andy didn't range his $(a,b)$ over $1\leq a\leq b\leq n$, instead of $a<b$. If he had, his proof would carry over just the same, and the summation would give $F(n)=S_1(n)=n(n+1)/2$.

The proof principle set up by Andy is valid for higher $j$, however, with plamenko's notation and [cond] returning 1 if cond is true and 0 otherwise, if you express $R_j(n)$ as $$R_j(n)=\sum_{1\leq a\leq b\leq n} [gcd(a,b)=1] b^j,$$ then you see that, by switching variables $a=mc,b=md$: $$R_j(\lfloor \frac nm\rfloor) = \sum_{1\leq c\leq d\leq \lfloor \frac nm\rfloor} [gcd(c,d)=1] d^j = \sum_{1\leq a\leq b\leq n} [gcd(a,b)=m] \frac{b^j}{m^j},$$ from which, summing over all $m$ like in Andy's proof, we get: $$\sum_{m=1}^n m^j R_j(\lfloor \frac nm\rfloor) = R_j(n) + \sum_{m=2}^n m^j R_j(\lfloor \frac nm\rfloor) = \sum_{m=1}^n \sum_{1\leq a\leq b\leq n} [gcd(a,b)=m] b^j = \sum_{1\leq a\leq b\leq n} b^j$$ which in the end resolves to $\sum_{1\leq b\leq n} b^{j+1} = S_{j+1}(n).$ So that proves plamenko's equation. The rest of plamenko's answer is unchanged, and very useful stuff indeed!

Answer 6

Actually I had a different idea. The number of elements in the Farey sequence of order $N$ is $1+\sum_{n=1}^N \phi(n)$. And one can recursively construct the entire Farey sequence in order (see the first displayed equation). So you might even be able to do this in time $O(N)$.