実装方針

\(A_i\) の倍数であるという条件を「条件 \(i\)」と呼ぶことにします．

\(N\) 個の条件からなる全体集合に対して部分集合 \(S \in \{1, \dots, N\}\) を考えます．

各部分集合 \(S\) について，\(S\) に属する条件は満たし，かつ他の条件は満たさないという条件を「一致条件 \(S\)」と呼ぶことにします．

このとき，すべての部分集合についてその一致条件を満たす \(1\) 以上 \(Y\) 以下の整数の個数が求まれば，\(|S| = M\) である部分集合 \(S\) すべてについてその合計を求めることで答えが求まります．

ここで，部分集合 \(T\) について，\(T\) に属する条件は満たし，他の条件は満たしても満たさなくてもよいという条件を「部分条件 \(T\)」と呼ぶことにします．

部分条件は一致条件よりも考えやすく，部分条件 \(\{i_1, \dots, i_k\}\) は「\(A_{i_1}, \dots, A_{i_k}\) の最小公倍数の倍数である」と簡潔に言い換えられます．この最小公倍数を \(x\) とおくと，この部分条件を満たす \(1\) 以上 \(Y\) 以下の整数の個数は \(\lfloor \frac{Y}{x} \rfloor\) で得られます．最小公倍数を求める際，\(Y\) を上回るならば \(Y+1\) として扱う，できる限り乗算ではなく除算を使って判定・計算するといった工夫をすることで，オーバーフローの危険を避けられます．

各部分条件を満たすものの個数が求まったら，これをゼータ変換の逆変換 (メビウス変換といいます) にかけることによって各一致条件を満たすものの個数が求まり，よって答えが求まります．これは，簡潔にいえば「条件 \(1, \dots, k\) については一致条件，条件 \(k+1, \dots, N\) については部分条件」というような条件の持ち方をし，部分条件から \(1\) つずつ一致条件に移していく形で計算することができます．

実装例 (C++)

#include <iostream>
using std::cin;
using std::cout;

typedef long long int ll;


ll n, m, y;
ll a[21];
ll zeta[1LL << 21];

ll mygcd (ll l, ll r) {
	if (r == 0) return l;
	return mygcd(r, l % r);
}

int main (void) {
	cin >> n >> m >> y;
	for (ll i = 0; i < n; i++) {
		cin >> a[i];
	}

	for (ll b = 0; !(b >> n); b++) {
		ll prod = 1;
		for (ll j = 0; j < n; j++) {
			if (b & (1LL << j)) {
				ll g = mygcd(prod, a[j]);
				if ((prod / g) > y / a[j]) {
					prod = y+1;
					break;
				} else {
					prod = (prod / g) * a[j];
				}
			}
		}
		zeta[b] = y / prod;
	}

	// zeta[01] <- mobius(zeta[*1])
	for (ll bi = 0; bi < n; bi++) {
		for (ll b = 0; !(b >> n); b++) {
			// [0] <- [*] - [1]
			if (!(b & (1LL << bi))) {
				zeta[b] = zeta[b] - zeta[b | (1LL << bi)];
			}
		}
	}

	ll ans[n+1];
	for (ll i = 0; i <= n; i++) {
		ans[i] = 0;
	}
	for (ll b = 0; !(b >> n); b++) {
		ll cnt = 0;
		for (ll j = 0; j < n; j++) {
			if (b & (1LL << j)) {
				cnt++;
			}
		}
		ans[cnt] += zeta[b];
	}

	cout << ans[m] << "\n";

	
	return 0;
}