Codeforces CodeTON Round 3 D. Count GCD【状压、容斥原理计数】

D. Count GCD

题意

给定一个长度为 $n$ 的正整数数组 $a$，$\forall 1\le i\le n，a_i\le m$

先要要统计满足以下条件的数组 $b$ 的数量：

• $\forall 1\le i\le n，b_i\le m$
• $gcd(b_1,b_2,b_3,...b_i)=a_i$

思路

可以发现在当前的 $i$，$gcd(a_1,a_2,...a_{i-1})=a_{i-1}$，那么只要 $gcd(a_{i-1},b_i)=a_i$ 就可以了，那么我们可以得出结论：$a_i\mid a_{i-1}$，也即是 $a_{i-1}$ 是 $a_i$ 的因子，同时 $a_i\mid b_i$

这样子答案就是每个位置可选的数量累乘，问题是如何快速求出位置 $i$ 可行的 $b_i$ 数量

$b_i$ 一定得是 $a_i$ 的倍数，我们可以假设 $k\cdot a_i=b_i$，那么不同的 $k$ 的数量就是 $b_i$ 的数量
那么 $gcd(a_{i-1},b_i)=a_i\Rightarrow gcd(\frac{a_{i-1}}{a_i},k)=1$，也就是 $\frac{a_{i-1}}{a_i}$ 中出现的质因子不能在 $k$ 中出现

为了快速得到 $\frac{a_{i-1}}{a_i}$ 的所有质因子，我们发现 $\frac{a_{i-1}}{a_i}$ 的质因子一定是 $a_1$ 的质因子。这是因为 $a_i\mid a_{i-1}\mid ....\mid a_2\mid a_1$；并且最小的 $10$ 个质数相乘就已经大于 $10^9$ 了！所以我们只需要最多 $9$ 个不同的质数（不一定是最小的 $9$ 个）就可以表示小于等于 $10^9$ 的数了

我们只需要在 $O(m)$ 下预处理 得出$a_1$ 的质因子，然后对于每个 $i>2$，枚举 $a_1$ 的质因子，看看 $a_i$ 是否被其整除，即可快速得出 $a_i$ 的所有质因子，质因子数量都不超过 $9$ 个！我们这里就可以考虑状压 + 容斥原理来计数了

容斥原理：
集合 $S$ 的子集 $A_1$ 具有性质 $P_1$，子集 $A_2$ 具有性质 $P_2$，… $A_n$ 具有性质 $P_n$


（1）集合 $S$ 中不具有性质 $P_1,P_2,...,P_n$ 的对象个数为：
$|\overline{A_1}\cap \overline{A_2},...,\cap \overline{A_n}|=|S|-\sum |A_i|+\sum |A_i\cap A_j|-\sum |A_i\cap A_j\cap A_k|+...+(-1{)}^n|A_1\cap A_2\cap ...\cap A_n|$


（2）集合 $S$ 中至少具有性质 $P_1,P_2,...P_n$ 之一的对象个数为：
$|A_1\cup A_2,...,\cup A_n|=\sum |A_i|-\sum |A_i\cap A_j|+\sum |A_i\cap A_j\cap A_k|+...+(-1{)}^{n+1}|A_1\cap A_2\cap ...\cap A_n|$

因为性质很少，最多只有 $9$ 个（分别被不同的 $9$ 个质数整除），我们要求出与 $\frac{a_{i-1}}{a_i}$ 互质的 $k$ 的数量，也就是 $k$ 不能拥有 $\frac{a_{i-1}}{a_i}$ 的质因子，也就是 $k$ 不具有性质 $P_1,P_2,...P_k$（$k$ 为 $\frac{a_{i-1}}{a_i}$ 质因子数），只需要用容斥原理简单地通过整除得到即可，注意 $k\le \frac{m}{a_i}$

时间复杂度：$O(\sqrt{m}+9\cdot 2^9\cdot \log m)$

最后那个 $\log m$ 是因为：$a_1$ 的每个质因子数量一定非递减，如果某个质因子在 $\frac{a_{i-1}}{a_i}$ 中出现，说明这个质因子数量在这里至少 $-1$ 了，那么从 $a_1$ 减少到 $a_n$，质因子最多就 $\log m$ 个，很快就没了

#include<bits/stdc++.h>
#define fore(i,l,r)	for(int i=(int)(l);i<(int)(r);++i)
#define fi first
#define se second
#define endl '\n'
#define ull unsigned long long
#define ALL(v) v.begin(), v.end()
#define Debug(x, ed) std::cerr << #x << " = " << x << ed;

const int INF=0x3f3f3f3f;
const long long INFLL=1e18;

