EXPGCD - Editorial

PROBLEM LINKS :

Contest : Division 1

Contest : Division 2

Practice

Setter : Yusuf Kharodawala

Tester : Istvan Nagy

Editorialist : Anand Jaisingh

DIFFICULTY :

Medium

PREREQUISITES :

Number Theory, Inclusion Exclusion Formula, Basic Combinatorics

PROBLEM :

Given 2 integer N and K, we need to find the summation of the GCD of each sub sequence of length K of the sequence 1,2,...N . We need to keep K fixed, and find the answer for each N in the range 1,2,....200000

QUICK EXPLANATION :

The expected value of the number to be calculated shall be \frac{X}{ \binom{N}{K} } . We now consider the subtask of trying to find X.

The number of numbers in the range 1,2,...N that are divisible by some number i is \left\lfloor{\frac{N}{i}}\right\rfloor, so the number of sequences of length K whose GCD is divisible by i is \binom{\left\lfloor{\frac{N}{i}}\right\rfloor}{K} . Then, we can use the inclusion exclusion method to find the number of sequences whose GCD is divisible by i and at least one prime even after all numbers are divided by i. We can speed up this inclusion exclusion using the Mobius Function. This will help us get the number of sequences whose GCD is exactly i. This gives us an efficient method to calculate the answer for a number N in O(\sum_{i=1}^{N} divisors(i)) time. Then, we can prove that we can calculate the answer for some N from the answer N-1 in O(divisors(N)) time by using the fact that \left\lfloor{\frac{N}{i}}\right\rfloor \ne \left\lfloor{\frac{N-1}{i}}\right\rfloor if and only if i is a divisor of N.

So, we can overall calculate the answer over all N in O(\sum_{i=1}^{N} divisors(i)) time.

EXPLANATION :

Keeping N fixed, let’s see how we can solve the problem :

The expected value of the number to be calculated shall be \frac{X}{ \binom{N}{K} } . We now consider the sub task of trying to find X.

Let’s try and invert the given subtask a little bit. Instead of finding the summation of GCD’s of each subsequence of length K, we’ try and find \forall i \hspace{0.2cm} 1 \le i \le N , the number of subsequences of length K having GCD equal to i. For some i, let’s call the so found number f(i).

Then, we can calculate X as \sum_{i=1}^{N} i \cdot f(i) . The main part of this problem now boils down to calculating f(i).

For some i, the number of numbers in the range 1,2,...N divisible by i is exactly \left\lfloor{ \frac{N}{i}} \right \rfloor . These numbers are 1 \cdot i , 2\cdot i , … \left\lfloor{ \frac{N}{i}} \right \rfloor \cdot i

So, assume we want to find the number of subsequenes whose GCD is divisible by i. That number shall be \binom{ \left \lfloor{\frac{N}{i}}\right\rfloor}{K} . If we select any subsequence of length K consisting of only multiples of the number i, then it’s obvious , the GCD of any such subsequence is also divisible by i. Let’s call this G(i)

Now, we an deduce that if the GCD of a subsequence is divisible by i, then the GCD of the subsequence is a multiple of i. . Now, If we subtract from G(i), the number of subsequences of length K such that their GCD is a multiple of i but \ne i . Then we can easily get f(i) from G(i).

So, we need to calculate the number of sub sequences of length K, such that their GCD is divisible by i and is also divisible by at least one prime number even after dividing by i .

Now, you need to read about the Mobius Function here and come back. We make a small change to that function, let \mu(1)=0 . The Mobius function is a faster way of performing inclusion exclusion.

Mathematically speaking , we get :

f(i) = G(i) + \sum_{j=1}^{\left\lfloor{\frac{N}{i}}\right\rfloor} \mu(j) \cdot G(i \cdot j )

And the number X is :

X = \sum_{i=1}^{N} i \cdot f(i) , So

