PROBLEM LINK:
Author: Piyush Kumar
Tester: Kevin Atienza
Editorialist: Kevin Atienza
PREREQUISITES:
Euler’s totient function, Euclidean algorithm, gcd
PROBLEM:
Given A, B and N satisfying A > B \ge 1, \text{gcd}(A,B) = 1, and N \ge 1, compute the following:
QUICK EXPLANATION:
The sum is equivalent to:
Let
Then the answer is F_A(N) - F_B(N), and we also have the following:
Such a sum can be split into two parts, where \lfloor N/g \rfloor \ge s and \lfloor N/g \rfloor < s, where s = \lfloor \sqrt{N} \rfloor, so each sum can be computed in O(N^{1/2 + \varepsilon}) time.
By precomputing S(n) for n = 1 \ldots s and (S(\lfloor N/g \rfloor) for g = s \ldots 1, F_a(N) (and thus the answer) can be computed in O(N^{3/4}) time. With the help of a sieve for smaller values, this can be reduced to O(N^{2/3} (\log \log N)^{1/3}).
EXPLANATION:
Let’s study \text{gcd}(A^i - B^i, A^j - B^j) first, because obviously we won’t be able to reach anywhere without simplifying this somehow. (The numbers involved are very large!) Without loss of generality, assume i \ge j.
You might notice that (A - B) is a common factor of A^i - B^i and A^j - B^j because the following is true for any n:
.
However, (A - B) is not necessarily the gcd! For one, if both i and j are even, then by the same principle (but replaceing A by A^2 and B by B^2), A^2 - B^2 is a common factor too. In fact, we can generalize this:
If x divides y, then A^x - B^x divides A^y - B^y. Thus, we find that if g is the greatest common divisor of i and j, then A^g - B^g is a common divisor of A^i - B^i and A^j - B^j. But is this the greatest?
Let’s run the Euclidean algorithm on \text{gcd}(A^i - B^i, A^j - B^j) to find out. The main property that makes this algorithm true is the following:
Nice property: For any x, y and s, we have \text{gcd}(x,y) = \text{gcd}(x-sy,y)
This holds true even for nonpositive x, y and s, although one has to be careful with interpreting \text{gcd}(x,y). (For example, define \text{gcd}(x,y) as the positive gcd. There are other elegant ways to handle this.)
Euclidean’s algorithm is just an instance of this, with s = \left\lfloor x/y \right\rfloor. Noting that x \bmod y = x - \left\lfloor x/y \right\rfloor y, we have:
The second one is just the base case.
So let’s use the nice property:
What can we say about B^j and A^j - B^j? Note that we can use the nice property again to get \text{gcd}(B^j, A^j - B^j) = \text{gcd}(B^j, (A^j - B^j) + B^j) = \text{gcd}(B^j, A^j). But it’s given that A and B are coprime, so we know that \text{gcd}(B^j, A^j - B^j) = \text{gcd}(B^j,A^j) = 1. Thus, B^j and A^j - B^j have no common factors other than 1, so we find that \text{gcd}(B^j (A^{i-j} - B^{i-j}), A^j - B^j) is just \text{gcd}(A^{i-j} - B^{i-j}, A^j - B^j). Thus, we have the following identity:
Applying this rule s times, we get:
Notice the similarity with the nice property! If we let s = \left\lfloor i / j \right\rfloor, we get the following:
It feels like we’re taking the gcd of the exponents instead. In fact, A^0 - B^0 = 0, so even the base case holds: \text{gcd}(A^i - B^i, A^0 - B^0) = \text{gcd}(A^i - B^i, 0) = A^i - B^i, which means that, if we let F(x) = A^x - B^x, then:
It really feels like we’re doing the Euclidean gcd on the exponents! In fact, using these two, we can prove by induction that:
Thus, we have just shown that if g is the gcd of i and j, then A^g - B^g is the greatest common divisor of A^i - B^i and A^j - B^j!
Now, we can proceed with the sum.
We can separate this into two sums:
Let’s define F_a(N) = \sum_{1 \le i, j \le N} a^{\text{gcd}(i,j)}. Thus, the answer is F_A(N) - F_B(N), so we can just focus on computing F_a(N) for a given a and N.
Computing F_a(N) quickly
To compute, F_a(N) = \sum_{1 \le i, j \le N} a^{\text{gcd}(i,j)}, we can instead loop among all possible gcds g, and count the number of pairs 1 \le i, j \le N whose gcd is g. in other words:
But if \text{gcd}(i,j)=g, then i and j are multiples of g, so let i = i'g and j = j'g, to get:
Define two new functions:
So that F_a(N) is just \sum_{g=1}^N a^g S(\left\lfloor N/g \right\rfloor). Thus, we want to compute the $S(\left\lfloor N/g \right\rfloor)$s quickly.
The former, S(N), is just the number of coprime pairs in the [1,N]\times[1,N] grid of lattice points, while T(N) is just the number of pairs all in all, which is just N^2. However, by again looping among all possible gcds g, we also have:
This gives us the following recurrence:
Amazingly, we can convert this into an algorithm to compute S(N) in o(N) (sublinear) time! To do so, notice that there aren’t very many distinct values possible for \left\lfloor N/g \right\rfloor. Indeed, when g \ge \sqrt{N}, N/g becomes \le \sqrt{N}, so after the first \sqrt{N} terms, there are only \sqrt{N} distinct values left for \left\lfloor N/g \right\rfloor to take. Thus, there are at most 2\sqrt{N} distinct values of \left\lfloor N/g \right\rfloor!
To use this to simplify the sum, we group together the summands with the same \left\lfloor N/g \right\rfloor value. First, we split it into two halves. Let s = \left\lfloor \sqrt{N} \right\rfloor. Then:
The first summation only has at most \sqrt{N} terms, so it can be computed naïvely. Let’s focus on the second sum:
We can instead iterate over all distinct values of \left\lfloor N/g \right\rfloor.
Which g s satisfy \left\lfloor N/g \right\rfloor = v? Simple manipulation gives the following list:
There are \left\lfloor \frac{N}{v} \right\rfloor - \left\lfloor \frac{N}{v+1} \right\rfloor such numbers. Therefore, the sum above is equal to:
Such a sum can be computed in O(s) = O(\sqrt{N}) time! Thus, S(N) can be computed in O(\sqrt{N}) time assuming we have the necessary smaller values of S.
By first computing S(n) for n = 1, 2, \ldots, \sqrt{N} and then S(\left\lfloor N/g \right\rfloor) for g = \left\lfloor \sqrt{N} \right\rfloor, \left\lfloor \sqrt{N} \right\rfloor -1, \ldots, 1 using this sum, one can compute S(N) quickly. (And indeed, all S(\left\lfloor N/g \right\rfloor) together with it.) One can work out that the running time of this algorithm is O(N^{3/4}), which is sublinear in N!
What about F_a(N) itself? Remember that:
We can group the summands into two cases similarly as above, but this time, the second sum is more complicated:
If we proceed similarly as above, we will get the following equivalent sum:
Thankfully, the sequence of g s satisfying \left\lfloor N/g \right\rfloor = v are consecutive integers, as we have seen above. And summing a^g for a list of consecutive integers can be done with a simple geometric series formula:
Such a number can be computed with fast modular exponentiation and computing the modular inverse of (a-1). (If a = 1, then a-1 is not invertible, but in that case we simply have F_a(N) = N^2.) Thus, F_a(N) can be computed in O(\sqrt{N} \log N) time! (Assuming we have already computed the $S(\left\lfloor N/g \right\rfloor)$s.)
The overall running time is O(N^{3/4}), dominated by the computation of the $S(\left\lfloor N/g \right\rfloor)$s.
Sieve
Once can improve this running time by sieving the small values of S. First, remember that:
Now, consider T(N) - T(N-1). Let’s look at T(N-1) first:
We have allowed g to go up to N because \left\lfloor \frac{N-1}{g} \right\rfloor will become 0 anyway. Thus:
Now, \left\lfloor \frac{N-1}{g} \right\rfloor is usually equal to \left\lfloor \frac{N}{g} \right\rfloor. The only times this doesn’t happen are when g is a divisor of N, in which case we have \left\lfloor \frac{N-1}{g} \right\rfloor = \left\lfloor \frac{N}{g} \right\rfloor - 1. Thus:
We can extract S(N) - S(N-1) from that last equality to get:
This equality allows us to create a sieve to compute S(n) for n = 1 \ldots U quickly, as shown below (noting that T(N) - T(N-1) is just 2N - 1):
S = array[1..U]
for n in 1..U:
S[n] = 2*n - 1
for g in 1..U:
for n in 2g..U by g:
S[n] -= S[g]
S[g] += S[g-1]
After running this code, S[n] now contains S(n) for n = 1 \ldots U. And since for each g we are only iterating through multiples of g, the running time is O(U/1 + U/2 + U/3 + \ldots + U/U) = O(U \log U).
To compute all S(\left\lfloor N/g \right\rfloor) quicker than O(N^{3/4}), we first choose a value for U, then compute S(n) for n = 1 \ldots U using the sieve above, and finally compute S(\left\lfloor N/g \right\rfloor) for \left\lfloor N/g \right\rfloor > U using our original O(\sqrt{N/g}) algorithm above. It turns out that the optimal choice for U is \Theta\left(\left(\frac{N}{\log N}\right)^{2/3}\right), which gives us an overall running time of O(N^{2/3} (\log N)^{1/3}).
Finally, we mention that in fact there’s a slightly faster sieve than the above, which runs in O(U \log \log U) time rather than O(U \log U). With such a sieve, and with a different choice for U, the running time becomes O(N^{2/3} (\log \log N)^{1/3}). We leave finding this sieve to the reader. (Hint: Express S(n) using Euler’s totient function.)
Time Complexity:
O(N^{3/4}), O(N^{2/3} (\log N)^{1/3}), or O(N^{2/3} (\log \log N)^{1/3})