YETGCD - Editorial [Code Uncode 5.0]

editorial
, , ,
1

PROBLEM LINK

Contest
Practice

Setter: dedsec_29
Tester: satyam_343
Editorialist: dedsec_29

PREREQUISITES

Number theory, Euler Totient Function

PROBLEM

Find the value of the following summation modulo 10^9+7

\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{k=1}^{i}\sum_{l=1}^{j}\gcd (k,i)\gcd(l,j)\gcd(i,j)

EXPLANATION

I will introduce some frequently used notations and lemmas first.

  • [P] refers to the boolean expression, i.e. [P] = 1 when P is true, and 0 otherwise.
  • d|n means that d can divide n (without a remainder).
  • Let us define f\displaystyle(x) = \sum_{i=1}^{x}\gcd (x,i)
  • Harmonic sum lemma:

We can rewrite the problem as
\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}f(i)f(j)\gcd(i,j)

Using the divisor sum property \sum_{d|n} \phi(d) = n, we can say that
\displaystyle \gcd(i,j) = \sum_{d|\gcd(i,j)} \phi(d) \displaystyle= \sum_{d=1}^{N} [d|i][d|j]\phi(d)

Let’s rewrite the summation using this simplification now.
\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}f(i)f(j)\gcd(i,j)
= \displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}f(i)f(j)\sum_{d=1}^{N} [d|i][d|j]\phi(d)
= \displaystyle\sum_{d=1}^{N}\phi(d)\sum_{i=1}^{N}\sum_{j=1}^{N}[d|i]f(i)[d|j]f(j)
= \displaystyle\sum_{d=1}^{N}\phi(d)\sum_{i=1}^{N}[d|i]f(i)\sum_{j=1}^{N}[d|j]f(j)
= \displaystyle\sum_{d=1}^{N}\phi(d)(\sum_{i=1}^{N}[d|i]f(i))^2
= \displaystyle\sum_{d=1}^{N}\phi(d)(\sum_{d|i}f(i))^2
= \displaystyle\sum_{d=1}^{N}\phi(d)(\sum_{i=1}^{\lfloor{\frac{n}{d}}\rfloor}f(i\cdot d))^2

which can be calculated using two loops in O(nlogn) time using the harmonic sum lemma, given that we can precompute values of f(i).

Now it only remains to show how to calculate f(n) = \sum_{i=1}^{n}\gcd (n,i) values efficiently.
Let’s count the number of i, 1 \leq i \leq n such that [\gcd(n,i) = d]
[\gcd(n,i) = d]\implies \gcd(n/d,i/d)=1
The count of such i is \phi(n/k), therefore f(n) = \sum_{d|n} d \cdot\phi(n/d)
\phi(n) can be precomputed in O(nloglogn) time, and f(n) for all n can be computed by iterating over all divisors in O(nlogn) time. Hence overall time complexity is O(nlogn).

Calculation of phi(n) 
for (int i = 1; i <= n; i++)
        phi[i] = i;
for (int i = 2; i <= n; i++) {
        if (phi[i] == i) {
            for (int j = i; j <= n; j += i) 
                phi[j] -= phi[j]/i;
        } 
    }
Calculation of f(n) 
for (ll i = 1; i <= n; i++) {
        for (int j = i; j <= n; j += i)
            f[j] = ((phi[j/i] * i) % mod + f[j]) % mod;
    }

This problem can also be solved using Principle of Inclusion-Exclusion (refer tester’s solution)

Setter's solution 
#include <bits/stdc++.h>
using namespace std;
#pragma GCC optimize "trapv"
#define ll long long

ll const mod = 1e9+7;
int const siz = 1e6;
ll f[siz+1], phi[siz+1];

void precomputation(int n) {
    for (int i = 1; i <= n; i++) 
        f[i] = phi[i] = 0;
    phi[0] = 1;
    for (int i = 1; i <= siz; i++)
        phi[i] = i;
    for (int i = 2; i <= siz; i++) {
        if (phi[i] == i) {
            for (int j = i; j <= siz; j += i) 
                phi[j] -= phi[j]/i;
        } 
    }
    for (ll i = 1; i <= siz; i++) {
        for (int j = i; j <= siz; j += i)
            f[j] = ((phi[j/i] * i) % mod + f[j]) % mod;
    }
}