X = \sum_{i=1}^{N} i \cdot G(i) + \sum_{i=1}^{N} \sum_{j=1}^{\left\lfloor{\frac{N}{i}}\right\rfloor} i \cdot \mu(j) \cdot G(i \cdot j )

We can find the first part of the summation easily. So, the second part can be inverted to :

\sum_{z=1}^{N} G(z) \cdot \sum_{d|z } d \cdot \mu( \frac{z}{d}) .

And finally, we get

X = \sum_{i=1}^{N} i \cdot G(i) + \sum_{i=1}^{N} G(i) \cdot \sum_{d|i} d \cdot \mu(\frac{i}{d})

If we precompute the divisor of each function and the mobius function, then we know that for a fixed z , the sum \sum_{d|z} d \cdot \mu(\frac{z}{d}) can be calculated easily. We know G(i) = \binom{\left\lfloor{\frac{N}{i}}\right\rfloor}{i}.

So, after some pre computations, we can use the above formula to calculate the answer for fixed N easily in O(N)

The last part we still have to deal with is calculating the answer for each N in the range 1,2,...200000 . If we can precompute each of these too, then all queries can be answered in O(1).

Let’s call for some z, \sum_{d|z} d \cdot \mu(\frac{z}{d}) = sum[z] . This number remains unaltered regardless of N. Now, we prove by induction that if we have calculated the answer for N-1, then we can calculate the answer for N in O(divisors(N)) time.

We can re-write X as :

X = \sum_{i=1}^{N} i \cdot G(i) + \sum_{i=1}^{N} G(i) \cdot sum[i] . So,

X = \sum_{i=1}^{N} (i+sum[i]) \cdot G(i) .

A crucial observation we make here is that the value \left\lfloor{\frac{N}{i}}\right\rfloor \ne \left \lfloor{\frac{N-1}{i}}\right\rfloor if and only if N is a multiple of i.

So, for some i, G(i) changes only at values of N when N is a multiple of i. So, if we know the answer for N-1, then we can use its answer, and replace G(i) for each divisor of N to find the new answer.

There are many many problems before that can be used to learn how to compute the Mobius Function, the divisors of a number and all factorials and inverse factorials fast. You need to consider these to be pre-requisites, as they are quite common by now.

That’s it, Thank you

Your comments are welcome !

COMPLEXITY ANALYSIS :

Time Complexity : O ( ( \sum_{i=0}^{MaxN} divisors[i] ) +Q ) , where MaxN \approx 200000

Space Complexity : : O ( \sum_{i=0}^{MaxN} divisors[i] ) , where MaxN \approx 200000

SOLUTION LINKS :

