I found this problem when I was trying to find a hard one related to summative functions, but surpresingly I don’t expect that I’ll really need to apply so many math and programming concepts here. My first tought was to solve it with previous methods that I use for other problems, but none of them works. The intuition to solve this is simple, and the recursive formula it’s easy for simple cases. I find a sub-linear algorithm to solve this problem, but I don’t know how to calculate a better bound for my time complexity.
Statement
Let $S(n,m) = \sum \varphi(n\times i)$ for $1 \leq i \leq m$. ($\varphi$ is Euler’s totient function)
You are given that $S(510510,10^6)=45480596821125120$.
Find $S(510510, 10^{11})$. Give the last $9$ digits of your answer.
Solution
Naive solution
Let’s try with the simplest and naive solution, just doing a for:
long long phi(long long n){
long long x = n;
for(int i=2; i*i<=x; i++){
if(x%i == 0){
while(x%i == 0) x /= i;
n -= n/i;
}
}
if(x != 1) n -= n/x;
return n;
}
long long S(int n, long long m){
long long x = 0;
for(long long k=1; k<=1e11; k++){
x += phi(n*k);
}
return x;
}
Let’s analyze the time complexity of each function, calculate $\varphi(n)$ in a naive way could take in the worst case up to $\mathcal{O}(\sqrt{n})$. The function to calculate $S(n,m)$ it’s just a simple for:
\begin{aligned} \mathcal{O}(S(n, m)) &= \sum_{k=1}^{m} \mathcal{O}(\varphi(n k))\cr &= \sum_{k=1}^{m} \mathcal{O}(\sqrt{nk})\cr &= \mathcal{O}(m\sqrt{nm}) \end{aligned}
As you can see, it’s a terrible time complexity for the constraints of the problem. You can try to compute it, but it will take you at least a few years to finish it.
You can try to improve the complexity calculating efficiently the function $\varphi(n)$, but even if you can do that in constant time, your complexity will be $\mathcal{O}(m)$. So, a common thinking for this problem is don’t iterate over $m$ to avoid the linear behavior.
Sub-linear method
I already solved some kind of these problems before, so my first tought was to apply the hyperbola method with the dirichlet notation to see if I can reduce to something like to $\mathcal{O}(n^{2/3})$, but I can’t find a nice way to do that.
My second attempt was using the Mobius function, I have more progress with this, but in some point I get stuck and don’t know how to procedure, so again, this idea was useless.
Complex problems requires simple problems, it makes some sense. I mean, let’s try with simpler cases.
$n$ is prime
Take a look to the function $S(n, m)$ when $n$ is a prime number, to simplify the notation I will use the letter $p$:
\begin{aligned} S(p, m) &= \sum_{k=1}^{m} \varphi(pk)\cr &= \varphi(p) + \varphi(2p) + \ldots + \varphi(pm)\cr &= \varphi(p)\Big(\varphi(1) + \varphi(2) + \ldots + \varphi(p) + \ldots + \varphi(m) \Big)\cr &+ (\varphi(p^2)-\varphi(p)\varphi(p)) + \ldots + (\varphi(p\cdot p \cdot \lfloor m/p\rfloor)-\varphi(p)\varphi(p\cdot \lfloor m/p\rfloor))\cr &= \varphi(p) \Phi(m) + \varphi(p) + \varphi(2p) + \ldots + \varphi(p\cdot \lfloor m/p \rfloor)\cr &= \varphi(p)\Phi(m) + S(p, \lfloor m/p \rfloor) \end{aligned}
This is a really good recursive definition. Let’s analyze its complexity:
\begin{aligned} \mathcal{O}(S(p, m)) &= \mathcal{O}(\Phi(m)) + \mathcal{O}(S(p, m/p))\cr &= \log_{p}(m) \mathcal{O}(\Phi(m))\cr &= \mathcal{O}(m^{2/3}\log_{p} m) \end{aligned}
I explained in a previous post how to calculate $\Phi(n)$ in $n^{2/3}$. As you can see, this is a sub-linear solution and could end the task in few seconds, but sadly, $n$ is not a prime number.
$n$ is a prime power
Let’s analyze a more general case, when $n = p^{k}$, for some prime number $p$ and for some power $k>1$:
\begin{aligned} S(p^{k}, m) &= \sum_{i=1}^{m} \varphi(p^{k} i)\cr &= \varphi(p^{k}) + \varphi(2p^{k}) + \ldots + \varphi(p^{k}m)\cr &= \varphi(p^{k})\Big(\varphi(1) + \varphi(2) + \ldots + \varphi(p) + \ldots + \varphi(m) \Big)\cr &+ (\varphi(p^{k+1})-\varphi(p^{k})\varphi(p)) + \ldots + (\varphi(p^{k}\cdot p \cdot \lfloor m/p\rfloor)-\varphi(p^{k})\varphi(p\cdot \lfloor m/p\rfloor))\cr &= \varphi(p^{k}) \Phi(m) + \varphi(p^{k}) + \varphi(2p^{k}) + \ldots + \varphi(p^{k}\cdot \lfloor m/p \rfloor)\cr &= \varphi(p)\Phi(m) + S(p^{k}, \lfloor m/p \rfloor) \end{aligned}
the recursion is similar to the previous one, but this formula could work for numbers like $2^{37}$, which is greater than $10^{11}$, but it’s not what the problem requires, so we need to continue searching the general formula.
$n$ is a product of two primes
Now, let’s analyze the case when $n = p_{1}p_{2}$, where $p_{1}$ and $p_{2}$ are different prime numbers.
\begin{aligned} S(p_{1}p_{2}, m) &= \sum_{i=1}^{m} \varphi(p_{1}p_{2}i) \cr &= \varphi(p_{1}p_{2})\Big( \varphi(1) + \varphi(2) + \ldots + \varphi(m) \Big) \cr &+ (\varphi(p_{1}^{2}p_{2}) - \varphi(p_{1})\varphi(p_{1}p_{2})) + (\varphi(p_{1}p_{2}^{2}) - \varphi(p_{2})\varphi(p_{1}p_{2})) - (\varphi(p_{1}^{2}p_{2}^{2}) - \varphi(p_{1}p_{2})\varphi(p_{1}^{2}p_{2}^{2})) + \ldots\cr &= \varphi(p_{1}p_{2})\Phi(m) + S(p_{1}p_{2}, \lfloor m/p_{1} \rfloor) + S(p_{1}p_{2}, \lfloor m/p_{2} \rfloor) - S(p_{1}p_{2}, \lfloor m/p_{1}p_{2} \rfloor) \end{aligned}
General case
I don’t have a proof for the next formula, but since the case for two primes correspond to apply the inclusion - exclusion principle, I think that makes sense to think that in the general case it also happens:
\begin{aligned} S(n, m) &= \varphi(n)\Phi(m) + \sum S(n, \lfloor m/p_{i}\rfloor) - \sum S(n, \lfloor m/p_{i}p_{j}\rfloor) + \ldots \end{aligned}
Let’s analyze the time complexity:
\begin{aligned} \mathcal{O}(S(n, m)) &= \mathcal{O}(m^{2/3}) + \mathcal{O}(\sum S(n, m/p_{i}))\cr &= \mathcal{O}(m^{2/3}) \times \mathcal{O}(f(\log(n)/\log\log(n), m)) \end{aligned}
Now, let me explain what is $f$. To calculate $S(n,m)$, I’m iterating over each possible subset of the prime factors of $n$, a nice upper bound for this concept is $\log(n)/\log \log(n)$. So, for a numer $n$, it’s expected to have $\log(n)/\log \log(n)$ number of distinct prime factors. Iterate over them takes $2^{size}$, but applying memoization to the function $S(n,m)$ means that it will only iterate over the distinct values of $\lfloor m/P\rfloor$, where $P$ is the product of some subset of the prime factors. It’s hard to say what is the time complexity, but to simplify, we can consider that we’re iterating over all possible subsets in each recursive call:
\begin{aligned} \mathcal{O}(S(n, m)) &< \mathcal{O}(2^{\log(n)/\log\log(n)} \log(m) m^{2/3}) \end{aligned}
my code finds the solution for the problem in less than 2 seconds.
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using t128 = __int128;
#define N 30000000
#define MOD 1000000000
int phi[N];
int lpf[N];
int Phi[N];
vector <int> f;
vector <int> v;
map <ll, ll> Phi_cache;
map <ll, ll> S_cache;
void cphi(){
phi[1] = 1;
for(int i=2; i<N; i++){
if(!phi[i]){
phi[i] = i-1;
lpf[i] = i;
for(ll j=2*i; j<N; j+=i){
if(phi[j] == 0) phi[j] = j;
phi[j] = phi[j]/i*(i-1);
if(lpf[j] == 0) lpf[j] = i;
}
}
}
for(int i=1; i<N; i++) Phi[i] = (Phi[i-1] + phi[i])%MOD;
}
vector <int> prime_factors(int n){
vector <int> ans;
while(n > 1){
int d = lpf[n];
while(n%d == 0) n /= d;
ans.push_back(d);
}
return ans;
}
ll T(ll n){
t128 ans = n; ans *= (n+1);
ans /= 2;
ans %= MOD;
return ans;
}
ll calc_Phi(ll n){
if(n < N) return Phi[n];
if(Phi_cache.find(n) != Phi_cache.end()) return Phi_cache[n];
ll ans = T(n);
ll sqrt_n = floor(sqrt(n));
for(int m=1; m < n/sqrt_n; m++){
ans -= (1LL*Phi[m]*(n/m - n/(m+1)))%MOD;
ans %= MOD;
}
for(int m=2; m<=sqrt_n; m++){
ans -= calc_Phi(n/m);
ans %= MOD;
}
ans += MOD; ans %= MOD;
Phi_cache[n] = ans;
return ans;
}
ll S(int n, ll m){
if(m == 0) return 0;
if(m == 1) return phi[n];
if(S_cache.find(m) != S_cache.end()) return S_cache[m];
int sz = f.size();
ll ans = (1LL*phi[n]*calc_Phi(m))%MOD;
for(int i=1; i<(1<<sz); i++){
int c = __builtin_popcount(i);
if(c%2 == 0) ans -= S(n, m/v[i]);
else ans += S(n, m/v[i]);
ans %= MOD;
}
if(ans < 0) ans += MOD;
S_cache[m] = ans;
return ans;
}
void fill_mask(){
int sz = f.size();
v.resize((1<<sz));
v[0] = 1;
for(int i=1; i<(1<<sz); i++){
int previous_mask = v[i&(i-1)];
int pos = __builtin_ctz(i&(-i));
int val = previous_mask*f[pos];
v[i] = val;
}
}
int main(){
cphi();
int n; ll m; cin >> n >> m;
f = prime_factors(n);
fill_mask();
cout << S(n, m) << '\n';
return 0;
}