typedef long long ll;

template<class T>
constexpr T power(T a, ll b){
    T res = 1;
    while(b){
        if(b&1) res = res * a;
        a = a * a;
        b >>= 1;
    }
    return res;
}

constexpr ll mul(ll a,ll b,ll mod){ //快速乘，避免两个long long相乘取模溢出
    ll res = a * b - ll(1.L * a * b / mod) * mod;
    res %= mod;
    if(res < 0) res += mod; //误差
    return res;
}

template<ll P>
struct MLL{
    ll x;
    constexpr MLL() = default;
    constexpr MLL(ll x) : x(norm(x % getMod())) {}

    static ll Mod;
    constexpr static ll getMod(){
       if(P > 0) return P;
       return Mod;
    }

    constexpr static void setMod(int _Mod){
       Mod = _Mod;
    }
    constexpr ll norm(ll x) const{
       if(x < 0){
           x += getMod();
       }
       if(x >= getMod()){
           x -= getMod();
       }
       return x;
    }
    constexpr ll val() const{
       return x;
    }
    explicit constexpr operator ll() const{ 
       return x; //将结构体显示转换为ll类型： ll res = static_cast<ll>(OBJ)
    }
    constexpr MLL operator -() const{ //负号，等价于加上Mod
       MLL res;
       res.x = norm(getMod() - x);
       return res;
    }
    constexpr MLL inv() const{
       assert(x != 0);
       return power(*this, getMod() - 2); //用费马小定理求逆
    }
    constexpr MLL& operator *= (MLL rhs) & { //& 表示“this”指针不能指向一个临时对象或const对象
       x = mul(x, rhs.x, getMod()); //该函数只能被一个左值调用
       return *this;
    }
    constexpr MLL& operator += (MLL rhs) & {
       x = norm(x + rhs.x);
       return *this;
    }
    constexpr MLL& operator -= (MLL rhs) & {
       x = norm(x - rhs.x);
       return *this;
    }
    constexpr MLL& operator /= (MLL rhs) & {
       return *this *= rhs.inv();
    }
    friend constexpr MLL operator * (MLL lhs, MLL rhs){
       MLL res = lhs;
       res *= rhs;
       return res;
    }
    friend constexpr MLL operator + (MLL lhs, MLL rhs){
       MLL res = lhs;
       res += rhs;
       return res;
    }
    friend constexpr MLL operator - (MLL lhs, MLL rhs){
       MLL res = lhs;
       res -= rhs;
       return res;
    }
    friend constexpr MLL operator / (MLL lhs, MLL rhs){
       MLL res = lhs;
       res /= rhs;
       return res;
    }
    friend constexpr std::istream& operator >> (std::istream& is, MLL& a){
       ll v;
       is >> v;
       a = MLL(v);
       return is;
    }
    friend constexpr std::ostream& operator << (std::ostream& os, MLL& a){
       return os << a.val();
    }
    friend constexpr bool operator == (MLL lhs, MLL rhs){
       return lhs.val() == rhs.val();
    }
    friend constexpr bool operator != (MLL lhs, MLL rhs){
       return lhs.val() != rhs.val();
    }
};

const ll mod = 998244353;
using Z = MLL<mod>;

std::vector<int> get_d(int x){
    std::vector<int> d;
    for(int i = 2; i * i <= x; ++i){
        if(x % i ) continue;
        d.push_back(i);
        while(x % i == 0) x /= i;
    }
    if(x > 1) d.push_back(x);
    return d;
}

void solve(){
    int n, m;
    std::cin >> n >> m;
    std::vector<int> a(n + 1);
    fore(i, 1, n + 1) std::cin >> a[i];
    fore(i, 2, n + 1)
        if(a[i - 1] % a[i]){
            std::cout << "0\n";
            return;
        }

    auto div = get_d(a[1]);
    Z ans = 1;
    fore(i, 2, n + 1){
        int val = a[i - 1] / a[i];
        std::vector<int> left;
        for(auto d : div)
            if(val % d == 0)
                left.push_back(d);

        int sz = left.size();
        ll K = m / a[i]; //k的上界
        ll tot = K;
        fore(mask, 1, 1 << sz){
            int cnt = __builtin_popcount(mask); //容斥原理系数
            int now = 1;
            fore(j, 0, sz)
                if(mask >> j & 1)
                    now *= left[j];
            if(cnt & 1) tot -= K / now;
            else tot += K / now;
        }
        
        ans *= tot;
    }
    std::cout << ans << endl;
}

int main(){
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    std::cout.tie(nullptr);
    int t;
    std::cin >> t;
    while(t--){
        solve();
    }
    return 0;
}