\(M=1\) のときは \(A_i=1\) となる \(i\) の個数を \(c\) として答えは \(2^c-1\) です．そうでない場合を考えます．

\(A_i\) の要素のうち \(M\) の約数のもののみを考えればよいです．以下では，すべての \(A\) の要素が \(M\) の約数であると仮定します．

\(M\) の素因数分解を \(p_1^{e_1}p_2^{e_2}\dots p_K^{e_K}\) とします．また，\(M\) の約数 \(x\) に対して長さ \(K\) の bit 列 \(f(x)\) を以下のように定めます．

• 各 \(i=1,\dots,K\) に対して，\(x\) が \(p_i^{e_i}\) の倍数ならば \(i\) bit 目は \(1\) ，そうでないならば \(i\) bit 目は \(0\) である．

（例えば \(M=12=2^23^1\) とき，\((f(1),f(2),f(3),f(4),f(6),f(12))=(00,00,10,01,10,11)\) です．\(6\) は \(2^2\) の倍数ではなく \(3^1\) の倍数であるため \(2\) bit 目のみが立ちます．）

このとき \(A\) の部分列 \(B\) について，以下の \(5\) つは同値です．

• \(B\) の最小公倍数は \(M\) である．
• 各 \(i=1,\ldots,K\) に対し，\(B\) の最小公倍数は \(p_i^{e_i}\) の倍数である．
• 各 \(i=1,\ldots,K\) に対し，\(B\) は \(p_i^{e_i}\) の倍数であるような要素を少なくとも一つ含む．
• 各 \(i=1,\ldots,K\) に対し，\(f(B_1),f(B_2),\dots,f(B_{|B|})\) の bitwise or の \(i\) bit 目は \(1\) である．
• \(f(B_1),f(B_2),\dots,f(B_{|B|})\) の bitwise or は \(2^K-1(=\overbrace{11\ldots11}^{K\ \mathrm{bit}})\) である．

\(1\) つ目と \(5\) つ目が同値であることより，次のような問題に帰着することができます．

\(0\) 以上 \(2^K\) 未満の整数からなる長さ \(N\) の整数列 \(C=(f(A_1),f(A_2),\ldots,f(A_N))\) がある．\(C\) の部分列であって，各要素の bitwise or が \(2^K-1\) になるものはいくつあるか？

この問題の解法を \(2\) 通り紹介します．

1. dp

\(C\) に含まれる \(i\) の個数を \(\mathrm{cnt}(i)\) とします．\(i=0,1,\ldots,2^K-1\) の順に，部分列に \(i\) を一つ以上含めるかどうかを決めていきます．\(\mathrm{dp}[i][j ]\) で \(i\) まで見たときに bitwise or が \(j\) であるような部分列の個数とします．このとき， \(i\) を部分列に含む場合は \(\mathrm{dp}[i][i\ \mathrm{or}\ j]\leftarrow \mathrm{dp}[i-1][j]\times (2^{\mathrm{cnt}(i)}-1)\)，含まない場合は \(\mathrm{dp}[i][j]\leftarrow \mathrm{dp}[i-1][j]\) という遷移をたどります．この dp は状態数，遷移ともに \(O(2^K)\) なので計算量は \(O(4^K)\) です．

2. ゼータ・メビウス変換

\(C\) に含まれる \(i\) の個数を \(\mathrm{cnt}(i)\), 各要素の bitwise or が \(i\) になるような \(C\) の部分列の個数を \(g(i)\) とします．求めたいものは \(g(2^K-1)\) です．\(g\) のゼータ変換を \(h\) とすると \(h(i)=\sum_{j\subset i}g(j)=2^{\sum_{j\subset i}\mathrm{cnt}(i)}\) が成り立ちます．

（例えば \(K=3\) のとき \(h(5)=g(0)+g(1)+g(4)+g(5)\) は bitwise or が \(0,1,4,5\) のいずれかであるような部分列の個数を表しています．これは \(C\) の要素のうち \(0,1,4,5\) であるものからいくつか選ぶ方法の通り数と等しいので \(2^{\mathrm{cnt}(0)+\mathrm{cnt}(1)+\mathrm{cnt}(4)+\mathrm{cnt}(5)}\) です．）

\(\mathrm{cnt}\) のゼータ変換によって \(h\) を求めることができ，\(h\) から \(g\) へはメビウス変換（ゼータ変換の逆）によって得ることができます．（\(h\) から \(g(2^K-1)\) だけを求めるのであれば包除原理で \(O(2^K)\) でも求められます．）ゼータ変換やメビウス変換は \(O(K2^K)\) などの計算量で行うことができます．

素因数分解は \(O(\sqrt{M})\)，各 \(i\) に対して \(f(A_i)\) は \(O(K)\) で求められるので，この問題を \(O(\sqrt{M}+NK+K2^K)\) や \(O(\sqrt{M}+NK+4^K)\) で解くことができます．\(M\leq 10^{16}\) の制約の下では \(K\leq 13\) が成り立つので適切に実装をすれば間に合います．

実装例 (Python)

• dp
• ゼータ・メビウス変換