void solve() {
    int n; cin >> n;
    precomputation(n);
    ll ans = 0;
    for (int i = 1; i <= n; i++) {
        ll sum = 0;
        for (int j = i; j <= n; j += i) 
            sum = (sum + f[j]) % mod;
        sum = sum * sum % mod;
        ans = (ans + (sum * phi[i]) % mod ) % mod;
    }
    cout << ans << "\n";
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);
    solve();
    return 0;
}
Tester's solution 
#pragma GCC optimize("O3")
#pragma GCC target("popcnt")
#pragma GCC target("avx,avx2,fma")
#pragma GCC optimize("Ofast,unroll-loops")
#include <bits/stdc++.h>   
#include <ext/pb_ds/tree_policy.hpp>
#include <ext/pb_ds/assoc_container.hpp>
using namespace __gnu_pbds;
using namespace std;
#define ll long long 
const ll INF_MUL=1e13;  
const ll INF_ADD=1e18;
#define pb push_back               
#define mp make_pair        
#define nline "\n"                           
#define f first                                          
#define s second                                             
#define pll pair<ll,ll> 
#define all(x) x.begin(),x.end()   
#define vl vector<ll>           
#define vvl vector<vector<ll>>    
#define vvvl vector<vector<vector<ll>>>          
#ifndef ONLINE_JUDGE    
#define debug(x) cerr<<#x<<" "; _print(x); cerr<<nline;
#else
#define debug(x);  
#endif     
void _print(ll x){cerr<<x;}  
void _print(char x){cerr<<x;}   
void _print(string x){cerr<<x;}     
mt19937 rng(chrono::steady_clock::now().time_since_epoch().count()); 
template<class T,class V> void _print(pair<T,V> p) {cerr<<"{"; _print(p.first);cerr<<","; _print(p.second);cerr<<"}";}
template<class T>void _print(vector<T> v) {cerr<<" [ "; for (T i:v){_print(i);cerr<<" ";}cerr<<"]";}
template<class T>void _print(set<T> v) {cerr<<" [ "; for (T i:v){_print(i); cerr<<" ";}cerr<<"]";}
template<class T>void _print(multiset<T> v) {cerr<< " [ "; for (T i:v){_print(i);cerr<<" ";}cerr<<"]";}
template<class T,class V>void _print(map<T, V> v) {cerr<<" [ "; for(auto i:v) {_print(i);cerr<<" ";} cerr<<"]";} 
typedef tree<ll, null_type, less<ll>, rb_tree_tag, tree_order_statistics_node_update> ordered_set;
typedef tree<ll, null_type, less_equal<ll>, rb_tree_tag, tree_order_statistics_node_update> ordered_multiset;
typedef tree<pair<ll,ll>, null_type, less<pair<ll,ll>>, rb_tree_tag, tree_order_statistics_node_update> ordered_pset;
//--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
const ll MOD=1e9+7;   
const ll MAX=1000000;
vector<int> phi(MAX+5);
void phi_1_to_n(int n) {
    for (int i = 0; i <= n; i++)
        phi[i] = i;

    for (int i = 2; i <= n; i++) {
        if (phi[i] == i) {
            for (int j = i; j <= n; j += i)
                phi[j] -= phi[j] / i;  
        }
    }
}
void solve(){ 
    ll n; cin>>n;
    vector<ll> use(n+5,0);
    phi_1_to_n(MAX);
    for(ll i=1;i<=n;i++){
        for(ll j=i;j<=n;j+=i){
            use[j]+=i*phi[j/i];  
        }
    }
    ll ans=0;
    vector<ll> anot(n+5,0);
    for(ll i=n;i>=1;i--){
        ll opt=0;
        for(ll j=i;j<=n;j+=i){
            opt+=use[j];
        }
        opt%=MOD;
        ll now=opt*opt; now%=MOD;
        for(ll j=i;j<=n;j+=i){  
            now-=anot[j]; now%=MOD;
        }   
        now=(now+MOD)%MOD;
        anot[i]=now;
        ans+=now*i; ans%=MOD;
    }
    cout<<ans;
    return;                     
}                                      
int main()                                                                              
{  
    ios_base::sync_with_stdio(false);                         
    cin.tie(NULL);                           
    #ifndef ONLINE_JUDGE                 
    freopen("input.txt", "r", stdin);                                              
    freopen("output.txt", "w", stdout);  
    freopen("error.txt", "w", stderr);                        
    #endif         
    ll test_cases=1;                  
    //cin>>test_cases;
    while(test_cases--){     
        solve();       
    }
    cout<<fixed<<setprecision(10);
    cerr<<"Time:"<<1000*((double)clock())/(double)CLOCKS_PER_SEC<<"ms\n"; 
}
2

There’s a somewhat more straightforward solution. Note that once you have f(n), you just need to compute the GCD convolution of the array f with itself and it can be done in time O(n \log n) (for instance, Library Checker). To compute f(n), note that it is the Dirichlet convolution of \phi with \text{id}, so it is a multiplicative function. A closed form solution for prime powers can be easily computed, and using a linear sieve, this can be done in O(n \log n).

1 Like