\(M\) の正の約数 \(d\) 全てについて最小公倍数が \(d\) になる部分列の個数を計算量 \(O(\sqrt{M}+(N+Kd(M))\log d(M))\) 程度で列挙することができます．ここで \(K\) は \(M\) の素因数の個数， \(d(M)\) は \(M\) の正の約数の個数です．

前計算として，\(M\) の約数および素因数の列挙をしておきます．共に \(O(\sqrt{M})\) でできます．

以下で出てくる数列の添字は基本的に \(M\) の正の約数のみを考えます．

• 数列 \(X\) のゼータ変換 \(Z(X)\) を \(Z(X)_n=\sum_{d\mid n}X_d\) で定めます．
• 数列 \(Y\) に対し \(Z(X)=Y\) をみたす \(X\) を \(Z^{-1}(Y)\) と表し，\(Y\) のメビウス変換と呼びます．
• 数列 \(X,Y\) の各点積 \(X\cdot Y\) を \((X\cdot Y)_n=X_nY_n\) で定めます．

ゼータ変換は次のようにして計算できます（\(\leftarrow\) は代入を表します）．

• \(M\) の各素因数 \(p\) に対し次を行う．
  • \(M\) の正の約数 \(d\) に対し昇順に次を行う．
    • \(X_{pd}\leftarrow X_{pd}+X_d\)

基本的に連想配列などで管理することになるので計算量は \(O(Kd(M)\log d(M))\) になるはずです．またメビウス変換も同様に計算できます．
（追記：このまま実装すると \(pd\) がオーバーフローする可能性があることに注意が必要です．また通常の配列で管理し尺取り法をすることで \(O(Kd(M))\) にすることもできます．）

さて，数列 \(X,Y\) に対し次が成り立ちます．

\[Z^{-1}(Z(X)\cdot Z(Y))_n=\sum_{\mathrm{lcm}(i,j)=n}X_iY_j\]

これを用いると，数列 \(X_i\) を \((X_i)_1=(X_i)_{A_i}=1,(X_i)_j=0\,(j\neq 1,A_i)\) で定めたとき次が求められればよいです（空部分列の考慮は必要になります）．

\[Z^{-1}(Z(X_1)\cdot Z(X_2)\cdot\cdots\cdot Z(X_N))\]

ここで \(Z(X_i)\) を具体的に考えると，\(Z(X_i)_n\) は \(A_i\mid n\) のとき \(2\)，それ以外のとき \(1\) になります．従って \(A_i\mid n\) となる \(i\) の個数を \(c\) として \((Z(X_1)\cdot\cdots\cdot Z(X_N))_n=2^c\) となるので，ゼータ変換を応用すれば \(O((N+Kd(N))\log d(N))\) で計算できます．

従って全体で \(O(\sqrt{M}+(N+Kd(N))\log d(N))\) です．

参考になりそうなページ

• 👊パンチ👊ゼータ変換👊パンチ👊 - ganmo::cout
• [数学] 畳み込み入門：Dirichlet積とゼータ変換・メビウス変換 - maspyのHP