Setter
    import java.util.*;
    import java.io.*;
    class Solution{
     
    	final static long mod = (long)1e9+7;
        static long[] fac;
     
        static long combi(int n, int m) {
            long ans = fac[n];
            ans *= power(fac[m], mod - 2, mod);
            ans %= mod;
            ans *= power(fac[n - m], mod - 2, mod);
            return ans %= mod;
        }
     
        static long power(long x, long y, long p){
            long res = 1;
            x %= p;
            while (y > 0) {
                if((y & 1)==1){
                    res *= x;
                    res %= p;
                }
                y >>= 1;
                x *= x;
                x %= p;
            }
            return res;
        }
     
    	public static void main(String[] args){
    		
    		FastReader s=new FastReader(System.in);
    		PrintWriter w=new PrintWriter(System.out);
     
    		int MXSIZE = 200001;
    		int t = s.nextInt(), k = s.nextInt();
    		long[] coprimes = new long[MXSIZE], ans = new long[MXSIZE], p = new long[MXSIZE], q = new long[MXSIZE];
            fac = new long[MXSIZE];
            fac[0] = 1;
            for(int i = 1; i < MXSIZE; i++){		//Calculate factorials
                fac[i] = fac[i - 1] * i;
                fac[i] %= mod;
            }
            ArrayList<Integer>[] factors = new ArrayList[MXSIZE];
            for(int i = 1; i < MXSIZE; i++){
                factors[i] = new ArrayList<>();
                factors[i].add(1);
            }
            for(int i = 2; i < MXSIZE; i++){		//Store Factors
                for(int j = i << 1; j < MXSIZE; j += i) {
                    factors[j].add(i);
                }
            }
            for(int i = k; i < MXSIZE; i++){		//Actual dp solution
                for(int f: factors[i]){	
                    coprimes[i] -= coprimes[i / f]; //Subtract sets whose gcd is NOT 1
                    coprimes[i] += mod;
                    coprimes[i] %= mod;
                }
                coprimes[i] += combi(i - 1, k - 1); //Adding all possible combitions. This is added at last because 1 is a factor, and i don't want it to get subtracted by  
                coprimes[i] %= mod;			 	   //itself in previous loop, which will result in its value becoming zero
                for(int f: factors[i]){
                    p[i] += coprimes[i / f] * f;		//Calculate p and q, p being the sum of gcd and q is total sets
                    q[i] += coprimes[i / f];		//Note although we are finding gcd of all permutations, we need not multiply factorial(k) as 
                    p[i] %= mod;					//it will be multiplied in both p and q, thus effectively cancelling out :)
                    q[i] %= mod;
                }
            }
            //Now coprime[i] contains total sets which have integer 'i' in it and have gcd zero
         	for(int i = k; i < MXSIZE; i++){
         		p[i] += p[i - 1];	//Calculate sum of all GCDs till i
         		q[i] += q[i - 1];	//Calculate sum of all sets till i
         		p[i] %= mod;
         		q[i] %= mod;
         		ans[i] = (p[i] * power(q[i], mod - 2, mod))%mod; //perform modular division
         	}
     
            while(t-->0){
            	w.println(ans[s.nextInt()]); //display the required answer
            }
     
    		w.close();
    	}
     
    	static class FastReader { 
            BufferedReader br; 
            StringTokenizer st; 
      
            public FastReader(InputStream i) { 
                br = new BufferedReader(new InputStreamReader(i)); 
            } 
      
            String next() { 
                while (st == null || !st.hasMoreElements()) { 
                    try{ 
                        st = new StringTokenizer(br.readLine()); 
                    } 
                    catch (IOException  e) { 
                        e.printStackTrace(); 
                    } 
                } 
                return st.nextToken(); 
            } 
      
            int nextInt() { 
                return Integer.parseInt(next()); 
            } 
     
        } 
    } 
Tester

#include <bits/stdc++.h>

#define all(x) (x).begin(), (x).end()
#define rall(x) (x).rbegin(), (x).rend()
#define forn(i, n) for (int i = 0; i < (int)(n); ++i)
#define for1(i, n) for (int i = 1; i <= (int)(n); ++i)
#define ford(i, n) for (int i = (int)(n) - 1; i >= 0; --i)
#define fore(i, a, b) for (int i = (int)(a); i <= (int)(b); ++i)
 
template<class T> bool umin(T &a, T b) { return a > b ? (a = b, true) : false; }
template<class T> bool umax(T &a, T b) { return a < b ? (a = b, true) : false; }
 
using namespace std;
 
const int64_t MOD = 1e9 + 7;
 
int inverse(int a)
{
	int p = MOD - 2;
	int64_t res = 1;
	while (p)
	{
		if (p & 1)
		{
			res = (res * a) % MOD;
		}
		a = (static_cast<int64_t>(a) * a) % MOD;
		p >>= 1;
	}
	return res;
}
 
