PROBLEM LINK:
Author: Kirill Shchetinin
Tester: Istvan Nagy
Editorialist: Xellos
DIFFICULTY:
EASY
PREREQUISITES:
Factorisation using the sieve of Erathostenes.
PROBLEM:
You’re given many values of N; for each of them, compute \sum_{K=1}^N \frac{N}{gcd(N,K)}.
QUICK EXPLANATION:
Compute the first few terms by hand, find a simple formula in OEIS. Precompute the smallest prime divisors, use them to factorise N and use that formula to compute the answer.
EXPLANATION:
The most efficient solution to our problem uses Google Search Algorithm. Just compute the first few values by hand and search for this sequence in OEIS. And look, it’s here! There are two things to take from the page, which we’ll call Result 1 and Result 2:
-
the answer is the sum of d\cdot\varphi(d) for all divisors d of N
-
the answer is the product of \frac{p^{2e+1}+1}{p+1}=\frac{(p^e)^2\cdot p+1}{p+1} for all maximum prime powers p^e dividing N
But it sucks to just state those results and provide a link to proofs - let’s prove it all here as well!
Result 1 can be seen quite easily: the number of terms with N/gcd(N,j)=d (each term must be a divisor of N, of course) is the Euler totient \varphi(d). That’s because we can rewrite this equality to N/d=k=gcd(N,j)=gcd(kd,j), which can only hold if j is a multiple of k - that is, j=lk for 0 < l \le d. We can divide both sides by k and get 1=gcd(d,l); the number of possible l-s coprime to d and \le d (which is the same as the number of j-s for which N/d=gcd(N,j)) is \varphi(d).
In order to prove Result 2, we’ll use Result 1. If N=\prod_j p_j^{e_j} (the product goes over prime divisors of N), then any divisor d can be written as \prod_j p_j^{\alpha_j} for any 0 \le \alpha_j\le e_j. The formula for \varphi(d) is \prod_j p_j^{\alpha_j-1}(p_j-1) for \alpha_j > 0; it can also be written as \prod_j \varphi(p_j^{\alpha_j}), where the factor for \alpha_j=0 is \varphi(1)=1.
Let’s take some number N > 1 (for N=1, the answer is a product of nothing, which is 1). It’ll be divisible by some prime p; let’s cut out the max. power of this prime p^e which divides N and write N=kp^e. Then, we can take the sum over divisors to be the sum over exponents \alpha of p^\alpha and divisors of k:
where the first sum is
The second sum is just the answer for N/p^e (note that p^{2e+1} can overflow, for example if N=p\approx10^7). By continually cutting out primes this way, we eventually arrive at N=1 and obtain the factors of the answer as in Result 2.
It’s not very hard to precompute divisors of all numbers up to N=10^7 (there are about N\log N divisors of numbers up to N), then precompute all their totients \varphi and simulate the sum over divisors, but with our limits, it’s too slow. In fact, we can just barely manage to compute one prime divisor of each number using a modified sieve of Erathostenes - instead of just marking numbers as composite as in the original sieve, we’re marking one prime for which they were found to be composite (or primes as their own prime divisors):
D[1..N] ={fill with zeroes} # one prime divisor of each number
for i in 2..N:
if D[i] == 0: # prime
# mark it as a divisor for all its multiples
for j in 1..N/i:
D[i*j] = i
This works in O(N\log{N}) time as well, but with a good constant factor.
With Result 2, however, all we need to do is decompose any N into prime powers. With the precomputed primes, that’s easy - as long as N > 1, we can just keep taking one prime divisor p of N, divide N by it as many times as possible and we’ve got p^e, with which we’ll compute one factor of the answer. Since N decreases quickly when we do this - we always divide it by at least 2, so we’ll reach N=1 in at most O(\log{N}) divisions - and computing the factors of our answer as per Result 2 is really fast, this is sufficient to solve the problem.
I’m not sure if I’d be able to make a fast enough solution using only Result 1. Sometimes, optimisation isn’t a good idea…
ALTERNATIVE SOLUTION:
Could contain more or less short descriptions of possible other approaches.
AUTHOR’S AND TESTER’S SOLUTIONS:
The author’s solution can be found here.
The tester’s solution can be found here.
The editorialist’s solution can be found here.
