If \(M=1\), then the answer is \(2^c-1\), where \(c\) is the number of \(i\) such that \(A_i=1\). We assume the other cases.

We can safely ignore the elements \(A_i\) that are not divisors of \(M\). We now assume that every element of \(A\) is a divisor of \(M\).

Suppose that \(M\) factorizes into \(p_1^{e_1}p_2^{e_2}\dots p_K^{e_K}\). For a divisor \(x\) of \(M\), define a length-\(K\) bit sequence \(f(x)\) as follows:

• for each \(i=1,\dots,K\), the \(i\)-th bit is \(1\) if \(x\) is a multiple of \(p_i^{e_i}\), and \(0\) otherwise.

(For example, if \(M=12=2^23^1\), then \((f(1),f(2),f(3),f(4),f(6),f(12))=(00,00,10,01,10,11)\). Since \(6\) is not a multiple of \(2^2\) but of \(3^1\), only the \(2\)-nd bit is set.)

Then, for a subsequence \(B\) of \(A\), the following five are equivalent:

• The Least Common Multiple (LCM) of \(B\) is \(M\).
• For all \(i=1,\ldots,K\), the LCM of \(B\) is a multiple of \(p_i^{e_i}\).
• For all \(i=1,\ldots,K\), \(B\) contains at least one element that is a multiple of \(p_i^{e_i}\).
• For all \(i=1,\ldots,K\), the \(i\)-th bit of the bitwise OR of \(f(B_1),f(B_2),\dots,f(B_{|B|})\) is \(1\).
• The bitwise OR of \(f(B_1),f(B_2),\dots,f(B_{|B|})\) is \(2^K-1(=\overbrace{11\ldots11}^{K\ \mathrm{bit}})\).

Since the first and fifth propositions are equivalent, the problem is boiled down to the following:

You are given a length-\(N\) sequence \(C=(f(A_1),f(A_2),\ldots,f(A_N))\) consisting of integers between \(0\) (inclusive) and \(2^K\) (exclusive). How many subsequences of \(C\) have a bitwise OR of \(2^K-1\)?

We introduce two ways to solve this problem.

1. DP (Dynamic Programming)

Let \(\mathrm{cnt}(i)\) be the number of occurrences of \(i\) in \(C\). For \(i=0,1,\ldots,2^K-1\), we determine if the subsequence should contain at least one \(i\). Let \(\mathrm{dp}[i][j ]\) be the number of subsequences of the first \(i\) elements whose bitwise OR is \(j\). Then, if \(i\) is contained in the subsequence, it follows the transition \(\mathrm{dp}[i][i\ \mathrm{or}\ j]\leftarrow \mathrm{dp}[i-1][j]\times (2^{\mathrm{cnt}(i)}-1)\); if it is not, it follows \(\mathrm{dp}[i][j]\leftarrow \mathrm{dp}[i-1][j]\). This DP has \(O(2^K)\) states and each transition costs \(O(2^K)\) time, so the complexity is \(O(4^K)\).

2. Zeta and Möbius Transforms

Let \(\mathrm{cnt}(i)\) be the number of occurrences of \(i\) in \(C\), and \(g(i)\) be the number of subsequences of \(C\) whose bitwise OR is \(i\). What we want is \(g(2^K-1)\). If we denote by \(h\) the Zeta transform of \(g\), then \(h(i)=\sum_{j\subset i}g(j)=2^{\sum_{j\subset i}\mathrm{cnt}(i)}\).

(For example, if \(K=3\), the count \(h(5)=g(0)+g(1)+g(4)+g(5)\) represents the number of subsequences whose bitwise OR is \(0\), \(1\), \(4\), or \(5\). This is equivalent to the number of ways to choose elements of \(C\) that are equal to one of \(0\), \(1\), \(4\), and \(5\), so it is \(2^{\mathrm{cnt}(0)+\mathrm{cnt}(1)+\mathrm{cnt}(4)+\mathrm{cnt}(5)}\).)

\(h\) can be found as the Zeta transform of \(\mathrm{cnt}\), and \(g\) as the Möbius transform (inverse of Zeta transform) of \(h\). (One can also apply the inclusion-exclusion principle to obtain \(g(2^K-1)\) from \(h\) in \(O(2^K)\) time.) Zeta and Möbius transforms can be done in like \(O(K2^K)\) time.

Prime factorization costs \(O(\sqrt{M})\) time, and \(f(A_i)\) can be evaluated in \(O(K)\) time for each \(i\), so this problem can be solved in a total of \(O(\sqrt{M}+NK+K2^K)\) or \(O(\sqrt{M}+NK+4^K)\) time. Under the constraints of \(M\leq 10^{16}\), we have \(K\leq 13\), so it fits in the time limit with appropriate implementation.

Sample code (Python)

• DP
• Zeta and Möbius transforms