状压 FWT P5366题解

· · 题解

第一次看见这题的时候觉得挺牛逼的，然后想了半个小时发现其实还是挺简单（

首先：A=\max_{n=1}^{10^8}\sigma(n)=768,B=\max_{n=1}^{10^8}\omega(n)=8。

首先 G,L 都没啥必要，因为可以变成 1,n。

那么 X/G 只能是 n 的因数。问题变成了在 n 的因数中钦定某个数必选，然后要使 \gcd=1,\operatorname{lcm}=n。

观察一下可以知道 \gcd,\operatorname{lcm} 都是套在表面上的。设 n=\prod p_i^{k_i}，把每个数状压一下，设这个数为 \prod P_i^{K_i}，第 i 位为 0 表示 0<K_i<k_i，2 表示 K_i=k_i，1 表示 K_i=0。然后设 f_{S1,S2} 表示选择一部分数，其 \gcd 的状态为 S1，\operatorname{lcm} 的状态为 S2 的方案数。

然后这玩意儿看上去像是在对“是否拥有一个完整的 p_i^{k_i}”做 Or 背包的同时对“是否完全不包含 p_i”做一个 Or 背包。如果不钦定是可以做到 A4^B 的。

注意到这个 A4^B 的代价我们是可以接受的。所以我们考虑做一个“逆元”的操作，把钦定的数字从其中删去，这个时候得到的相当于是钦定这个数不选的结果。然后用之前的减去这个结果就行了。（就是用所有方案数减去强行删掉这个数后的方案数）

这玩意儿好像不太可做。但是好像把 S1,S2 拼接在一起到一起是没问题的。。。

所以重新设状态，f_{S1\times 2^k+S2} 表示上面的 f_{S1,S2}。

然后这里就是纯纯 Or 卷积了。并且可以发现逆元是将这个数做 FWT 后所有位置取逆元，能够发现显然是 2^{-1}。

需要注意的是要以点值形式保存，然后做一个子集反演（等价于一次 IFWT）。复杂度是单次 O(4^B),可以接受。

不过嘛，可以考虑计算每一个数对最开始背包的贡献。做一个子集和差不多了。

复杂度 O(A4^B)。

然后我们发现一件事：

统计不包含某个数的方案数时，实际上是将超集全部除以了 2。

然后使用最开始的值减去这一部分，就只剩下了超集的部分除以 2。

由于我们是每次做 IFWT 后取在 0 处的位置，等价于每个数带了一个权值。

于是全部除以 2 后做一个超集求和即可。复杂度是十分优秀的 O(B4^B+A)。

#include<cstdio>
#include<cctype>
typedef unsigned ui;
const ui M=1<<8,mod=1e9+7,inv2=500000004;
ui Q,G,L,X,N,pw2[1024];ui m,p[10],k[10],t[10];double invG,x[10],y[10];ui f[M*M];bool ppc[M*M];
char _INPUT[1<<21|1],_OUTPUT[1<<21|1],*_P1=_INPUT,*_P2=_OUTPUT;
inline void DFS(const ui&id,const ui&S1,const ui&S2,ui x){
    if(id==m)return void(++f[S1<<m|S2]);
    for(ui i=0;x<=N&&i<=k[id];x*=p[id],++i)DFS(id+1,S1|(i==k[id])<<id,S2|(i==0)<<id,x);
}
inline void Divide(ui N){
    for(ui i=2;i*i<=N;++i)if(!(N%i)){
        p[m]=i;while(!(N%i))++k[m],N/=i;++m;
    }
    if(N^1)p[m]=N,k[m]=1,++m;
    for(ui i=m-1;~i;--i){
        x[i]=1./p[i]+1e-15;t[i]=y[i]=1;for(ui j=1;j<=k[i];++j)y[i]*=x[i],t[i]*=p[i];
    }
}
inline ui read(){
    ui n(0);char s;while(!isdigit(s=*_P1++));while(n=n*10+(s&15),isdigit(s=*_P1++));return n;
}
inline void write(ui n){
    static char s[10];ui top(0);while(s[++top]=n%10^48,n/=10);while(*_P2++=s[top],--top);*_P2++=10;
}
signed main(){
    fread(_INPUT,1,1<<21,stdin);
    N=read();G=read();L=read();Q=read();invG=1./G+1e-15;
    if(L%G){
        while(Q--)putchar(48),putchar(10);return 0;
    }
    N*=invG;L*=invG;Divide(L);DFS(0,0,0,1);pw2[0]=1;for(ui i=1;i^1024;++i)pw2[i]=2*pw2[i-1]%mod;
    const ui&M=1<<m*2;
    for(ui len=1;len^M;len<<=1)for(ui k=0;k^M;k+=len<<1)for(ui i=0;i^len;++i)f[i|k|len]+=f[i|k];
    for(ui S=0;S^M;++S)f[S]&&--f[S],f[S]=S&&(ppc[S]=ppc[S>>1]^(S&1))?mod-pw2[f[S]]:pw2[f[S]];
    for(ui len=1;len^M;len<<=1)for(ui k=0;k^M;k+=len<<1)for(ui i=0;i^len;++i)f[i|k]=(f[i|k]+f[i|k|len])%mod;
    while(Q--){
        X=read();
        if(ui(X*invG)*G^X){
            *_P2++=48;*_P2++=10;continue;
        }
        if(X*=invG,L%X||X>N){
            *_P2++=48;*_P2++=10;continue;
        }
        ui S1(0),S2(0);
        for(ui i=0;i^m;++i)S1|=(ui(X*y[i])*t[i]==X)<<i,S2|=(ui(X*x[i])*p[i]!=X)<<i;
        write(f[S1<<m|S2]);continue;
    }
    fwrite(_OUTPUT,1,_P2-_OUTPUT,stdout);
}