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 :
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.