PGCDS - EDITORIAL

PROBLEM LINK:

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

Author: Kimi Wu
Tester: Anshu Garg
Editorialist: Istvan Nagy

DIFFICULTY:

MEDIUM

PREREQUISITES:

Number theory, Inclusion-Exclusion principle, Combinatorics

PROBLEM:

We define the GCD set of an array A of length N as the set containing the GCD of all prefixes of A. More formally, the GCD set G of A contains the values \gcd([A_1]), \gcd([A_1, A_2]), \gcd([A_1, A_2, A_3]), \ldots, \gcd([A_1, A_2, \ldots, A_n]), where \gcd(B) is the largest integer d such that all the numbers in array B are divisible by d. Note that all values in G are distinct.

For example, \gcd([6]) = 6, \gcd([6, 12]) = 6, \gcd([6, 12, 9]) = 3, and \gcd([6, 12, 9, 2]) = 1. Therefore, the GCD set of the array [6, 12, 9, 2] is \{1, 3, 6\}.

You are given an integer set S of size M and two integers N, MOD, where MOD is a prime. You need to calculate the probability that the GCD set of P is equal to S, where P is a permutation of values 1 to N chosen uniformly randomly.

It can be shown that the answer is in the form of \frac{P}{Q} where P and Q are coprime non-negative integers and Q \neq 0. You need to find the value P \cdot Q^{-1} \pmod{MOD}.

EXPLANATION

During the calculation of the gcd-s in the following order: \gcd([A_1]), gcd([A_1, A_2]), \gcd([A_1, A_2, A_3]), \ldots, \gcd([A_1, A_2, \ldots, A_n]) we can observe few properties:

  • In every step the result cannot exceed the previous one, moreover
  • In every step the result is always a divisor of the previous one because gcd([a,b,c])=gcd(gcd(a,b),c)
  • At the last step gcd result is 1: because the 1 occurs in the array
  • A_1=S_M

We also can realize if a permutation generate S set then the first item of the permutation is S_M. It seems to calculate the next item of a good permutation we have to calculate how much integer exists in the range 1\leq x \leq N to get the second item of the set: gcd(S_M,x)=S_{M-1} or more generally:

Count the number of 1\leq x \leq N where gcd(a,x)=b.

This equivalent with the following equation: 1\leq x \leq N/b where gcd(a/b,x)=1. So firstly we should count a relative primes to a given number c in a range [1,K].

Its easy to see its enough to solve the special case when c is a product of different prime numbers, because if a number is relative prime to a prime number then it is also relative prime any power of that prime number.
First consider a simple case when c is a prime then the non-relative primes in the range : c, 2\cdot c,..,(K/c)\cdot c, so the number of solution is K-K/c. In the case c is product of two different prime numbers p_1\times p_2 then the number of non-relative numbers K-K/p_1 + K-K/p_2 but we subtract in both the following numbers: p_1\cdot p_2, 2\cdot p_1\cdot p_2,...,\frac{K}{c}\cdot p_1\cdot p_2 so we get the number of relative primes K-\frac{K}{p_1}-\frac{K}{p_2}+\frac{K}{p_1\cdot p_2}. We can solve it generally when c is a product of different prime factors with using the [inclusion-exclusion principle] (Inclusion–exclusion principle - Wikipedia). Because of the N is less than 1 billion, the number of different prime factors cannot exceed 9, so we can use the above principle it has at most 2^9 term.

Lets consider the following partition of the [1, N] range: {C_1, C_2, .., C_M} where :

  • C_M=\{S_M\}
  • C_{M-1}=\{x: gcd(x, S_{M})=S_{M-1}\}
  • …
  • C_i=\{x: gcd(x, S_{i+1})=S_{1}\}

If we want to create a permutation what generates the S set then we have to produce the gcd's in decreasing order so lets consider the case when the last generated gcd value is S_i, and we want to change gcd to S_{i-1}. In this case we have to choose a number from C_{i-1}.

  • If we choose a number from a larger index partition then gcd won’t change
  • If we choose a number from a smaller index partition then we get smaller gcd than S_{i-1}, after that we cannot generate the S set.

