素因数分解を使わない解法です。

まず \(A_i\) が \(M\) の約数でなければ部分列に含めてはいけないので、そのような \(A_i\) は除いておきます。\(M\) の素因数分解を \(M = \prod_{j = 1} ^ K p_j ^ {e_j} \) とします。 今興味があるのはLCMを取った結果 \(M\) になるかどうかなので、\(A_i\) が素因数に \(p_j\) を \(e_j\) 個持たなければ、その素因数は \(0\) 個と見なしてしまってよいです。 つまり、\(A_i\) の代わりに以下の \(A'_i\) で置き換えても答えは変わりません。

\[ A'_i = \prod_{j=1}^{K} p_j^{e'_j}, \text{ where } e'_j = \begin{cases} e_j & \text{if } p_j^{e_j} \mid A_i \\ 0 & \text{otherwise} \end{cases} \]

この変換は以下のように素因数分解せずに行えます。

long long remove_redundant_primes(long long a, long long m) {
  assert(m % a == 0);
  // もし a が素因数 p_j を e_j 個持っていなければ m/a は素因数 p_j を持つ
  long long g = gcd(m / a, a);
  // g に含まれる p_j の個数に着目すると、以下の while は O(max_j log e_j) 回しか回らない
  while (g != 1) {
    a /= g;
    g = gcd(m / a, a);
  }
  return a;
}

これにより \(A_i\) は \(p_j ^ {e_j}\) をいくつか選んで掛け合わせたものだけになります。 素因数分解が行えればこのような \(p_j ^ {e_j}\) を列挙することで本解と同様に解くことができます。 ここでは素因数分解の代わりに以下を満たす \(f = (f_1, f_2, \ldots, f_{K'})\) を考えることにします。

• \(f_i\) は \(2\) 以上の整数 (\(1 \le i \le {K'}\))
• \(f_i\) と \(f_j\) は互いに素 (\(1 \le i \lt j \le {K'}\))
• 各 \(A_i\) について、ある \(S \subseteq \{1, 2, \ldots, {K'}\}\) が存在して \(A_i = \prod_{j \in S} f_j\) と表せる

このような \(f\) を \(1\) つ求めることができれば、LCMの分布を高速ゼータ変換・高速メビウス変換で求めることできます。 これを満たす \(f\) は以下のように\(A_i\) を一つ一つ見ていくと同時に、必要があればさらに細分化したり要素を追加したりすることで構築することが可能です。

vector<long long> compute_factors(const vector<long long>& a) {
  vector<long long> f;
  for (auto x : a) {
    for (auto fi : f) {
      if (x % fi == 0) {
        x /= fi;
      }
    }
    // 現在のfでxが表現できるようになっていればこの時点でx==1
    if (x != 1) {  // x!=1となるのは高々Mの素因数の種類数回
      for (int i = ssize(f) - 1; i >= 0; i--) {
        long long g = gcd(f[i], x);
        if (g != 1) {
          x /= g;
          // f[i]をgとf[i]/gに細分化
          f.push_back(g);
          f[i] /= g;
        }
      }
      if (x != 1) f.push_back(x);
    }
  }
  return f;
}

まとめると以下のようなコードで解くことができます。計算量は \(M\) の素因数の種類数を \(K\) 、\(M\) の約数の個数を \(d(M)\) として \( O(\min(N,d(M)) \log M \log \log M + K 2^K)\) です。

#include <bits/stdc++.h>
#include <atcoder/modint>
using namespace std;
using mint = atcoder::modint998244353;

long long remove_redundant_primes(long long a, long long m) {
  assert(m % a == 0);
  // もし a が素因数 p_j を e_j 個持っていなければ m/a は素因数 p_j を持つ
  long long g = gcd(m / a, a);
  // g に含まれる p_j の個数に着目すると、以下の while は O(max_j log e_j) 回しか回らない
  while (g != 1) {
    a /= g;
    g = gcd(m / a, a);
  }
  return a;
}

auto get_freq(const auto& a, long long m) {
  // freq_divs : mの約数の頻度
  // freq_reduced : 不要な素因数を削除した上での頻度
  unordered_map<long long, int> freq_divs, freq_reduced;
  for (auto x : a) if (m % x == 0) freq_divs[x]++;
  for (auto [k, v] : freq_divs) freq_reduced[remove_redundant_primes(k, m)] += v;
  return freq_reduced;
}

vector<long long> compute_factors(const vector<long long>& a) {
  vector<long long> f;
  for (auto x : a) {
    for (auto fi : f) {
      if (x % fi == 0) {
        x /= fi;
      }
    }
    // 現在のfでxが表現できるようになっていればこの時点でx==1
    if (x != 1) {  // x!=1となるのは高々Mの素因数の種類数回
      for (int i = ssize(f) - 1; i >= 0; i--) {
        long long g = gcd(f[i], x);
        if (g != 1) {
          x /= g;
          // f[i]をgとf[i]/gに細分化
          f.push_back(g);
          f[i] /= g;
        }
      }
      if (x != 1) f.push_back(x);
    }
  }
  return f;
}

vector<mint> lcm_dist(vector<int> a) {
  int n = a.size();
  int k = countr_zero(1U * n);
  for (int i = 0; i < k; i++) {
    for (int s = 0; s < n; s++) {
      if (~s >> i & 1) a[s | 1 << i] += a[s];
    }
  }
  vector<mint> d(n);
  for (int s = 0; s < n; s++) d[s] = mint(2).pow(a[s]) - 1;
  for (int i = 0; i < k; i++) {
    for (int s = 0; s < n; s++) {
      if (~s >> i & 1) d[s | 1 << i] -= d[s];
    }
  }
  return d;
}

int main() {
  int n;
  long long m;
  cin >> n >> m;
  vector<long long> a(n);
  for (auto& x : a) cin >> x;
  // 不要な要素・素因数を取り除いた上で、それらの頻度を求める
  unordered_map<long long, int> freq = get_freq(a, m);
  // fを求める
  vector<long long> b;
  for (auto [k, v] : freq) b.push_back(k);
  vector<long long> f = compute_factors(b);
  {
    // もしfの総積にMなっていなければ全部選んでも不可
    long long prod = 1;
    for (auto fi : f) prod *= fi;
    if (prod != m) {
      cout << 0 << endl;
      return 0;
    }
  }
  // 2^|f|通りの値がそれぞれいくつあるかを表す配列に直す
  int sz = f.size();
  vector<int> freq_arr(1 << sz);
  for (int s = 0; s < 1 << sz; s++) {
    long long k = 1;
    for (int i = 0; i < sz; i++) {
      if (s >> i & 1) k *= f[i];
    }
    if (auto it = freq.find(k); it != freq.end()) {
      freq_arr[s] = it->second;
    }
  }
  // LCMの分布を高速ゼータ変換・高速メビウス変換を用いて求める
  auto d = lcm_dist(move(freq_arr));
  mint ans = d.back();
  cout << ans.val() << endl;
}

提出(69 ms)