洛谷P5158 【模板】多项式快速插值
https://www.luogu.org/problemnew/show/P5158
题解：https://www.cnblogs.com/zzqsblog/p/7923192.html
版本1：基于版本1
1 #prag\ 2 ma GCC optimize(2) 3 #include<cstdio> 4 #include<algorithm> 5 #include<cstring> 6 #include<vector> 7 #include<cmath> 8 using namespace std; 9 #define fi first 10 #define se second 11 #define mp make_pair 12 #define pb push_back 13 typedef long long ll; 14 typedef unsigned long long ull; 15 const int md=998244353; 16 const int N=262144; 17 #define delto(a,b) ((a)-=(b),((a)<0)&&((a)+=md)) 18 inline int del(int a,int b) 19 { 20 a-=b; 21 return a<0?a+md:a; 22 } 23 int rev[N]; 24 void init(int len) 25 { 26 int bit=0,i; 27 while((1<<(bit+1))<=len) ++bit; 28 for(i=1;i<len;++i) 29 rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); 30 } 31 ull poww(ull a,ull b) 32 { 33 ull ans=1; 34 for(;b;b>>=1,a=a*a%md) 35 if(b&1) 36 ans=ans*a%md; 37 return ans; 38 } 39 int inv[300011],pw_1[300011],pw_2[300011]; 40 void dft(int *a,int len,int idx)//要求len为2的幂 41 { 42 int i,j,k,t1,t2;ull wn,wnk; 43 for(i=0;i<len;++i) 44 if(i<rev[i]) 45 swap(a[i],a[rev[i]]); 46 for(i=1;i<len;i<<=1) 47 { 48 wn=idx==1?pw_1[i]:pw_2[i]; 49 //wn=poww(idx==1?3:332748118,(md-1)/(i<<1)); 50 for(j=0;j<len;j+=(i<<1)) 51 { 52 wnk=1; 53 for(k=j;k<j+i;++k,wnk=wnk*wn%md) 54 { 55 t1=a[k];t2=a[k+i]*wnk%md; 56 a[k]+=t2; 57 (a[k]>=md)&&(a[k]-=md); 58 a[k+i]=t1-t2; 59 (a[k+i]<0)&&(a[k+i]+=md); 60 } 61 } 62 } 63 if(idx==-1) 64 { 65 ull ilen=inv[len]; 66 for(i=0;i<len;++i) 67 a[i]=a[i]*ilen%md; 68 } 69 } 70 void p_inv(int *f,int *g,int len)//g=f^(-1);f,g数组的长度不小于2len(需要足够长用于临时存放元素);要求len是2的幂 71 { 72 static int t1[N],t2[N]; 73 g[0]=poww(f[0],md-2); 74 for(int i=2,j;i<=len;i<<=1) 75 { 76 memcpy(t1,f,sizeof(int)*i); 77 memcpy(t2,g,sizeof(int)*(i>>1)); 78 memset(t2+(i>>1),0,sizeof(int)*(i>>1)); 79 init(i); 80 dft(t1,i,1);dft(t2,i,1); 81 for(j=0;j<i;++j) 82 t1[j]=ull(t1[j])*t2[j]%md; 83 dft(t1,i,-1); 84 for(j=0;j<(i>>1);++j) 85 t1[j]=t1[j+(i>>1)]; 86 memset(t1+(i>>1),0,sizeof(int)*(i>>1)); 87 dft(t1,i,1); 88 for(j=0;j<i;++j) 89 t1[j]=ull(t1[j])*t2[j]%md; 90 dft(t1,i,-1); 91 for(j=i>>1;j<i;++j) 92 g[j]=md-t1[j-(i>>1)]; 93 } 94 } 95 inline void p_de(int *f,int len)//derivative求导;f=f' 96 { 97 for(int i=0;i<len-1;++i) 98 f[i]=ull(i+1)*f[i+1]%md; 99 f[len-1]=0; 100 } 101 inline void p_in(int *f,int len)//integral积分;f=?f 102 { 103 for(int i=len-1;i>=1;--i) 104 f[i]=ull(f[i-1])*inv[i]%md; 105 f[0]=0; 106 } 107 void p_ln(int *f,int len)//要求len为2的幂，f[0]=1 108 { 109 static int t3[N]; 110 p_inv(f,t3,len);p_de(f,len); 111 init(len<<1); 112 dft(f,len<<1,1);dft(t3,len<<1,1); 113 for(int i=0;i<(len<<1);++i) 114 f[i]=ull(f[i])*t3[i]%md; 115 dft(f,len<<1,-1);p_in(f,len); 116 } 117 void p_exp(int *f,int *g,int len)//要求len为2的幂，f[0]=0 118 { 119 static int t1[N],t2[N]; 120 g[0]=1; 121 for(int i=2,j;i<=len;i<<=1) 122 { 123 memcpy(t1,g,sizeof(int)*(i>>1)); 124 memset(t1+(i>>1),0,sizeof(int)*(i>>1)); 125 p_ln(t1,i); 126 for(j=0;j<(i>>1);++j) 127 t1[j]=del(f[j+(i>>1)],t1[j+(i>>1)]); 128 memset(t1+(i>>1),0,sizeof(int)*(i>>1)); 129 init(i); 130 dft(t1,i,1); 131 memcpy(t2,g,sizeof(int)*(i>>1)); 132 memset(t2+(i>>1),0,sizeof(int)*(i>>1)); 133 dft(t2,i,1); 134 for(j=0;j<i;++j) 135 t1[j]=ull(t1[j])*t2[j]%md; 136 dft(t1,i,-1); 137 for(j=i>>1;j<i;++j) 138 g[j]=t1[j-(i>>1)]; 139 } 140 } 141 void p_div(int *a,int *b,int *c,int n,int m)//c=a/b;deg(a)=n,deg(b)=m,deg(c)=n-m;a,b无前导0;n>=m 142 { 143 reverse(a,a+n+1);reverse(b,b+m+1); 144 int x=n-m+1,t=1; 145 for(;t<x;t<<=1); 146 memset(b+m+1,0,sizeof(int)*max(t-m-1,0)); 147 p_inv(b,c,t); 148 memset(c+x,0,sizeof(int)*((t<<1)-x)); 149 memset(a+x,0,sizeof(int)*((t<<1)-x)); 150 init(t<<1); 151 dft(a,t<<1,1);dft(c,t<<1,1); 152 for(int i=0;i<(t<<1);++i) 153 c[i]=ull(c[i])*a[i]%md; 154 dft(c,t<<1,-1); 155 memset(c+(n-m+1),0,sizeof(int)*((t<<1)-n+m-1)); 156 reverse(c,c+x); 157 } 158 void p_divmod(int *a,int *b,int *c,int *d,int n,int m)//c=a/b,d=a%b,deg(d)=(<=)m-1;其余同上 159 { 160 static int t1[N]; 161 memcpy(d,a,sizeof(int)*(m+1)); 162 int x=n+1,t=1; 163 for(;t<x;t<<=1); 164 memcpy(t1,b,sizeof(int)*(m+1)); 165 memset(t1+m+1,0,sizeof(int)*max(t-m-1,0)); 166 p_div(a,b,c,n,m); 167 memcpy(a,c,sizeof(int)*(n-m+1)); 168 memset(a+n-m+1,0,sizeof(int)*(t-n+m-1)); 169 init(t); 170 dft(a,t,1);dft(t1,t,1); 171 for(int i=0;i<t;++i) 172 t1[i]=ull(t1[i])*a[i]%md; 173 dft(t1,t,-1); 174 for(int i=0;i<=m;++i) 175 delto(d[i],t1[i]); 176 } 177 namespace P_me 178 { 179 int *ta[N];//用线段树的方法给递归的每一层一个编号，ta[i]表示编号为i的层的P函数的各项系数 180 int data[N*40],*tp;//内存池 181 int *a,*x,*y; 182 #define LC (u<<1) 183 #define RC (u<<1|1) 184 int mt1[N]; 185 const int T=200;//小范围暴力阀值 186 void _p_me1(int l,int r,int u)//计算(x-x_l)(x-x_{l+1})..(x-x_r)并存下来 187 { 188 if(r-l<=T) 189 { 190 int i,j; 191 tp[0]=1; 192 for(i=l;i<=r;++i) 193 { 194 tp[i-l+1]=tp[i-l]; 195 for(j=i-l;j>=1;--j) 196 { 197 tp[j]=(ull(tp[j])*(md-x[i])+tp[j-1])%md; 198 } 199 tp[0]=ull(tp[0])*(md-x[i])%md; 200 } 201 ta[u]=tp;tp+=r-l+2; 202 return; 203 } 204 int mid=(l+r)>>1; 205 _p_me1(l,mid,LC);_p_me1(mid+1,r,RC); 206 int x=r-l+2,t=1;//x=(mid-l+1)+(r-mid)+1 207 for(;t<x;t<<=1); 208 memcpy(mt1,ta[LC],sizeof(int)*(mid-l+2)); 209 memset(mt1+mid-l+2,0,sizeof(int)*(t-mid+l-2)); 210 memcpy(tp,ta[RC],sizeof(int)*(r-mid+1)); 211 memset(tp+r-mid+1,0,sizeof(int)*(t-r+mid-1)); 212 init(t); 213 dft(mt1,t,1);dft(tp,t,1); 214 for(int i=0;i<t;++i) 215 tp[i]=ull(tp[i])*mt1[i]%md; 216 dft(tp,t,-1); 217 ta[u]=tp;tp+=r-l+2; 218 } 219 int mt2[N],mt3[N]; 220 void _p_me2(int *a,int n,int l,int r,int u)//a是A的系数,deg(A)<=n；求A(x_l)到A(x_r),放入y_l到y_r 221 { 222 if(r-l<=T) 223 { 224 int t,i,j; 225 for(i=l;i<=r;++i) 226 { 227 t=a[n]; 228 for(j=n-1;j>=0;--j) 229 t=(ull(t)*x[i]+a[j])%md; 230 y[i]=t; 231 } 232 return; 233 } 234 int x=(n+1)<<1,t=1; 235 for(;t<x;t<<=1); 236 int mt4[t];//根据需要改成new? 237 int mid=(l+r)>>1,n1; 238 memcpy(mt1,a,sizeof(int)*(n+1)); 239 for(n1=n;n1>=0 && mt1[n1]==0;) --n1; 240 if(n1<0) 241 { 242 memset(y+l,0,sizeof(int)*(r-l+1)); 243 return; 244 } 245 memcpy(mt2,ta[LC],sizeof(int)*(mid-l+2)); 246 if(n1<mid-l+1) 247 { 248 memcpy(mt4,mt1,sizeof(int)*(n1+1)); 249 _p_me2(mt4,n1,l,mid,LC); 250 } 251 else 252 { 253 p_divmod(mt1,mt2,mt3,mt4,n1,mid-l+1); 254 _p_me2(mt4,mid-l,l,mid,LC); 255 } 256 memcpy(mt1,a,sizeof(int)*(n+1)); 257 for(n1=n;n1>=0 && mt1[n1]==0;) --n1; 258 memcpy(mt2,ta[RC],sizeof(int)*(r-mid+1)); 259 if(n1<r-mid) 260 { 261 memcpy(mt4,mt1,sizeof(int)*(n1+1)); 262 _p_me2(mt4,n1,mid+1,r,RC); 263 } 264 else 265 { 266 p_divmod(mt1,mt2,mt3,mt4,n1,r-mid); 267 _p_me2(mt4,r-mid-1,mid+1,r,RC); 268 } 269 } 270 void p_multieval(int *a0,int *x0,int *y0,int n,int m)//deg(a)=n,x有m个数 271 { 272 tp=data; 273 a=a0;x=x0;y=y0; 274 _p_me1(0,m-1,1); 275 _p_me2(a,n,0,m-1,1); 276 } 277 int mi1[N],mi2[N],mi3[N]; 278 void _p_in1(int *ans,int l,int r,int u) 279 { 280 if(r-l<=T) 281 { 282 int i,j;int *tau=ta[u]; 283 memset(ans,0,sizeof(int)*(r-l+1)); 284 for(i=l;i<=r;++i) 285 { 286 mi3[r-l]=tau[r-l+1]; 287 for(j=r-l-1;j>=0;--j) 288 mi3[j]=(ull(mi3[j+1])*x[i]+tau[j+1])%md; 289 for(j=0;j<=r-l;++j) 290 ans[j]=(ans[j]+ull(mi3[j])*mi2[i])%md; 291 } 292 return; 293 } 294 int mid=(l+r)>>1; 295 int x=r-l+1,t=1,i; 296 for(;t<x;t<<=1); 297 int mia[t],mib[t];//视情况改为new? 298 _p_in1(mia,l,mid,LC); 299 memset(mia+mid-l+1,0,sizeof(int)*(t-mid+l-1)); 300 _p_in1(mib,mid+1,r,RC); 301 memset(mib+r-mid,0,sizeof(int)*(t-r+mid)); 302 init(t); 303 dft(mia,t,1);dft(mib,t,1); 304 memcpy(mi3,ta[LC],sizeof(int)*(mid-l+2)); 305 memset(mi3+mid-l+2,0,sizeof(int)*(t-mid+l-2)); 306 dft(mi3,t,1); 307 for(i=0;i<t;++i) 308 ans[i]=ull(mib[i])*mi3[i]%md; 309 memcpy(mi3,ta[RC],sizeof(int)*(r-mid+1)); 310 memset(mi3+r-mid+1,0,sizeof(int)*(t-r+mid-1)); 311 dft(mi3,t,1); 312 for(i=0;i<t;++i) 313 ans[i]=(ans[i]+ull(mia[i])*mi3[i])%md; 314 dft(ans,t,-1); 315 } 316 void p_inter(int *x0,int *y0,int *a0,int n) 317 { 318 x=x0;tp=data; 319 _p_me1(0,n-1,1); 320 memcpy(mi1,ta[1],sizeof(int)*(n+1)); 321 p_de(mi1,n+1); 322 y=mi2; 323 _p_me2(mi1,n-1,0,n-1,1); 324 for(int i=0;i<n;++i) 325 mi2[i]=ull(poww(mi2[i],md-2))*y0[i]%md; 326 _p_in1(a0,0,n-1,1); 327 } 328 } 329 using P_me::p_multieval; 330 using P_me::p_inter; 331 int a[N],x[N],y[N]; 332 int n; 333 int main() 334 { 335 int i; 336 inv[1]=1; 337 for(i=2;i<=300000;++i) 338 inv[i]=ull(md-md/i)*inv[md%i]%md; 339 for(i=1;i<300000;i<<=1) 340 { 341 pw_1[i]=poww(3,(md-1)/(i<<1)); 342 pw_2[i]=poww(332748118,(md-1)/(i<<1)); 343 } 344 scanf("%d",&n); 345 for(i=0;i<n;++i) 346 scanf("%d%d",x+i,y+i); 347 p_inter(x,y,a,n); 348 for(i=0;i<n;++i) 349 printf("%d ",a[i]); 350 return 0; 351 }