With these observations we can define the a permutation is good if and only if the first appearance of an item from C_i then all at least one item from the C_j where j>i already appeared, and no item occured from C_k where k<i.

Calculate the final probability from the C_i's : If we consider the N slots, firstly we have |C_M| option. The remaining integers from C_M we can assign any slots so we have (N-1)\cdot (N-2) \cdot ...\cdot(N-|C_M|) possibility. Now check the C_{M-1} set: the first item has to go to the first non used slot, we have |C_M-1| choice which item put there. The other items can assign all the non used slots so we have $(N-|C_M|-1)\cdot (N-|C_M|-2) \cdot …\cdot(N-|C_M|-|C_{M-1}|) options.

So we get the final number of correct permutation is : \frac{N!\cdot \prod_{i=1}^K C_i}{\prod_{i=1}^K(\sum_{j=1}^i C_j)}. Because of the task asks probability of the correct permutation we can get rid off the N! term :

\frac{\prod_{i=1}^K C_i}{\prod_{i=1}^K(\sum_{j=1}^i C_j)}

TIME COMPLEXITY:

O(sqrt(N)+M) per test case

SOLUTIONS:

Setter's Solution
#include <bits/stdc++.h>
using namespace std;

long long modpow(long long a, int b, int mod) {
    long long ans = 1;
    for (; b; b >>= 1, a = a * a % mod) {
        if (b & 1) {
            ans = ans * a % mod;
        }
    }
    return ans;
}

int get(int a, int b, int n) {
    // count the number 1 <= x <= n and gcd(a, x) = b
    n /= b, a /= b, b = 1;
    vector <int> p;
    int tmp = a;
    for (int i = 2; i * i <= tmp; ++i) if (tmp % i == 0) {
        p.push_back(i);
        while (tmp % i == 0) tmp /= i;
    }
    if (tmp > 1) p.push_back(tmp);
    int m = p.size();
    int ans = 0;
    for (int s = 0; s < 1 << m; ++s) {
        int res = 1;
        for (int i = 0; i < m; ++i) if (s >> i & 1) {
            res *= p[i];
        }
        ans += (n / res) * (__builtin_popcount(s) & 1 ? -1 : 1);
    }
    return ans;
}

void solve() {
    int n, m, mod;
    cin >> n >> m >> mod;
    vector <int> a(m);
    for (int i = 0; i < m; ++i) {
        cin >> a[i];
    }
    reverse(a.begin(), a.end());
    for (int i = 1; i < m; ++i) if (a[i - 1] % a[i] != 0) {
        cout << 0 << endl;
        return;
    }
    if (a[m - 1] != 1) {
        cout << 0 << endl;
        return;
    }
    long long ans = 1;
    for (int i = 1; i < m; ++i) {
        ans = ans * get(a[i - 1], a[i], n) % mod;
    }
    vector <int> cnt(m, 0);
    for (int i = 0; i < m; ++i) {
        cnt[i] = n / a[i];
        if (i) cnt[i] -= n / a[i - 1];
    }
    for (int i = 0; i < m; ++i) {
        cnt[i]--;
    }
    int pre = 0;
    for (int i = m - 1; ~i; --i) {
        pre += cnt[i] + 1;
        ans = ans * modpow(pre, mod - 2, mod) % mod;
    }
    cout << ans << endl;
}

int main () {
  ios::sync_with_stdio(false);
  cin.tie(0);
    int t;
    cin >> t;
    while (t--) {
        solve();
    }
}
Tester's Solution
// Permutation and GCD Set
#include<bits/stdc++.h>
using namespace std ;

#define ll              long long 
#define pb              push_back
#define all(v)          v.begin(),v.end()
#define sz(a)           (ll)a.size()
#define F               first
#define S               second
#define INF             2000000000000000000
#define popcount(x)     __builtin_popcountll(x)
#define pll             pair<ll,ll>
#define pii             pair<int,int>
#define ld              long double

