[HAOI2018]奇怪的背包

题目

暴力(dp)好有道理啊

于是我们来个反演吧

考虑一个体积序列({v_1,v_2,...v_n})能凑成(w)的条件

显然是

[v_1x_1+v_2x_2+...+v_nx_nequiv w(mod P) ]

根据贝祖定理,我们知道上面的同余方程有解的条件是

[gcd(v_1,v_2...v_n,P)|w ]

现在题目转化成了求有多少个子集满足(gcd(v_1,v_2..v_n,P)|w)

显然一个暴力(dp)记录一下当前(gcd)转移就好了,由于(gcd)必然是(P)的约数,(P)的约数是(sqrt{P})级别,我们我们只存这些有用的状态,之后每次暴力转移就好了,复杂度(O(nsqrt{P}))就没了

但是我们来反演一波吧

(F(n))表示有多少个子集的(gcd)(n)的倍数,(f(n))表示有多少个子集的(gcd)(n)

非常显然存在

[F(n)=sum_{n|d}f(d) ]

反演可得

[f(n)=sum_{n|d}F(d)mu(frac{d}{n}) ]

(F(n))我们随便搞一搞就好了,若有(s)个数都有(n)这个约数,那么(F(n)=2^s-1)

之后反演得到(f),这里直接枚举(P)的约数,复杂度是(O(sigma^2(P))),但是我们还需要求一个(mu),配合线筛预处理还是挺快的

最后我们的答案是

[Ans(n)=sum_{d|n}f(d) ]

我们累加一遍就好了

复杂度很玄学,毕竟有一个(mu)需要求,可能是(O(sigma^2(P)sqrt{P}))

代码

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#define re register
#define LL long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
const int maxn=1e6+5;
const int mod=1e9+7;
const int M=2e5+5;
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;
}
int d[M],F[M],ans[M],f[M],g[M];
int pw[maxn],a[maxn];
int n,m,P,Sqr,tot;
int is[M],p[M>>1],mu[M],id2[M],id1[M];
inline int find(int x) {return (x<=Sqr)?id1[x]:id2[P/x];}
inline int gcd(int a,int b) {return !b?a:gcd(b,a%b);}
inline void Init() {
    is[1]=1;mu[1]=1;
    for(re int i=2;i<=Sqr;i++) {
        if(!is[i]) p[++p[0]]=i,mu[i]=-1;
        for(re int j=1;j<=p[0]&&p[j]*i<=Sqr;j++) {
            is[p[j]*i]=1;
            if(i%p[j]==0) break;
            mu[p[j]*i]=-1*mu[i];
        }
    }
    for(re int i=1;i<=tot;i++) 
    if(d[i]<=Sqr) id1[d[i]]=i;
        else id2[P/d[i]]=i;
}
inline int getMu(int x) {
    if(x<=Sqr) return mu[x];
    int now=1,t=x;
    for(re int i=1;p[i]*p[i]<=t&&i<=p[0];i++)
    if(x%p[i]==0) {
        x/=p[i];
        if(x%p[i]==0) return 0;
        now=-1*now;
        if(x==1) break;
    }
    if(x!=1) now=-1*now;
    return now;
}
int main() {
    n=read(),m=read(),P=read();
    for(re int i=1;i*i<=P;i++)
    if(P%i==0) {
        d[++tot]=i;
        if(i!=P/i) d[++tot]=P/i;
    }
    pw[0]=1;
    for(re int i=1;i<=n;i++) pw[i]=(pw[i-1]+pw[i-1])%mod;
    for(re int i=1;i<=n;i++) a[i]=read();
    std::sort(d+1,d+tot+1);
    Sqr=std::sqrt(P);Init();
    for(re int i=1;i<=n;i++) {
        int k=gcd(a[i],P);
        g[find(k)]++;
    }
    for(re int i=1;i<=tot;i++)
        for(re int j=i;j<=tot;j++) {
            if(d[j]%d[i]) continue;
            F[i]+=g[j];
        }
    for(re int i=1;i<=tot;i++) F[i]=pw[F[i]]-1;
    for(re int i=1;i<=tot;i++)
        for(re int j=i;j<=tot;j++) {
            if(d[j]%d[i]) continue;
            f[i]+=F[j]*getMu(d[j]/d[i]);
            f[i]%=mod;
            f[i]=(f[i]+mod)%mod;
        }
    for(re int i=1;i<=tot;i++)
        for(re int j=i;j<=tot;j++) {
            if(d[j]%d[i]) continue;
            ans[j]=(ans[j]+f[i])%mod;
        }
    for(re int i=1;i<=m;i++) {
        int x=read();
        int k=gcd(x,P);
        printf("%d
",ans[find(k)]);
    }
    return 0;
}
原文地址:https://www.cnblogs.com/asuldb/p/10592662.html