GCDLIM - Editorial

editorial
, , ,
1

PROBLEM LINK:

Practice
Contest: Division 1
Contest: Division 2
Contest: Division 3
Contest: Division 4

Author: raysh07
Tester: sushil2006
Editorialist: raysh07

DIFFICULTY:

Easy-Medium

PREREQUISITES:

Harmonic Sum

PROBLEM:

Let P be the prefix gcd array of A. f(A) is defined as number of indices such that P_i \ne P_{i - 1}.

g(N) is defined as expected value of f(A) when A is selected from arrays of size N and 1 \le A_i \le M. Find \lim_{N \to \infty} g(N).

EXPLANATION:

First of all, we only need to care about arrays such that \gcd(A_1, A_2, \ldots, A_N) = 1. To prove this, we can see that if the array ever contains 1, it automatically has a gcd of 1; and the probability of not containing 1 is (\frac{M - 1}{M})^N which tends to 0 as N \to \infty.

Now, we can simplify our computation. Note that because the final gcd is always 1, we can translate f(A) to the number of distinct prefix gcds except for 1. Thus, using Linearity of expectation, we can say f(A) = \sum_{x = 2}^{M}f(x) where f(x) denotes the probability there exists a prefix gcd = x.

To have a prefix gcd of x, obviously x | A_1 is a necessary condition. Suppose that holds, then it is only possible to not have a prefix gcd of x when the prefix gcd jumps from a multiple of x (say y) to not a multiple of x. Note that the probability of this event is \dfrac{M - \frac{M}{y}}{M - \frac{M}{x}}.

Further the “skipping gcd = x” case for all multiples are disjoint events, so we can just add them all and subtract the probability from the probability that x | A_1.

Hence, f(x) = P(x | A_1) - \sum_{y : x | y} f(y) \cdot \dfrac{M - \frac{M}{y}}{M - \frac{M}{x}}.

This can be computed in O(\frac{M}{x}) per x, which sums to O(M log M) by harmonic series.

TIME COMPLEXITY:

\mathcal{O}(M log M) per testcase.

CODE:

Editorialist's Code (C++) 
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF (int)1e18

mt19937_64 RNG(chrono::steady_clock::now().time_since_epoch().count());

// Removed Mint class for readability 

const int mod = 998244353;

void Solve() 
{
    int m; cin >> m;
    
    vector <mint> dp(m + 1);
    for (int i = m; i >= 2; i--){
        dp[i] = mint(m / i) / m;
        
        for (int j = 2 * i; j <= m; j += i){
            int tot = m - (m / j);
            int bad = m - (m / i);
            
            dp[i] -= dp[j] * bad / tot;
        }
    }
    
    mint ans = 0;
    for (int i = 2; i <= m; i++){
        ans += dp[i];
    }
    cout << ans << "\n";
}

int32_t main() 
{
    auto begin = std::chrono::high_resolution_clock::now();
    ios_base::sync_with_stdio(0);
    cin.tie(0);
    int t = 1;
    // freopen("in",  "r", stdin);
    // freopen("out", "w", stdout);
    
    cin >> t;
    for(int i = 1; i <= t; i++) 
    {
        //cout << "Case #" << i << ": ";
        Solve();
    }
    auto end = std::chrono::high_resolution_clock::now();
    auto elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(end - begin);
    cerr << "Time measured: " << elapsed.count() * 1e-9 << " seconds.\n"; 
    return 0;
}
2

this question can also be done using precomputed mobius values and then running a state dp that stores expected value of gcd changes for current prefix gcd=x.about the dp transition think it as dp=1+summationoverk(such that k<x) of{conditional probability of gcd=k given gcd not equal to x * dp[k]}. this simplifies to dp=1+A/(m-gif(m/x)) where x runs 2 to m and A is contribution of cnt*dp[x/u] here cnt is no. of coprimes to u less then lim where u is a divisor of x and lim is gif(m/(x/u) . this no. of coprimes can be found using mobius transform … learn how to compute this from codeforces blog by [Nisiyama_Suzune’s blog] on mobius transform.Now the answer is simply sum over all dps over m mod…
attaching sample code
`#include <bits/stdc++.h>
using namespace std;
using ll=long long;
const int MOD = 998244353;
const int MAX = 200005;
vector divisors[MAX];
int mu[MAX];
//fmexp
ll modpow(ll b, ll e) {
ll res = 1;
while (e) {
if (e & 1) res = res * b % MOD;
b = b * b % MOD;
e >>= 1;
}
return res;
}
//min
ll minv(ll n) {
return modpow(n, MOD - 2);
}
//precompmu
void pre() {
mu[1] = 1;
vector primes;
vector isPrime(MAX, true);

for (int i = 2; i < MAX; i++) {
    if (isPrime[i]) {
        primes.push_back(i);
        mu[i] = -1;  
    }
    for (int p : primes) {
        if ((ll)i * p >= MAX) break;
        isPrime[i * p] = false;
        if (i % p == 0) {
            mu[i * p] = 0;
            break;
        } else {
            mu[i * p] = -mu[i];
        }
    }
} 
for (int i = 1; i < MAX; i++) {
    for (int j = i; j < MAX; j += i) {
        divisors[j].push_back(i);
    }
}

}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
pre();
int t;
cin >>t;
while (t–) {
int M;
cin>>M;
vector dp(M+1,0);
ll tot=0;
for (int x=2; x<=M;x++) {
ll num=0;
for (int u:divisors) {
if (u==1) continue;
int k=x/u;
int lim=M/k;
ll cnt=0;
for (int g:divisors[u]) {
if (mu[g]!=0){
cnt+=mu[g](lim/g);
}
}
num+=cntdp[k];
}

        ll den= M-M/x;
        dp[x]=(1+num%MOD*minv(den)%MOD+MOD)%MOD;
        tot=(tot+dp[x])%MOD;
    }
    ll asn=tot*minv(M)%MOD;
    cout<<asn<<"\n";
}
return 0;

}`

3

i think the transition probability is inverse of what is actually written in the editorial