【精進2】ABC349-F

問題

F - Subsequence LCM
長さ$N$の正整数列$A$と正整数$M$が与えられる。$A$の空でない部分列であって、列の要素の最小公倍数が$M$になるようなものの個数を求めよ。

解法

$M$の約数でないような要素は選ばないので予め取り除いておく。また、空でないという条件を無視する($M=1$のときのみ空部分列を数えてしまうから、答えから$1$引く)。 $M = \prod_{i=1}^K p_i^{e_i}$と素因数分解する。ここで$M\le 10^{16}$であるから$K\le 12$である。部分列に対する条件は各$i \in[K]$について列の最小公倍数が$p_i^{e_i}$を約数として含むことになる。以前の問題同様、包除原理で求めていく。$i$について上記の条件を満たすような部分列全体の集合を$B_i$とすると求めるのは$|\bigcap_{i \in [K]}B_i| $である。包除原理から

\[\left| \bigcap_{i\in[K]}B_i \right| = \sum_{S \subseteq [K]}(-1)^{|S|}\left| \bigcap_{i\in S}\overline{B_i} \right|\]

と書ける。

これを \(\sum_{S \subseteq [K]}(-1)^{K-|S|}\left| \bigcap_{i\in [K]\setminus S}\overline{B_i} \right|\)と変形しておく。$\bigcap_{i\in [K]\setminus S}\overline{B_i}$は任意の$i\not\in S$について、$p_i^{e_i}\nmid \mathrm{lcm} $となるような部分列全体の集合である。 $\mathrm{cnt}[S]$を任意の$i \not\in S$について$p_i^{e_i}\nmid A_x $となる$x$の個数とすると、$\left|\bigcap_{i\in [K]\setminus S}\overline{B_i}\right|=2^{\mathrm{cnt}[S]}$となる。従って$\mathrm{cnt}$が求まれば、 $\mathcal{O}(2^K)$で計算できる。 $\mathrm{c}[S] $を$i\in S \Leftrightarrow p_i^{e_i}\mid A_x$となる$x$の個数とすると、$\mathrm{cnt}[S] = \sum_{T\subseteq S}\mathrm{c}[T] $である。$\mathrm{c}$は各$A_x$についてちょうど$1$つの$S$が対応するから計算することができ、$\mathrm{c}$から$\mathrm{cnt} $は高速ゼータ変換で$\mathcal{O}(K 2^K)$で計算できる。$K$が小さいので普通に$\mathcal{O}(3^K)$で計算してもよさそう。

#include "/home/y_midori/cp/library/template.hpp"
#include "math/pollard_rho.hpp"
#include <atcoder/modint>
using mint = atcoder::modint998244353;
void solve() {
    LL(n, m);
    auto fc = factor_count(m);
    ll k = fc.size();
    vl pc;
    for(auto [p, e] : fc) {
        pc.push_back(1);
        rep(i, e) pc.back() *= p;
    }
    vl c(1 << k);
    rep(i, n) {
        LL(ai);
        if(m % ai)
            continue;
        ll s = 0;
        rep(j, k) {
            if(ai % pc[j] == 0) {
                s |= 1 << j;
            }
        }
        c[s]++;
    }
    // 高速ゼータ変換 c[s] <- Σ[t⊆s]c[t]
    rep(j, k) {
        rep(i, 1 << k) {
            if((i >> j) & 1) {
                c[i] += c[i ^ (1 << j)];
            }
        }
    }
    mint ans = 0;
    rep(s, 1 << k) {
        mint tmp = mint(2).pow(c[s]);
        if((k - __popcount(s)) & 1)
            tmp *= -1;
        ans += tmp;
    }
    if(m == 1)
        ans--;
    print(ans.val());
}

int main() {
    ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    cout << std::setprecision(16);
    int t = 1;
    rep(_, t) {
        solve();
    }
}

別解

ユーザー解説通りにやってみる。 $ Z^{-1}(Z(X)\cdot Z(Y))_n= \sum _{\mathrm{lcm}(i,j)=n}X_i Y_j $ が重要。$(Z(X_1)\cdot Z(X_2)\cdot\dots \cdot Z(X_N))_n = 2^{c_n} $となる$c$をゼータ変換で求めて左辺の列をメビウス変換する。
(追記)ついでにlibrary checkerのgcd convolution,lcm convolutionを解いてライブラリに追加しておいた。

#include "/home/y_midori/cp/library/template.hpp"
#include "data_structure/hash-map-variable-length.hpp"
#include "math/pollard_rho.hpp"
#include <atcoder/modint>
using mint = atcoder::modint998244353;
void solve() {
    LL(n, m);
    vl d = divisors(m);
    HashMap<ll> c;
    rep(i, n) {
        LL(a);
        if(m % a)
            continue;
        c[a]++;
    }
    auto fc = factor_count(m);
    // ゼータ変換
    for(auto [p, e] : fc) {
        rep(i, d.size()) {
            if(d[i] % p == 0) {
                c[d[i]] += c[d[i] / p];
            }
        }
    }
    HashMap<mint> x;
    rep(i, d.size()) {
        x[d[i]] = mint(2).pow(c[d[i]]);
    }

    // メビウス変換
    for(auto [p, e] : factor_count(m)) {
        for(int i = ssize(d) - 1; i >= 0; i--) {
            if(d[i] % p == 0) {
                x[d[i]] -= x[d[i] / p];
            }
        }
    }
    auto ans = x[m];
    if(m == 1)
        ans--;
    cout << ans.val() << endl;
}

int main() {
    ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    cout << std::setprecision(16);
    int t = 1;
    rep(_, t) {
        solve();
    }
}