首先(G,L)肯定会满足(G|L)，否则直接全部输出(0)
之后我们考虑一下能用到的质因数最多只有(8)个
同时我们能选择的数(x)肯定是(L)的约数，还得是(G)的倍数，还不能超过(N)，感性理解一下这样的(x)显然不多，我们直接(dfs)出来
对于每一个数我们把它压成一个(2 imes 8)的二进制数，第(2i-1,2i)位分别表示第(i)个质因子是否到上界/下界，我们直接一个(dp[i][j])表示前(i)个数选择状态为(j)的方案数，同时我们还能处理出一个后缀的(dp)值来
对于一个必须选择的(x)，我们用前缀和后缀做一下或(FWT)，求所有和这个数的状态或起来为全集的位置的和就好了
复杂度看起来不是很科学，但是常数比较小，就跑过去了
代码
#include <bits/stdc++.h>
#define re register
const int mod = 1e9 + 7;
inline int read() {
char c = getchar();
int x = 0;
while (c < '0' || c > '9') c = getchar();
while (c >= '0' && c <= '9') x = (x << 3) + (x << 1) + c - 48, c = getchar();
return x;
}
std::map<int, int> ma;
int N, G, L, Q, M, n, len;
int f[50005], p[20005];
int cnt[2], b[10], a[10][2], l[10], r[10], beg[10];
int c[850], st[850], ans[850], vis[850];
int dp[850][65537], tp[850][65537], d[65537];
void dfs(int x, int now, int S) {
if (x == cnt[0] + 1) {
c[++n] = now, st[n] = S;
return;
}
int t = beg[x];
for (re int i = l[x]; i <= r[x]; ++i)
if (now <= N / t) {
int s = S;
if (i == l[x])
s |= (1 << (2 * x - 2));
if (i == r[x])
s |= (1 << (2 * x - 1));
dfs(x + 1, now * t, s);
t = t * a[x][0];
} else
break;
}
inline int ksm(int a, int b) {
int S = 1;
for (; b; b >>= 1, a = a * a)
if (b & 1)
S = S * a;
return S;
}
inline int qm(int x) { return x < 0 ? x + mod : x % mod; }
inline void Fwt(int *f, int o) {
for (re int i = 2; i <= len; i <<= 1)
for (re int ln = i >> 1, l = 0; l < len; l += i)
for (re int x = l; x < l + ln; ++x) f[x + ln] = qm(f[x + ln] + o * f[x]);
}
int main() {
N = read(), G = read(), L = read(), Q = read();
M = std::ceil(std::sqrt(L));
if (L % G) {
for (re int i = 1; i <= Q; i++) puts("0");
return 0;
}
for (re int i = 2; i <= M; i++) {
if (!f[i])
p[++p[0]] = i;
for (re int j = 1; j <= p[0] && p[j] * i <= M; j++) {
f[p[j] * i] = 1;
if (i % p[j] == 0)
break;
}
}
std::swap(G, L);
for (re int i = 1; i <= p[0]; i++) {
if (G % p[i])
continue;
a[++cnt[0]][0] = p[i];
while (G % p[i] == 0) G /= p[i], ++r[cnt[0]];
}
if (G > 1)
a[++cnt[0]][0] = G, r[cnt[0]] = 1;
for (re int i = 1; i <= p[0]; i++) {
if (L % p[i])
continue;
a[++cnt[1]][1] = p[i];
while (L % p[i] == 0) L /= p[i], ++b[cnt[1]];
}
if (L > 1)
a[++cnt[1]][1] = L, b[cnt[1]] = 1;
for (re int i = 1; i <= cnt[0]; i++)
for (re int j = 1; j <= cnt[1]; j++)
if (a[j][1] == a[i][0]) {
l[i] = b[j];
break;
}
for (re int i = 1; i <= cnt[0]; ++i) beg[i] = ksm(a[i][0], l[i]);
dfs(1, 1, 0);
for (re int i = 1; i <= n; i++) ma[c[i]] = i;
dp[0][0] = 1;
len = 1 << (2 * cnt[0]);
for (re int i = 1; i <= n; i++)
for (re int j = 0; j < len; j++) {
if (!dp[i - 1][j])
continue;
dp[i][j] = (dp[i][j] + dp[i - 1][j]) % mod;
dp[i][j | st[i]] = (dp[i][j | st[i]] + dp[i - 1][j]) % mod;
}
tp[n + 1][0] = 1;
for (re int i = n; i; --i)
for (re int j = 0; j < len; j++) {
if (!tp[i + 1][j])
continue;
tp[i][j] = (tp[i][j] + tp[i + 1][j]) % mod;
tp[i][j | st[i]] = (tp[i][j | st[i]] + tp[i + 1][j]) % mod;
}
while (Q--) {
int x = read(), pos = ma[x];
if (!pos) {
puts("0");
continue;
}
if (vis[pos]) {
printf("%d
", ans[pos]);
continue;
}
Fwt(dp[pos - 1], 1), Fwt(tp[pos + 1], 1);
for (re int i = 0; i < len; i++) d[i] = 1ll * dp[pos - 1][i] * tp[pos + 1][i] % mod;
Fwt(d, -1);
for (re int i = 0; i < len; i++)
if ((i | st[pos]) + 1 == len)
ans[pos] = (ans[pos] + d[i]) % mod;
vis[pos] = 1;
printf("%d
", ans[pos]);
}
return 0;
}