萌新求卡常

P3784 [SDOI2017] 遗忘的集合

```cpp #include<cstdio> #include<iostream> #include<cstring> #include<cmath> #include<algorithm> #include<bitset> #include<vector> #define N 524293 #define ll long long #define reg register using namespace std; const double pi = acos(-1); const int LIM = 1<<19; struct complex{ double x,y; complex(double x=0,double y=0):x(x),y(y){} inline complex operator + (const complex& b) const{ return complex(x+b.x,y+b.y); } inline complex operator - (const complex& b) const{ return complex(x-b.x,y-b.y); } inline complex operator * (const complex& b) const{ return complex(x*b.x-y*b.y,x*b.y+y*b.x); } inline complex operator / (const int& b) const{ return complex(x/b,y/b); } inline complex operator ~ () const{ return complex(x,-y); } }; struct poly{ int a[N]; int t; }; complex rt[N]; int lg2[N],rev[N],mu[N],pr[N],inv[N]; bitset<N> vis; inline void init(); inline void read(int &x); void print(int x); inline void print_poly(poly f); inline int add(int a,int b); inline int dec(int a,int b); inline int power(int a,int t); inline void FFT(complex *a,int type,int lim); poly multiply(poly A,poly B,int len); poly square(poly A,int len); poly inverse(poly f); inline poly derivative(poly f); inline poly integral(poly f); poly log(poly f); poly F; int a[N]; int n,p,ans; int main(){ int cnt = 0; read(n),read(p); init(); F.t = n,F.a[0] = 1; for(reg int i=1;i<=n;++i) read(F.a[i]); F = log(F); for(reg int i=1;i<=n;++i) F.a[i] = (ll)F.a[i]*i%p; for(reg int i=1;i<=n;++i){ cnt = 0; for(reg int j=i;j<=n;j+=i){ ++cnt; if(!mu[i]) continue; if(mu[i]>0) a[j] += F.a[cnt]; else a[j] -= F.a[cnt]; } } for(reg int i=1;i<=n;++i) if(a[i]) ++ans; print(ans),putchar('\n'); for(reg int i=1;i<=n;++i) if(a[i]) print(a[i]),putchar(' '); return 0; } poly log(poly f){ return integral(multiply(inverse(f),derivative(f),f.t)); } inline poly derivative(poly f){ for(reg int i=0;i<f.t;++i) f.a[i] = (ll)f.a[i+1]*(i+1)%p; f.a[f.t] = 0; return f; } inline poly integral(poly f){ for(reg int i=f.t;i;--i) f.a[i] = (ll)f.a[i-1]*inv[i]%p; f.a[0] = 0; return f; } poly inverse(poly f){ poly g,h,q; memset(g.a,0,sizeof(g.a)); int top = 0,n = f.t; int s[30]; while(n){ s[++top] = n; n >>= 1; } g.a[0] = power(f.a[0],p-2); while(top--){ g.t = n = s[top+1]; h = f,q = g; for(reg int i=n+1;i<=f.t;++i) h.a[i] = 0; g = multiply(square(g,n),h,n); for(reg int i=0;i<=n;++i) g.a[i] = dec(add(q.a[i],q.a[i]),g.a[i]); } return g; } poly multiply(poly A,poly B,int len){ complex f[N],g[N],h[N],q[N]; reg complex t,f0,f1,g0,g1; reg ll x,y,z; reg int lim = 1; while(lim<=A.t+B.t) lim <<= 1; int l = lg2[lim]-1; for(reg int i=1;i<=lim;++i) rev[i] = (rev[i>>1]>>1)|((i&1)<<l); for(reg int i=0;i<=lim;++i){ f[i] = complex(A.a[i]>>15,A.a[i]&32767); g[i] = complex(B.a[i]>>15,B.a[i]&32767); } FFT(f,1,lim),FFT(g,1,lim); for(reg int i=0;i<=lim;++i){ t = ~f[i?lim-i:0]; f0 = (f[i]-t)*complex(0,-0.5),f1 = (f[i]+t)*0.5; t = ~g[i?lim-i:0]; g0 = (g[i]-t)*complex(0,-0.5),g1 = (g[i]+t)*0.5; h[i] = f1*g1; q[i] = f1*g0 + f0*g1 + f0*g0*complex(0,1); } FFT(h,-1,lim),FFT(q,-1,lim); for(reg int i=0;i<=len;++i){ x = (ll)(h[i].x+0.5)%p<<30; y = (ll)(q[i].x+0.5)<<15; z = q[i].y+0.5; A.a[i] = (x+y+z)%p; } A.t = len; for(reg int i=len+1;i<=lim;++i) A.a[i] = 0; return A; } poly square(poly A,int len){ complex f[N],h[N],q[N]; reg complex t,f0,f1; reg ll x,y,z; int lim = 1; while(lim<=(A.t<<1)) lim <<= 1; int l = lg2[lim]-1; for(reg int i=1;i<=lim;++i) rev[i] = (rev[i>>1]>>1)|((i&1)<<l); for(reg int i=0;i<=lim;++i) f[i] = complex(A.a[i]>>15,A.a[i]&32767); FFT(f,1,lim); for(reg int i=0;i<=lim;++i){ t = ~f[i?lim-i:0]; f0 = (f[i]-t)*complex(0,-0.5),f1 = (f[i]+t)*0.5; h[i] = f1*f1; q[i] = f1*f0*2 + f0*f0*complex(0,1); } FFT(h,-1,lim),FFT(q,-1,lim); for(reg int i=0;i<=len;++i){ x = (ll)(h[i].x+0.5)%p<<30; y = (ll)(q[i].x+0.5)<<15; z = q[i].y+0.5; A.a[i] = (x+y+z)%p; } A.t = len; for(reg int i=len+1;i<=lim;++i) A.a[i] = 0; return A; } inline void FFT(complex *a,int type,int lim){ reg int r,l; for(reg int i=1;i<=lim;++i){ if(i>=rev[i]) continue; swap(a[i],a[rev[i]]); } reg complex w,x,y; l = LIM; for(reg int mid=1;mid<lim;mid<<=1){ r = mid<<1; for(reg int j=0;j<lim;j+=r){ for(reg int k=0;k<mid;++k){ w = rt[l*k]; if(type==-1) w = ~w; x = a[j|k]; y = w*a[j|k|mid]; a[j|k] = x+y; a[j|k|mid] = x-y; } } l >>= 1; } if(type==1) return; for(reg int i=0;i<=lim;++i) a[i] = a[i]/lim; } inline void init(){ int cnt = 0; inv[1] = 1; for(reg int i=2;i<262147;++i) inv[i] = (ll)(p-p/i)*inv[p%i]%p; for(reg int i=0;i<20;++i) lg2[1<<i] = i; rt[0] = complex(1,0); for(reg int i=1;i<N;++i) rt[i] = complex(cos(pi/LIM*i),sin(pi/LIM*i)); mu[1] = 1; for(reg int i=2;i<N;++i){ if(!vis[i]){ pr[++cnt] = i; mu[i] = -1; } for(reg int j=1;j<=cnt;++j){ if(i*pr[j]>=N) break; vis[i*pr[j]] = 1; if(i%pr[j]==0){ mu[i*pr[j]] = 0; break; } mu[i*pr[j]] = -mu[i]; } } } inline int add(int a,int b){ return a+b>=p?a+b-p:a+b; } inline int dec(int a,int b){ return a<b?a-b+p:a-b; } inline int power(int a,int t){ int res = 1; while(t){ if(t&1) res = (ll)res*a%p; a = (ll)a*a%p; t >>= 1; } return res; } inline void print_poly(poly f){ for(reg int i=0;i<=f.t;++i) print(f.a[i]),putchar(' '); putchar('\n'); } inline void read(int &x){ x = 0; char c = getchar(); while(c<'0'||c>'9') c = getchar(); while(c>='0'&&c<='9'){ x = (x<<3)+(x<<1)+(c^48); c = getchar(); } } void print(int x){ if(x>9) print(x/10); putchar(x%10+'0'); } ```
by NaCly_Fish @ 2019-07-01 22:55:41


# 萌新 FaKe
by αnonymous @ 2019-07-01 22:56:22


本来多项式求逆迭代一次需要8次FFT,我搞成了7次结果还是T飞了。。 而且8次FFT的总时间只要了12000ms
by NaCly_Fish @ 2019-07-01 22:57:17


@[NaCly_Fish](/space/show?uid=115864) 去你大爷的萌新
by xht @ 2019-07-01 22:58:27


~~去你的萌新~~
by lukelin @ 2019-07-01 22:59:42


~~这是常数问题吗~~
by lukelin @ 2019-07-01 23:00:51


您是萌新我连人都不配做了
by 犇犇犇犇 @ 2019-07-01 23:08:32


卡常什么的,还是要自己尝试
by Sai0511 @ 2019-07-01 23:12:18


您是萌新我连人都不配做了
by disangan233 @ 2019-07-01 23:34:12


您是萌新我连蒟蒻都不配做了
by lukelin @ 2019-07-01 23:39:54


| 下一页