题解
考虑把 $n/$$G$ 和 $L/$$G$ 得到新的 $n,L$ ，现在问题就是从 $[1,n]$ 中选出若干个数使得 $gcd=1,lcm=L$
考虑到 $le 10^8$ 的数最多分解出 $8$ 个质因子，我们设一个数的状态 $(s1,s2)$ 的第 $i$ 位表示第 $i$ 个质因子是否是最低次幂和是否是最高次幂，于是我们直接爆搜出 $L$ 的约数，然后得到每个状态的个数，发现有用的状态不会超过 $600$ 个，于是我们可以考虑状态之间 $dp$ ，转移显然
考虑一个数必选，那也就是剩下的 $dp$ 转移方式是一样的，如果我们能快速得到剩下状态的 $dp$ 值的话那就把这个数代表的状态 $u$ 转移进去就好了，具体实现就是设 $f,g$ 表示 $1-i$ 和 $i-m$ 的 $dp$ 值，于是我们就可以提取两个数组，他们俩之间的转移就直接用 $fwt_{or}$ 转移即可，然后我们得到最终 $dp$ 数组 $F$ ，那现在如果 $x|u$ 是全集的话就可以贡献答案，那对于 $u$ 来说就是 $x$ 的补集的超集，列出式子发现就是 $fwt_{and}$ 的形式，于是我们对 $F$ 再进行 $fwt_{and}$即可
代码
#include <bits/stdc++.h> using namespace std; const int N=605,M=1<<16,P=1e9+7; int n,G,L,Q,A,a[10],b[10],c,t[M],s[N],S,f[N][M],g[N][M]; int X(int x){return x>=P?x-P:x;} int K(int x,int y){ int z=1; for (;y;y>>=1,x=1ll*x*x%P) if (y&1) z=1ll*z*x%P; return z; } void dfs(int x,int v,int s1,int s2){ if (x>c){t[s1|(s2<<c)]++;return;} for (int u=v,i=0;i<=b[x];i++){ dfs(x+1,u,s1|((i==0)<<(x-1)),s2|((i==b[x])<<(x-1))); if (1ll*u*a[x]>n) return;u*=a[x]; } } void init(){ int v=L; for (int i=2;i*i<=v;i++) if (v%i==0){ a[++c]=i; while(v%i==0) v/=i,b[c]++; } if (v>1) a[++c]=v,b[c]=1; dfs(1,1,0,0); } void Fwt(int *a,int o){ for (int i=1;i<A;i<<=1) for (int j=0;j<A;j+=(i<<1)) for (int k=0;k<i;k++) if (o) a[i+j+k]=X(a[i+j+k]+a[j+k]); else a[i+j+k]=X(a[i+j+k]-a[j+k]+P); } void fwt(int *a){ for (int i=1;i<A;i<<=1) for (int j=0;j<A;j+=(i<<1)) for (int k=0;k<i;k++) a[j+k]=X(a[i+j+k]+a[j+k]); } int get(int x){ int u=0; for (int v,i=1;i<=c;i++){ v=0; while(x%a[i]==0) v++,x/=a[i]; if (v==0) u|=(1<<(i-1)); if (v==b[i]) u|=(1<<(i-1+c)); } return u; } int main(){ cin>>n>>G>>L>>Q; if (L%G){ for (;Q--;) puts("0"); return 0; } n/=G;L/=G;init();A=(1<<(c+c)); for (int i=0;i<A;i++) if (t[i]) s[++S]=i,t[i]=K(2,t[i])-1; f[0][0]=g[S+1][0]=1; for (int i=1;i<=S;i++){ for (int j=0;j<A;j++) f[i][j]=f[i-1][j]; for (int j=0;j<A;j++) f[i][j|s[i]]=X(f[i][j|s[i]]+1ll*f[i-1][j]*t[s[i]]%P); } for (int i=S;i;i--){ for (int j=0;j<A;j++) g[i][j]=g[i+1][j]; for (int j=0;j<A;j++) g[i][j|s[i]]=X(g[i][j|s[i]]+1ll*g[i+1][j]*t[s[i]]%P); } for (int i=0;i<=S;i++) Fwt(f[i],1); for (int i=1;i<=S+1;i++) Fwt(g[i],1); for (int i=0;i<S;i++) for (int j=0;j<A;j++) f[i][j]=1ll*f[i][j]*g[i+2][j]%P; for (int i=0;i<S;i++) Fwt(f[i],0),fwt(f[i]); for (int x,u,y;Q--;){ scanf("%d",&x); if (x%G){puts("0");continue;}x/=G; if (L%x || x>n){puts("0");continue;} u=get(x);y=lower_bound(s+1,s+S+1,u)-s; printf("%lld ",1ll*f[y-1][(A-1)^u]*(t[u]+1)%P*((P+1)>>1)%P); } return 0; }