# How to calculate these totient summation sums efficiently?

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?

• Is $\phi$ Euler's totient function? It might be a good idea to specify that in the question (even though it's implied in the title), because the symbol $\phi$ is not reserved exclusively for this function. Feb 28, 2013 at 0:26
• It is, I'll edit the OP Feb 28, 2013 at 0:40
• Even for $j=0$, I don't know of any way to calculate $\sum^N\phi(k)$ without finding all the numbers and adding them. There are asymptotic estimates for these sums --- are they of any use, or do you really need the exact numbers? Feb 28, 2013 at 0:43
• For what it's worth, $\sum k\phi(k)$ is tabulated at oeis.org/A011755 Feb 28, 2013 at 0:46
• @GerryMyerson But my question is if there is a fast way to get all the numbers. An analogy I might use is that it's faster to use a sieve to generate primes than it is to check each number if it's prime, etc. I am curious if there's some sort of "quicker shortcut" to get the same results (in perhaps less than O(N) time) Feb 28, 2013 at 1:11

$$\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.

• This is genius. Apr 30, 2019 at 13:56

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

• This seems to be taken from Project Euler problem 73 PDF.
– qwr
Oct 7 at 21:59

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.

• A Google search took me to this page: hackage.haskell.org/packages/archive/NumberSieves/0.1.2/doc/… ... I can't personally verify its effectiveness though. Feb 28, 2013 at 2:10
• I know there's a faster way than O(N log log N), but I am struggling to get there. I've already got the O(N log log N) method working. Feb 28, 2013 at 2:38
• If you had an $O(N\log\log N)$ method working, it would've been nice if you'd mentioned that in your post. Feb 28, 2013 at 6:58
• I'm inclined to believe that it's not necessary to keep track of a list of all primes. The product formula with $(1 - 1/p)$ can be calculated in-place.
– qwr
Apr 12, 2016 at 5:31
• I can confirm my $O(n)$ memory modification in my previous comment. Also, using integer division and subtraction, floating-point computations can be avoided.
– qwr
Apr 15, 2016 at 11:27

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.

• +1, Very nice, and quite efficient. How would you modify this algorithm if, for example, you only wanted the totient sum for odd numbers $\leq n$? May 13, 2015 at 13:35

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!

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)$.

• Wouldn't enumerating the Farey sequence using a recursive technique be of order of the length of the Farey Sequence, \sum_{n=1}^N \phi(n) and not of order N?
– Dan
Jun 26, 2014 at 21:52
• hmm, I think you're right. drat. Jun 27, 2014 at 2:03