vector<int> solver(int T, int K, int maxN)
{//vN sorted
	vector<int64_t> binomK(maxN + 1), su(maxN+1), fin(maxN + 1);
	//init binom
	if (K <= maxN)
	{
		binomK[K] = 1;
		for (int i = K + 1; i <= maxN; ++i)
		{
			binomK[i] = ((binomK[i - 1] * i) % MOD) * inverse(i - K) % MOD;
		}
	}
	vector<int32_t> res(maxN + 1);
	for (int i = K; i <= maxN; ++i)
	{
		fin[i] = binomK[i];
		for (int j = 2; static_cast<int64_t>(i / j) * (i/j) >=i ; ++j)
		{
			fin[i] = (fin[i] + MOD - fin[i / j])%MOD;
		}
		for (int j = 2; j *j < i; ++j)
		{
			fin[i] = (fin[i] + (MOD - (i/j -i/(j+1))) * fin[j]) % MOD;
		}
	}
 
	vector<vector<int> > divi(maxN + 1);
	for (int i = 2; i <= maxN; ++i)
	{
		for (int j = 2 * i; j <= maxN; j += i)
		{
			divi[j].push_back(i);
		}
	}
 
	for (int i = K; i <= maxN; ++i)
	{
		su[i] = (su[i - 1] + fin[i] - fin[i - 1] + MOD) % MOD;
 
		for(auto d: divi[i])
		{
			su[i] = (su[i] + d * (MOD + fin[i/d] - fin[(i-1)/d])) % MOD;
		}
		
		
		res[i] = (su[i] * inverse(binomK[i])) % MOD;
	}
	return res;
}
 
int main(int argc, char** argv)
{
#ifdef HOME
	if (IsDebuggerPresent())
	{
		freopen("../in.txt", "rb", stdin);
		freopen("../out.txt", "wb", stdout);
	}
#endif
	int T, N, K;
	scanf("%d %d", &T, &K);
	vector<int> w = solver(1, K, 200001);
	for (uint32_t i = 0; i < T; ++i)
	{
		scanf("%d", &N);
		printf("%d\n", w[N]);
	}
	return 0;
}
Editorialist
//#pragma GCC optimize("Ofast,no-stack-protector")
//#pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx")
//#pragma GCC target("avx,tune=native")
// Anand Jaisingh
 
#include<bits/stdc++.h>
 
using namespace std;
 
typedef complex<double> base;
typedef long double ld;
typedef long long ll;
 
#define pb push_back
#define pii pair<int,int>
#define pll pair< ll , ll >
#define vi vector<int>
#define vvi vector< vi >
const int maxn=(int)(2e5+5);
const ll mod=(ll)(1e9+7);
int fact[maxn],inv_fact[maxn],sum[maxn],mu[maxn],sum2[maxn],pre[maxn];
vector<int> divs[maxn];
bool pr[maxn];
 
inline int mul(ll a,ll b)
{
    return (a*b)%mod;
}
 
inline int add(int a,int b)
{
    int ret=a+b;
 
    if(ret>=mod)
    {
        ret-=mod;
    }
 
    return ret;
}
 
inline int sub(int a,int b)
{
    int ret=a-b;
 
    if(ret<0)
    {
        ret+=mod;
    }
 
    return ret;
}
 
 
inline int poww(ll a,ll b)
{
    int x=1,y=a;
 
    while(b>0)
    {
        if(b%2)
        {
            x=mul(x,y);
        }
 
        y=mul(y,y);b/=2;
    }
 
    return (int)(x%mod);
}
 
void build()
{
    fact[0]=1;
 
    for(int i=1;i<maxn;i++)
    {
        fact[i]=mul(fact[i-1],i);
    }
 
    inv_fact[maxn-1]=poww(fact[maxn-1],mod-2);
 
    for(int i=maxn-2;i>=0;i--)
    {
        inv_fact[i]=mul(inv_fact[i+1],i+1);
    }
 
    fill(mu,mu+maxn,1);mu[1]=0;
 
    for(int i=2;i<maxn;i++)
    {
        if(!pr[i])
        {
            for(int j=i;j<maxn;j+=i)
            {
                pr[j]=true;
 
                mu[j]=mul(mu[j],mod-1);
 
                ll now=i*1ll*i;
 
                if(j%now==0)
                {
                    mu[j]=0;
                }
            }
        }
    }
 
    for(int i=1;i<maxn;i++)
    {
        for(int j=i;j<maxn;j+=i)
        {
            divs[j].pb(i);
        }
    }
 
    for(int i=1;i<maxn;i++)
    {
        for(int x:divs[i])
        {
            sum[i]=add(sum[i],mul(mu[x],i/x));
        }
    }
 
    /*
 
    for(int i=1;i<maxn;i++)
    {
        sum2[i]=add(sum2[i-1],i);
    }
 
    for(int i=1;i<maxn;i++)
    {
        pre[i]=add(pre[i-1],sum[i]);
    }
 
     */
}
 
