FFT
板子
点此查看代码
#include <bits/stdc++.h>
using namespace std;
const int maxn=2097160;
const int G=3,md=998244353;
int n,m,i,j,k,u,v,p,t,le;
int f[maxn],g[maxn],r[maxn];
inline int Pow(int d,int z){
int r=1;
do{
if(z&1) r=1LL*r*d%md;
d=1LL*d*d%md;
}while(z>>=1); return r;
}
inline void Add(int &x,int y){x-=((x+=y)>=md)*md;}
void dft(int *f){
for(i=1;i<t;++i) if(i<r[i]) swap(f[i],f[r[i]]);
for(i=1,le=2;le<=t;i=le,le<<=1){
for(j=0,v=Pow(G,(md-1)/le);j<t;j+=le){
for(k=0,u=1;k<i;++k,u=1LL*u*v%md){
p=1LL*f[i+j+k]*u%md;
Add(f[i+j+k]=f[j+k],md-p);
Add(f[j+k],p);
}
}
}
}
int main(){
scanf("%d%d",&n,&m);
for(i=0;i<=n;++i) scanf("%d",f+i);
for(i=0;i<=m;++i) scanf("%d",g+i);
for(n+=m,t=1;t<=n;t<<=1);
for(i=1;i<t;++i) r[i]=(r[i>>1]>>1)|((i&1)*(t>>1));
dft(f); dft(g); for(i=0;i<t;++i) f[i]=1LL*f[i]*g[i]%md;
dft(f); reverse(f+1,f+t); p=Pow(t,md-2);
for(i=0;i<=n;++i) printf("%d ",1LL*f[i]*p%md);
return 0;
}
属于需要背的知识。原根表。常见的原根形如 \(998244353\) 原根为 \(3\),\(924844033\) 原根为 \(5\)。
分治 FFT
例题:P4721 【模板】分治 FFT
令 \(m=\left\lfloor\frac{l+r}2\right\rfloor\)。考虑在知道 \(f_l\sim f_m\) 时,其对 \(f_{m+1}\sim f_r\) 会造成怎样的贡献。此时 \(\forall x\in[m+1,r];f_x\gets\sum_{i=l}^mf_ig_{x-i}\),此时将 \(f_l\sim f_m\) 和 \(g_1\sim g_{r-l}\) 卷积可得对应的贡献,分治下去计算其他部分的贡献即可。
使用多项式求逆:考虑 \(f,g\) 的生成函数 \(F(x)=\sum_{i\ge 0}f_ix^i,G(x)=\sum_{i\ge 0}g_ix^i\),则由 \(\forall i>0,f_i=\sum_{j=0}^{i-1}f_jg_{i-j}\) 有
\]
由于 \(g_0\) 没有对 \(f\) 产生贡献,所以可以令 \(g_0=0\),则 \(F(x)=1+F(x)G(x)=\frac 1{1-G(x)}\)。可以多项式求逆解决。
代码
点此查看代码
#include <bits/stdc++.h>
using namespace std;
const int maxn=100010;
const int maxb=262150;
const int md=998244353,_g=3;
int n,i,j,k,u,v,t,p,le,f[maxn],g[maxn];
int r[maxb],F[maxb],G[maxb];
inline int Pow(int d,int z){
int r=1;
do{
if(z&1) r=1LL*r*d%md;
d=1LL*d*d%md;
}while(z>>=1); return r;
}
inline void Add(int &x,int y){x-=((x+=y)>=md)*md;}
void dft(int *f){
for(i=1;i<t;++i) if(i<r[i]) swap(f[i],f[r[i]]);
for(i=1,le=2;le<=t;i=le,le<<=1){
v=Pow(_g,(md-1)/le);
for(j=0;j<t;j+=le){
for(k=0,u=1;k<i;++k,u=1LL*u*v%md){
p=1LL*f[i+j+k]*u%md;
Add(f[i+j+k]=f[j+k],md-p);
Add(f[j+k],p);
}
}
}
}
void Solve(int l,int r){
if(l==r) return; int m=(l+r)>>1; Solve(l,m);
for(t=1,k=r-l+m-l-1;t<=k;t<<=1);
memset(F,0,t<<2); memset(G,0,t<<2);
memcpy(F,f+l,(m-l+1)<<2); memcpy(G,g+1,(r-l)<<2);
for(i=1;i<t;++i) ::r[i]=(::r[i>>1]>>1)|((i&1)*(t>>1));
dft(F); dft(G); for(i=0;i<t;++i) F[i]=(1LL*F[i]*G[i])%md;
dft(F); reverse(F+1,F+t); p=Pow(t,md-2);
for(i=m+1;i<=r;++i) Add(f[i],1LL*F[i-l-1]*p%md); Solve(m+1,r);
}
int main(){
scanf("%d",&n); --n; f[0]=1;
for(i=1;i<=n;++i) scanf("%d",g+i); Solve(0,n);
for(i=0;i<=n;++i) printf("%d ",f[i]);
return 0;
}
多项式牛顿迭代
引入:导数
导数表
\((a)'=0\)(其中 \(a\) 表示和 \(x\) 无关的常数,下同)
\((x^a)'=ax^{a-1}\)
\((a^x)'=a^x\ln a\)
\((\log_ax)'=\frac 1{x\ln a}\)
\((\sin x)'=\cos x\)
\((\cos x)'=-\sin x\)
导数/积分运算法则
乘法法则:\((f(x)g(x))'=f'(x)g(x)+f(x)g'(x)\)
除法法则:\(\left(\frac{f(x)}{g(x)}\right)'=\frac{f'(x)g(x)-f(x)g'(x)}{g^2(x)}\)
链式法则:\(\frac{d}{dx}\frac{dx}{dy}=\frac{d}{dy}\)
复合函数求导:\(g(f(x))=g'(f(x))f'(x)\)
分部积分:\(fg=\int f'gdx+\int fg'dx\)
洛必达法则:如果 \(\lim_{x\to a}f(x)=0,\lim_{x\to a}g(x)=0\),且 \(f,g\) 在 \(a\) 的某去心邻域内两者均可导且 \(g'(x)\ne 0\),则有 \(\lim_{x\to a}\frac{f(x)}{g(x)}=\lim_{x\to a}\frac{f'(x)}{g'(x)}\)。如果 \(\lim_{x\to a}f(x)=\infty,\lim_{x\to a}g(x)=\infty\),且 \(f,g\) 在 \(a\) 的某去心邻域内两者均可导且 \(g'(x)\ne 0\),则有 \(\lim_{x\to a^+}\frac{f(x)}{g(x)}=\lim_{x\to a^+}\frac{f'(x)}{g'(x)}\)。
引入:泰勒展开
在逼*某个函数 \(f(x)\) 时,可以考虑使用某个多项式,通过对 \(f(x)\) 在某个位置 \(x_0\) 的 \(1\sim n\) 阶导数求出在同样位置处 \(1\sim n\) 阶导数相同的多项式,从而逼* \(f\),由 \((\sum_{i=0}^na_ix^i)'=\sum_{i=1}^nia_ix^{i-1}\) 有 \(f(x)=\sum_{i\ge 0}\frac{f^{(i)}(x_0)}{i!}(x-x_0)^i\)。某些时候可以直接取 \(x_0=0\),则有 \(f(x)=\sum_{i\ge 0}\frac{f^{(i)}(0)}{i!}x^i\)。一些例子。
引入:牛顿迭代
在求解方程 \(f(x)=0\) 的解时,我们考虑使用靠*零点的切线逼*解。具体地,取函数的某个位置 \(x_0\),得出在 \(x_0\) 位置和 \(f(x)\) 相切的切线 \(y=f'(x_0)(x-x_0)+f(x_0)\)。此时取该切线的零点 \(x_1=x_0-\frac{f(x_0)}{f'(x_0)}\),然后再求出 \(x_1\) 位置和 \(f(x)\) 相切的切线 \(y=f'(x_1)(x-x_1)+f(x_1)\);就可以不断迭代下去以逼*该解(显然无法得到解的精确值)。
虽然某些时候不能用牛顿迭代法逼*解,但是这个算法有较为广泛的应用,例如 \(O(1)\) 求某个数的*似*方根。同时牛顿迭代法收敛速度(逼*正解速度)快,可以说每次迭代后有效数位均会上升一倍。
在已知多项式 \(g(x)\),求满足 \(g(f(x))\equiv 0\pmod {x^n}\) 的 \(f(x)\bmod x^n\) 时,同样可以考虑牛顿迭代法的思路。设已经求出了 \(\bmod x^{\lceil\frac n2\rceil}\) 部分的解 \(f_0(x)\),(注意需要单独求出常数项)则求 \(f(x)\) 时可以考虑将 \(g\) 在 \(f_0(x)\) 处泰勒展开得:
\]
然后考虑 \(f_0(x)\) 的最高次数为 \(\lceil\frac n2\rceil-1\),所以 \(f(x)-f_0(x)\) 的最低次数非零项的次数为 \(\lceil\frac n2\rceil\),则上式 \(\bmod x^n\) 等效于 \(g(f_0(x))+g'(f_0(x))(f(x)-f_0(x))\bmod x^n\),解得 \(f(x)=f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}\)。
多项式牛顿迭代可以用于求多项式经过某些运算后的值。例如求多项式 \(h(x)\) 的逆时可以构造 \(g(f(x))=\frac 1{f(x)}-h(x)\),开方时可以构造 \(g(f(x))=f^2(x)-h(x)\),求 \(\exp\) 时可以构造 \(g(f(x))=\ln f(x)-h(x)\)。
多项式求逆
构造 \(g(f(x))=\frac 1{f(x)}-h(x)\),在求出上述的 \(f_0(x)\) 后,由于 \(g'(f(x))=-\frac 1{f^2(x)}\),则有 \(f(x)=f_0(x)-\frac{\frac 1{f_0(x)}-h(x)}{-\frac 1{f_0^2(x)}}=2f_0(x)-f_0^2(x)h(x)\)。
代码
点此查看代码
#include <bits/stdc++.h>
using namespace std;
const int maxb=262150;
const int md=998244353,G=3;
int n,i,j,k,u,v,p,d,t,le;
int f[maxb],g[maxb],h[maxb],r[maxb],x[maxb];
inline int Pow(int d,int z){
int r=1;
do{
if(z&1) r=1LL*r*d%md;
d=1LL*d*d%md;
}while(z>>=1); return r;
}
inline void Add(int &x,int y){x-=((x+=y)>=md)*md;}
void dft(int *f){
for(i=1;i<t;++i) if(i<r[i]) swap(f[i],f[r[i]]);
for(i=1,le=2;le<=t;i=le,le<<=1){
for(j=0,v=Pow(G,(md-1)/le);j<t;j+=le){
for(k=0,u=1;k<i;++k,u=1LL*u*v%md){
p=1LL*f[i+j+k]*u%md;
Add(f[i+j+k]=f[j+k],md-p);
Add(f[j+k],p);
}
}
}
}
int main(){
scanf("%d",&n);
for(i=0;i<n;++i) scanf("%d",h+i); n<<=1;
f[0]=Pow(h[0],md-2); Add(x[0]=f[0],f[0]);
for(d=2,t=4;d<=n;d=t,t<<=1){
for(i=0;i<t;++i) r[i]=(r[i>>1]>>1)+(i&1)*d;
memcpy(g,h,d<<2); dft(g); dft(f);
for(i=0;i<t;++i) f[i]=1LL*f[i]*f[i]%md*g[i]%md;
dft(f); reverse(f+1,f+t); p=Pow(t,md-2);
for(i=0;i<d;++i){
Add(f[i]=md-1LL*f[i]*p%md,x[i]);
Add(x[i]=f[i],f[i]); f[i+d]=0;
}
}
for(n>>=1,i=0;i<n;++i) printf("%d ",f[i]);
return 0;
}
多项式 \(\ln\)
考虑 \(\ln'(x)=\frac 1x\)。此时我们需要求 \(g(x)=\ln f(x)\),两边同时求导有 \(g'(x)=(\ln f(x))'=\frac{f'(x)}{f(x)}\)。注意此时 \(g(x)=\int_0^x\frac{f'(t)}{f(t)}dt+\ln f(0)\),而 \(e\) 为超越数(不能作为某个方程 \(\sum_i a_ix^i=0\) 的根,其中 \(a_i\) 均为整数,所以不能在模意义下被表示),所以 \(\ln f(0)\) 在 \(f(0)\) 也就是 \(a_0\) 不为 \(1\) 时不能在模意义下被表示。
代码
点此查看代码
#include <bits/stdc++.h>
using namespace std;
const int maxb=262150;
const int md=998244353,G=3;
int n,i,j,k,u,v,d,t,p,le;
int f[maxb],g[maxb],h[maxb],x[maxb],r[maxb];
inline int Pow(int d,int z){
int r=1;
do{
if(z&1) r=1LL*r*d%md;
d=1LL*d*d%md;
}while(z>>=1); return r;
}
inline void Add(int &x,int y){x-=((x+=y)>=md)*md;}
void dft(int *f){
for(i=1;i<t;++i) if(i<r[i]) swap(f[i],f[r[i]]);
for(i=1,le=2;le<=t;i=le,le<<=1){
for(j=0,v=Pow(G,(md-1)/le);j<t;j+=le){
for(k=0,u=1;k<i;++k,u=1LL*u*v%md){
p=1LL*f[i+j+k]*u%md;
Add(f[i+j+k]=f[j+k],md-p);
Add(f[j+k],p);
}
}
}
}
int main(){
scanf("%d",&n);
for(i=0;i<n;++i) scanf("%d",h+i); n<<=1;
for(f[0]=1,d=x[0]=2,t=4;d<=n;d=t,t<<=1){
for(i=0;i<t;++i) r[i]=(r[i>>1]>>1)+(i&1)*d;
memcpy(g,h,d<<2); dft(f); dft(g);
for(i=0;i<t;++i) f[i]=1LL*f[i]*f[i]%md*g[i]%md;
dft(f); reverse(f+1,f+t); p=Pow(t,md-2);
for(i=0;i<d;++i){
Add(f[i]=md-1LL*f[i]*p%md,x[i]);
Add(x[i]=f[i],f[i]); f[i+d]=0;
}
}
for(n>>=1,i=1;i<n;++i) h[i-1]=1LL*h[i]*i%md; h[n]=0;
memset(f+n,0,(d-n)<<2); for(t=1;t<=(n-1)<<1;t<<=1);
for(i=0;i<t;++i) r[i]=(r[i>>1]>>1)+(i&1)*(t>>1);
dft(f); dft(h); for(i=0;i<t;++i) f[i]=1LL*f[i]*h[i]%md;
dft(f); reverse(f+1,f+t); p=Pow(t,md-2);
for(i=0;i<t;++i) f[i]=1LL*f[i]*p%md;
for(i=n-1;i;--i) f[i]=1LL*f[i-1]*Pow(i,md-2)%md;
for(f[0]=i=0;i<n;++i) printf("%d ",f[i]); return 0;
}
多项式 \(\exp\)
构造 \(g(f(x))=\ln f(x)-h(x)\),在求出上述的 \(f_0(x)\) 后,由于 \(g'(f(x))=\frac1{f(x)}\),则有 \(f(x)=f_0(x)-\frac{\ln f_0(x)-h(x)}{\frac1{f_0(x)}}=f_0(x)(1-\ln f_0(x)+h(x))\)。
代码
点此查看代码
#include <bits/stdc++.h>
using namespace std;
const int maxb=262150;
const int maxt=524300;
const int md=998244353,G=3;
int n,i,j,k,u,v,p,t,le,d0,t0,d1;
int f0[maxb],f1[maxt],g1[maxb],g[maxt],h[maxb],x[maxt],r[maxt];
inline int Pow(int d,int z){
int r=1;
do{
if(z&1) r=1LL*r*d%md;
d=1LL*d*d%md;
}while(z>>=1); return r;
}
inline void Add(int &x,int y){x-=((x+=y)>=md)*md;}
void dft(int *f){
for(i=1;i<t;++i) if(i<r[i]) swap(f[i],f[r[i]]);
for(i=1,le=2;le<=t;i=le,le<<=1){
for(j=0,v=Pow(G,(md-1)/le);j<t;j+=le){
for(k=0,u=1;k<i;++k,u=1LL*u*v%md){
p=1LL*f[i+j+k]*u%md;
Add(f[i+j+k]=f[j+k],md-p);
Add(f[j+k],p);
}
}
}
}
int main(){
scanf("%d",&n);
for(i=1;i<n;++i) scanf("%d",h+i); n<<=1;
for(f0[0]=1,d0=2,t0=4;d0<=n;d0=t0,t0<<=1){
for(i=1;i<d0;++i) g1[i-1]=1LL*f0[i]*i%md;
for(f1[0]=1,x[0]=d1=2,t=4;d1<=t0;d1=t,t<<=1){
for(i=1;i<t;++i) r[i]=(r[i>>1]>>1)|((i&1)*d1);
memcpy(g,f0,d1<<2); dft(f1); dft(g);
for(i=0;i<t;++i) f1[i]=1LL*f1[i]*f1[i]%md*g[i]%md;
dft(f1); reverse(f1+1,f1+t); p=Pow(t,md-2);
for(i=0;i<d1;++i){
Add(f1[i]=md-1LL*f1[i]*p%md,x[i]);
Add(x[i]=f1[i],f1[i]); f1[i+d1]=0;
}
}
memset(f1+d0,0,(d1-d0)<<2); for(t=1;t<=t0-2;t<<=1);
for(i=1;i<t;++i) r[i]=(r[i>>1]>>1)|((i&1)*(t>>1));
dft(f1); dft(g1); for(i=0;i<t;++i) f1[i]=1LL*f1[i]*g1[i]%md;
dft(f1); reverse(f1+1,f1+t); p=Pow(t,md-2);
memset(f1+d0,0,(t-d0)<<2); memset(g1,0,t<<2);
for(i=0;i<d0;++i) f1[i]=1LL*f1[i]*p%md;
for(i=d0-1;i;--i) f1[i]=1LL*f1[i-1]*Pow(i,md-2)%md;
for(f1[0]=i=1;i<d0;++i) f1[i]=(1LL*f1[i]*(md-1)+h[i])%md;
for(i=1,t=t0;i<t;++i) r[i]=(r[i>>1]>>1)|((i&1)*d0);
dft(f0); dft(f1); for(i=0;i<t;++i) f0[i]=1LL*f0[i]*f1[i]%md;
dft(f0); reverse(f0+1,f0+t); p=Pow(t,md-2);
for(i=0;i<d0;++i) f0[i]=1LL*f0[i]*p%md;
memset(f0+d0,0,d0<<2); memset(f1,0,t<<2); memset(g,0,d1<<2);
}
for(n>>=1,i=0;i<n;++i) printf("%d ",f0[i]); return 0;
}
多项式带余除法(待补)
多项式多点求值(待补)
拉格朗日插值(待补)
FMT/FWT(一种清奇(愚蠢?)的思路)
原文链接: https://www.cnblogs.com/Fran-CENSORED-Cwoi/p/17137654.html
欢迎关注
微信关注下方公众号,第一时间获取干货硬货;公众号内回复【pdf】免费获取数百本计算机经典书籍
原创文章受到原创版权保护。转载请注明出处:https://www.ccppcoding.com/archives/319404
非原创文章文中已经注明原地址,如有侵权,联系删除
关注公众号【高性能架构探索】,第一时间获取最新文章
转载文章受原作者版权保护。转载请注明原作者出处!