[LOJ2271] [SDOI2017] 遗忘的集合

时间:2022-11-07 04:29:42

题目链接

LOJ:https://loj.ac/problem/2271

洛谷:https://www.luogu.org/problemnew/show/P3784

BZOJ太伤身体死活卡不过还是算了吧...

Solution

为啥窝洛谷\(rk4\) \(\rm BZOJ\)死活跑不过啊...

技不如人,肝败吓疯...


题目并不难

设\(a_i\)表示\(i\)有没有出现在集合中,这是我们要求的答案。

那么把背包写成生成函数就是:

\[\prod_{i=1}^{n}\left(\sum_{j=0}^{+\infty}x^{ij}\right)^{a_i}=\prod_{i=1}^{n}\left(\frac{1}{1-x^i}\right)^{a_i}=F(x)
\]

两边取对数:

\[\sum_{i=1}^{n}a_i\ln \frac{1}{1-x^i}=\ln F(x)
\]

注意到:

\[\ln G(x)=\int \frac{G'(x)}{G(x)}dx
\]

把上面那个分式化一下就是:

\[\ln \frac{1}{1-x^i}=\int \frac{ix^{i-1}}{1-x^i} dx=i\int \sum_{j=0}^{+\infty}x^{ij+i-1}dx=\sum_{j=1}^{+\infty}\frac{1}{j}x^{ij}
\]

带一下:

\[\sum_{i=1}^{n}a_i\sum_{j=1}^{+\infty}\frac{1}{j}x^{ij} =\ln F(x)
\]

我们枚举\(ij\):

\[\sum_{T=1}^{+\infty}x^T\sum_{i|T}a_i\cdot \frac{i}{T}=\ln F(x)\\
\sum_{T=1}^{+\infty}x^T\sum_{i|T}a_i\cdot i=T\cdot \ln F(x)
\]

注意到左边是一个莫反的形式,我们对\(F(x)\)求\(\ln\),然后把系数搞出来,在套莫反公式就好了。

复杂度\(O(n\log n)\)。

#include<bits/stdc++.h>
using namespace std; void read(int &x) {
x=0;int f=1;char ch=getchar();
for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
} void print(int x) {
if(x<0) putchar('-'),x=-x;
if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('\n');} #define lf double
#define ll long long #define pii pair<int,int >
#define vec vector<int > #define pb push_back
#define mp make_pair
#define fr first
#define sc second const int maxn = 6e5+10;
const int inf = 1e9;
const lf eps = 1e-8;
const lf pi = acos(-1);
const int all = (1<<15)-1; struct cp {
lf r,i;
cp () {}
cp (lf _r,lf _i) {r=_r,i=_i;}
cp operator + (const cp &rhs) const {return cp(r+rhs.r,i+rhs.i);}
cp operator - (const cp &rhs) const {return cp(r-rhs.r,i-rhs.i);}
cp operator * (const cp &rhs) const {return cp(r*rhs.r-i*rhs.i,i*rhs.r+r*rhs.i);}
cp operator / (const int &rhs) const {return cp(r/(lf)rhs,i/(lf)rhs);}
void print() {printf("%lf %lf\n",r,i);}
}w[maxn],g[5][maxn]; cp conj(cp x) {return cp(x.r,-x.i);} int n,m,N,bit,mxn,mod;
int pos[maxn],a[maxn],b[maxn],c[maxn],inv[maxn],f[maxn],h[maxn]; void fft(cp *r) {
for(int i=1;i<N;i++) if(pos[i]>i) swap(r[i],r[pos[i]]);
for(int i=1,d=mxn>>1;i<N;i<<=1,d>>=1)
for(int j=0;j<N;j+=i<<1)
for(int k=0;k<i;k++) {
cp x=r[j+k],y=w[k*d]*r[i+j+k];
r[j+k]=x+y,r[i+j+k]=x-y;
}
} void ntt_get(int len) {
for(N=1,bit=0;N<=len;N<<=1,bit++);
for(int i=1;i<N;i++) pos[i]=pos[i>>1]>>1|((i&1)<<(bit-1));
} void poly_mul(int *r,int *s,int *t,int len) {
ntt_get(len);
for(int i=0;i<len;i++) g[0][i]=cp(r[i]&all,r[i]>>15),g[1][i]=cp(s[i]&all,s[i]>>15);
for(int i=len;i<N;i++) g[0][i]=g[1][i]=cp(0,0);
fft(g[0]),fft(g[1]);
for(int i=0;i<N;i++) {
int j=(N-i)&(N-1);
g[2][j]=(g[0][i]+conj(g[0][j]))*cp(0.5,0)*g[1][i];
g[3][j]=(g[0][i]-conj(g[0][j]))*cp(0,-0.5)*g[1][i];
}fft(g[2]),fft(g[3]);
for(int i=0;i<N;i++) g[2][i]=g[2][i]/N,g[3][i]=g[3][i]/N;
for(int i=0;i<N;i++) {
ll pp=g[2][i].r+0.5,x=g[2][i].i+0.5,y=g[3][i].r+0.5,z=g[3][i].i+0.5;
t[i]=((pp%mod+(((x+y)%mod)<<15)+((z%mod)<<30))%mod+mod)%mod;
}
} int add(int x,int y) {return x+y>mod?x+y-mod:x+y;}
int del(int x,int y) {return x-y<0?x-y+mod:x-y;}
int mul(int x,int y) {return 1ll*x*y-1ll*x*y/mod*mod;} int tmp[10][maxn]; int qpow(int aa,int x) {
int res=1;
for(;x;x>>=1,aa=mul(aa,aa)) if(x&1) res=mul(res,aa);
return res;
} void poly_inv(int *r,int *t,int len) {
if(len==1) return t[0]=qpow(r[0],mod-2),void();
poly_inv(r,t,len>>1);ntt_get(len);
poly_mul(t,r,tmp[2],len);
poly_mul(tmp[2],t,tmp[3],len);
for(int i=0;i<len;i++) t[i]=del(add(t[i],t[i]),tmp[3][i]);
} void poly_der(int *r,int *t,int len) {
ntt_get(len);
for(int i=1;i<len;i++) t[i-1]=mul(i,r[i]);
for(int i=len-1;i<N;i++) t[i]=0;
} void poly_int(int *r,int *t,int len) {
ntt_get(len);
for(int i=0;i<len;i++) t[i+1]=mul(inv[i+1],r[i]);t[0]=0;
for(int i=len+1;i<N;i++) t[i]=0;
} void poly_ln(int *r,int *t,int len) {
poly_der(r,tmp[4],len);
poly_inv(r,tmp[5],len);
poly_mul(tmp[4],tmp[5],tmp[6],len);
poly_int(tmp[6],t,len);
} int pri[maxn],tot,mu[maxn],A[maxn],vis[maxn]; void sieve() {
mu[1]=1;
for(int i=2;i<=mxn;i++) {
if(!vis[i]) mu[i]=-1,pri[++tot]=i;
for(int j=1;j<=tot&&i*pri[j]<=mxn;j++) {
vis[i*pri[j]]=1;
if(i%pri[j]==0) break;
mu[i*pri[j]]=-mu[i];
}
}
} int main() {
read(n),read(mod);ntt_get(n<<1);inv[0]=inv[1]=1;
for(mxn=1;mxn<N;mxn<<=1);
for(int i=0;i<=mxn;i++) w[i]=cp(cos(2.0*pi*i/mxn),sin(2.0*pi*i/mxn));
for(int i=2;i<N;i++) inv[i]=mul(mod-mod/i,inv[mod%i]);
for(int i=1;i<=n;i++) read(f[i]);f[0]=1;
poly_ln(f,h,N>>1);sieve();
for(int i=1;i<=n;i++) h[i]=mul(h[i],i);
for(int i=1;i<=n;i++)
for(int j=i;j<=n;j+=i)
A[j]+=h[j/i]*mu[i];
int ans=0;
for(int i=1;i<=n;i++) if(A[i]) ans++;
write(ans);
for(int i=1;i<=n;i++) if(A[i]) printf("%d ",i);puts("");
return 0;
}