const int M = 1000000007;
const int MM = 998244353;

template<typename T, typename U> static inline void amin(T &x, U y){ if(y<x) x=y; }
template<typename T, typename U> static inline void amax(T &x, U y){ if(x<y) x=y; }

#ifdef LOCAL
#define debug(...) debug_out(#__VA_ARGS__, __VA_ARGS__)
#else
#define debug(...) 2351
#endif

long long readInt(long long l,long long r,char end){
    long long x = 0;
    int cnt = 0;
    int first =-1;
    bool is_neg = false;
    while(true) {
        char g = getchar();
        if(g == '-') {
            assert(first == -1);
            is_neg = true;
            continue;
        }
        if('0' <= g && g <= '9') {
            x *= 10;
            x += g - '0';
            if(cnt == 0) {
                first = g - '0';
            }
            ++cnt;
            assert(first != 0 || cnt == 1);
            assert(first != 0 || is_neg == false);
            
            assert(!(cnt > 19 || (cnt == 19 && first > 1)));
        } 
        else if(g == end) {
            if(is_neg) {
                x = -x;
            }
            assert(l <= x && x <= r);
            return x;
        } 
        else {
            assert(false);
        }
    }
}
string readString(int l,int r,char end){
    string ret = "";
    int cnt = 0;
    while(true) {
        char g = getchar();
        assert(g != -1);
        if(g == end) {
            break;
        }
        ++cnt;
        ret += g;
    }
    assert(l <= cnt && cnt <= r);
    return ret;
}
long long readIntSp(long long l,long long r){
    return readInt(l,r,' ');
}
long long readIntLn(long long l,long long r){
    return readInt(l,r,'\n');
}
string readStringLn(int l,int r){
    return readString(l,r,'\n');
}
string readStringSp(int l,int r){
    return readString(l,r,' ');
}

    
int MOD=1000000007;
struct Mint {
    int val;
 
    Mint(long long v = 0) {
        if (v < 0)
            v = v % MOD + MOD;
        if (v >= MOD)
            v %= MOD;
        val = v;
    }
 
    static int mod_inv(int a, int m = MOD) {
        int g = m, r = a, x = 0, y = 1;
        while (r != 0) {
            int q = g / r;
            g %= r; swap(g, r);
            x -= q * y; swap(x, y);
        } 
        return x < 0 ? x + m : x;
    } 
    explicit operator int() const {
        return val;
    }
    Mint& operator+=(const Mint &other) {
        val += other.val;
        if (val >= MOD) val -= MOD;
        return *this;
    }
    Mint& operator-=(const Mint &other) {
        val -= other.val;
        if (val < 0) val += MOD;
        return *this;
    }
    static unsigned fast_mod(uint64_t x, unsigned m = MOD) {
           #if !defined(_WIN32) || defined(_WIN64)
                return x % m;
           #endif
           unsigned x_high = x >> 32, x_low = (unsigned) x;
           unsigned quot, rem;
           asm("divl %4\n"
            : "=a" (quot), "=d" (rem)
            : "d" (x_high), "a" (x_low), "r" (m));
           return rem;
    }
    Mint& operator*=(const Mint &other) {
        val = fast_mod((uint64_t) val * other.val);
        return *this;
    }
    Mint& operator/=(const Mint &other) {
        return *this *= other.inv();
    }
    friend Mint operator+(const Mint &a, const Mint &b) { return Mint(a) += b; }
    friend Mint operator-(const Mint &a, const Mint &b) { return Mint(a) -= b; }
    friend Mint operator*(const Mint &a, const Mint &b) { return Mint(a) *= b; }
    friend Mint operator/(const Mint &a, const Mint &b) { return Mint(a) /= b; }
    Mint& operator++() {
        val = val == MOD - 1 ? 0 : val + 1;
        return *this;
    }
    Mint& operator--() {
        val = val == 0 ? MOD - 1 : val - 1;
        return *this;
    }
    Mint operator++(int32_t) { Mint before = *this; ++*this; return before; }
    Mint operator--(int32_t) { Mint before = *this; --*this; return before; }
    Mint operator-() const {
        return val == 0 ? 0 : MOD - val;
    }
    bool operator==(const Mint &other) const { return val == other.val; }
    bool operator!=(const Mint &other) const { return val != other.val; }
    Mint inv() const {
        return mod_inv(val);
    }
    Mint power(long long p) const {
        assert(p >= 0);
        Mint a = *this, result = 1;
        while (p > 0) {
            if (p & 1)
                result *= a;
 
            a *= a;
            p >>= 1;
        }
        return result;
    }
    friend ostream& operator << (ostream &stream, const Mint &m) {
        return stream << m.val;
    }
    friend istream& operator >> (istream &stream, Mint &m) {
        return stream>>m.val;   
    }
};

