将 \(L\) 唯一分解为 \(p_1^{a_1}p_2^{a_2}...p_k^{a_k}\)

对 G 也分解为 \(p_1^{b_1}p_2^{b_2}...p_k^{b_k}\)。

称 \(a_i , b_i\) 分别为 \(p_i\) 这个质因子幂次的上下界。

显然为了满足 \(gcd\) 为 \(G\) 且 \(lcm\) 为 \(L\),对于每一个 \(p_i\) ,就至少要

有一个数触其上界,有一个数触其下界。

那么就可以拿一个长为 \(2^k\) 的二进制状态表示每个质因子是否已

触其上界,是否已触其下界。

然后前后缀对于 \(L\) 的约数且是 \(G\) 的倍数的数作一遍 \(DP\)

直接 \(FWT\) 合并成为答案

还有更加优秀的容斥做法

首先问题已经转化为给定一些集合，求出或为全集的方案数

设 \(f_s\) 表示或为 \(s\) 的方案数

直接求不方便

设 \(g_s=\sum_{i\subset s}f_i\) 表示或为 \(s\) 的子集的方案数

那么 \(g_s=2^{cnt_s}-1\)(\(cnt_s\) 表示 \(s\) 的子集个数)

那么容斥得到 \(f_s=\sum_{i\subset s}(-1)^{|s|-|i|}g_i\)

对于钦定了一定要选集合 \(x\)，只需要把 \(x\) 的超集 \(-1\) 然后再次做上面的容斥

记得 \(g_s=2^{cnt_s-1}\)

# include 
    
     using namespace std;typedef long long ll; namespace IO {    const int maxn(1 << 21 | 1);     char ibuf[maxn], *iS, *iT, c;    int f;     inline char Getc() {        return iS == iT ? (iT = (iS = ibuf) + fread(ibuf, 1, maxn, stdin), (iS == iT ? EOF : *iS++)) : *iS++;    }     template 
     
       inline void In(Int &x) {        for (f = 1, c = Getc(); c < '0' || c > '9'; c = Getc()) f = c == '-' ? -1 : 1;        for (x = 0; c >= '0' && c <= '9'; c = Getc()) x = (x << 1) + (x << 3) + (c ^ 48);        x *= f;    }} using IO :: In; const int maxn(805);const int maxm(16);const int mod(1e9 + 7); int n, g, l, bl, q, p[maxm], mx[maxm], mn[maxm], tot, st, cnt;int size[1 << maxm], ans[maxn], pw[1 << maxm], bit[1 << maxm]; struct Fac {    int x, y;     inline bool operator <(Fac b) const {        return x < b.x;    }} a[maxn]; inline void Inc(int &x, int y) {    x += y;    if (x >= mod) x -= mod;} inline void Pre_Work(int x) {    if (x > n || x % g) return;    register int i, v;    a[++cnt].x = x;    for (i = 0; i < tot; ++i) {        v = 0;        while (x % p[i] == 0) x /= p[i], ++v;        if (v == mn[i]) a[cnt].y |= 1 << i;        if (v == mx[i]) a[cnt].y |= 1 << (i + tot);    }} int main() {    In(n), In(g), In(l), In(q);    register int qq, i, j, x, y, k, tp;    if (l % g) {        for (qq = 1; qq <= q; ++qq) puts("0");        return 0;    }    x = l;    for (i = 2; i * i <= x; ++i)        if (x % i == 0) {            while (x % i == 0) x /= i;            p[tot++] = i;        }    if (x > 1) p[tot++] = x;    st = 1 << (tot << 1), x = g, y = l;    for (i = 0; i < tot; ++i) {        while (x % p[i] == 0) x /= p[i], ++mn[i];        while (y % p[i] == 0) y /= p[i], ++mx[i];    }    for (i = 1; i * i <= l; ++i)        if (l % i == 0) {            Pre_Work(i);            if (i * i != l) Pre_Work(l / i);        }    sort(a + 1, a + cnt + 1);    for (i = 1; i <= cnt; ++i) ++size[a[i].y];    for (i = 1; i < st; i <<= 1)        for (tp = i << 1, j = 0; j < st; j += tp)            for (k = 0; k < i; ++k) size[j + k + i] += size[j + k];    for (pw[0] = 1, j = size[st - 1], i = 1; i <= j; ++i) pw[i] = (pw[i - 1] + pw[i - 1]) % mod;    for (i = 1; i < st; ++i) bit[i] = bit[i >> 1] + (i & 1);    for (i = 1; i <= cnt; ++i)        for (j = 0; j < st; ++j)            if ((a[i].y & j) == a[i].y)                Inc(ans[i], (bit[(st - 1) ^ j] & 1) ? mod - pw[size[j] - 1] : pw[size[j] - 1]);    for (qq = 1; qq <= q; ++qq) {        In(x);        if (l % x || x % g) puts("0");        else {            y = lower_bound(a + 1, a + cnt + 1, (Fac){x, 0}) - a;            a[y].x != x ? puts("0") : printf("%d\n", ans[y]);        }    }    return 0;}