int C(int n,int r)
{
    if(n<r || r<0)
    {
        return 0;
    }
 
    return mul(fact[n],mul(inv_fact[r],inv_fact[n-r]));
}
 
int P(int n,int k)
{
    if(k<0 || n<k)
    {
        return 0;
    }
 
    return mul(fact[n],inv_fact[n-k]);
}
 
int solve2(int n,int k)
{
    int ret=0;
 
    for(int i=1;i<=n;i++)
    {
        int val1=add(i,sum[i]),val2=C(n/i,k);
 
        ret=add(ret,mul(val1,val2));
    }
 
    return ret;
}
 
int solve3(int n,int k)
{
    int ret=0;
 
    vector<int> vals;vals.resize(n+1);
 
    for(int i=1;i<=n;i++)
    {
        vals[i]=C(n/i,k);
    }
 
    for(int i=n;i>=1;i--)
    {
        for(int j=i+i;j<=n;j+=i)
        {
            vals[i]=sub(vals[i],vals[j]);
        }
    }
 
    for(int i=1;i<=n;i++)
    {
        ret=add(ret,mul(i,vals[i]));
    }
 
    return ret;
}
 
int main()
{
    build();
 
    int q,k;scanf("%d %d",&q,&k);
 
    vector<int> res,prev_val;res.resize(200001);prev_val.resize(200001);
 
    int ans=0;
 
    for(int i=k;i<=200000;i++)
    {
       for(int x:divs[i])
       {
            int now=add(x,sum[x]);
 
            ans=sub(ans,mul(prev_val[x],now));
 
            int new_coff=C(i/x,k);
 
            ans=add(ans,mul(new_coff,now));
 
            prev_val[x]=new_coff;
       }
 
       int den=mul(inv_fact[i],mul(fact[k],fact[i-k]));
 
       res[i]=mul(ans,den);
    }
 
    while(q>0)
    {
        int n;scanf("%d",&n);
 
        printf("%d\n",res[n]);
 
        q--;
    }
 
    return 0;
}
2 Likes

I little more math and we can get, i+sum[i] = phi[i]

Here we assume mu[1]=0,

instead we take the original mobius function, we can reduce i+sum[i] of previous case into sum[i].

Now we can use Dirichlet convolution and obtain sum[i] = idOmu;

Now we know phiOI = id and muOI = en;

therefore idOmu = phiOIOmu = phiOen = phi(n);

Here is my accepted solution.

Refernce : https://www.quora.com/profile/Surya-Kiran/Posts/A-Dance-with-Mobius-Function?fbclid=IwAR1fRllkPF6w4CqkYS9YKwgMOgcYwgOF9DYzFTMzcd5_EHg6-jCY3fgtJQQ

I learned about Mobius function in this challenge, so there might be a chance of errors.

1 Like

https://www.codechef.com/viewsolution/26479707
Can anyone help me to find out , why am I getting WA in basic subtask .

I have a doubt,
f(i) is the number of subsequences of length K having GCD equal to i.
G(i) is the number of subsequenes of length K whose GCD is divisible by i.
Hence, f(i) is a part of G(i), so G(i)>f(i). So how is the below formula valid?

According to this, f(i) will be greater than G(i).

Not at all, the mobius function takes negative values too, using the plus sign doesn’t mean we’re adding positive stuff, the G(i) +Z line, Z will be \le 0

https://codeforces.com/blog/entry/61053