bool prime(int n) {
    for(int i=2;i*i<=n;++i) {
        if(n % i == 0) {
            return false;
        }
    }
    return true;
}
int sumM = 0;

int _runtimeTerror_()
{
    int n, m;
    cin >> n >> m >> MOD;
    // n = readIntSp(1,1e9),m = readIntSp(1,2e4), MOD = readIntLn(1e9+1,2e9-1);
    vector<int> a(m);
    sumM += m;
    for(int i=0;i<m;++i) {
        cin >> a[i];
        continue;
        if(i < m - 1) {
            a[i] = readIntSp(1,n);
        }
        else {
            a[i] = readIntLn(1,n);
        }
    }
    assert(prime(MOD));
    auto get = [&](int x,int up) {
        vector<int> pr;
        for(int i=2;i*i<=x;++i) {
            if(x % i == 0) {
                while(x % i == 0) {
                    x /= i;
                }
                pr.push_back(i);
            }
        }
        if(x > 1) {
            pr.push_back(x);
        }
        long long ans = 0;
        int n = sz(pr);
        for(int i=0;i<(1 << n);++i) {
            int product = 1;
            for(int j=0;j<n;++j) {
                if(i & (1 << j)) {
                    product *= pr[j];
                }
            }
            ans += popcount(i) & 1 ? -up/product : up/product;
        }
        return ans;
    };
    assert(is_sorted(all(a)));
    reverse(all(a));
    Mint ans = 1;
    int cnt = 0;
    for(int i=0;i<m;++i) {
        if(i == 0) {
            ans /= n;
        }
        else {
            assert(a[i] != a[i-1]);
            if(a[i-1] % a[i]) {
                cout << "0\n";
                return 0;
            }
            auto u = get(a[i-1]/a[i], n/a[i]);
            ans *= u;
            ans /= n - (n/a[i-1]);
        }
    }   
    if(a.back() > 1) {
        ans = 0;
    }
    debug(n,a);
    cerr << "Prime: " << MOD << "\n";
    cout << ans << "\n";
    return 0;
}

int main()
{
    // freopen("PGCDS-0.in","r",stdin);
    ios_base::sync_with_stdio(0);cin.tie(0);cout.tie(0);
    #ifdef runSieve
        sieve();
    #endif
    #ifdef NCR
        initialize();
    #endif
    int TESTS = 1;
    cin >> TESTS;
    // TESTS = readIntLn(1,10);
    while(TESTS--)
        _runtimeTerror_();
    assert(getchar() == -1);
    cerr << "sumM: " << sumM << "\n";
    return 0;
}

If you have other approaches or solutions, let’s discuss in comments.

4 Likes

" With these observations we can define the a permutation is good if and only if the first appearance of an item from C_i then all at least one item from the C_j where j>i already appeared, and no item occured from C_k where k<i. "

How to do this math? From these sets how do we find the correct answer?

Sorry I just now finished that part of the editorial.

Wait! Wrong links have been posted in the editorial.

Editorial for Contest Page | CodeChef is already published. This is editorial of which problem?

Thanks, fixed.