标签:operator type pre line com inter while 向量 反转
形如 \(f(x)=a_0x^0+a_1x^1+a_2x^2+...+a_{n-1}x^{n-1}\)
点值表示法:通过代入\(n\)个不同的值\(x_0,x_1...x_{n-1}\)到\(f(x)\)中,得到\(y_0,y_1...y_{n-1}\),用\((x_0,y_0),(x_1,y_1)...(x_{n-1},y_{n-1})\)即可表示这个多项式
若用点值表示法,并且两个多项式的\(x\)对应相等,那么将\(y\)对应相乘,即可\(O(n)\)的得到它们乘积的点值表示法
从多项式的系数表示向点值表示的转化,称为求值,从多项式的点值表示向系数表示的转化,称为插值
快速傅里叶变换可以做到\(O(nlog\ n)\)进行求值和插值
复数相乘,模长相乘,幅角相加
\((a+bi)×(c+di)=(ac-bd)+(bc+ad)i\)
在复平面上,以原点为圆心,\(1\)为半径作圆,所得的圆叫单位圆。以圆点为起点,圆的\(n\)等分点为终点,做\(n\)个向量,设幅角为正且最小的向量对应的复数为\(ω_n\),称为\(n\)次单位根
\(ω_n^k=cos\frac{2πk}{n}+isin\frac{2πk}{n}\)
\(ω_{2n}^{2k}=ω_n^k\)
\(ω_{n}^{k+\frac{n}{2}}=-ω_n^k\)
\(f(x)=a_0x^0+a_1x^1+a_2x^2+a_3x^3+...+a_{n-2}x^{n-2}+a_{n-1}x^{n-1}\)
进行奇偶分类得
\(f1(x)=a_0+a_2x^1+a_4x^2+a_6x^3+...+a_{n-2}x^{\frac{n}{2}-1}\)
\(f2(x)=a_1+a_3x^1+a_5x^2+a_7x^3+...+a_{n-1}x^{\frac{n}{2}-1}\)
得\(f(x)=f1(x^2)+xf2(x^2)\)
设\(k<\frac{n}{2}\),代入\(ω_n^k\)得
\(f(ω_n^k)=f1(ω_n^{2k})+ω_n^kf2(ω_n^{2k})\)
再代入\(ω_n^{k+\frac{n}{2}}\)得
\(f(ω_n^{k+\frac{n}{2}})=f1(ω_n^{2k})-ω_n^kf2(ω_n^{2k})\)
通过这两个式子,我们就可以进行递归了,但应递归常数过大,我们通过迭代来实现
我们将原多项式系数重新排序,将每个系数都在它递归的最后位置,就可以用迭代来代替递归实现
发现原来在第\(a\)个位置的系数,在递归的最后位置为\(a\)的二进制反转以后的数
这称为蝴蝶操作
代入\(ω_n^{-k}\),可再次求得一系列值,将其除以\(n\)即可求得多项式相乘后的系数表示
\(code :\)
struct Complex
{
double x,y;
Complex(double a=0,double b=0)
{
x=a,y=b;
}
}A[maxn],B[maxn];
Complex operator +(const Complex &a,const Complex &b)
{
return Complex(a.x+b.x,a.y+b.y);
}
Complex operator -(const Complex &a,const Complex &b)
{
return Complex(a.x-b.x,a.y-b.y);
}
Complex operator *(const Complex &a,const Complex &b)
{
return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
void FFT(Complex *a,int type)
{
for(int i=0;i<lim;++i)
if(i<r[i])
swap(a[i],a[r[i]]);
for(int mid=1;mid<lim;mid<<=1)
{
Complex wn(cos(Pi/mid),type*sin(Pi/mid));
for(int j=0;j<lim;j+=(mid<<1))
{
Complex w(1,0);
for(int k=0;k<mid;++k,w=w*wn)
{
Complex x=a[j+k],y=t*a[j+k+mid];
a[j+k]=x+y,a[j+k+mid]=x-y;
}
}
}
}
......
while(lim<=n+m) lim<<=1,l++;
for(int i=0;i<lim;++i)
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
FFT(A,1),FFT(B,1);
for(int i=0;i<=lim;++i) A[i]=A[i]*B[i];
FFT(A,-1);
for(int i=0;i<=n+m;++i)
printf("%d ",(int)(A[i].x/lim+0.5));
若\(a,p\)互素,且\(p>1\),对于\(a^n≡1\ (mod\ p)\)最小的\(n\),称为\(a\)模\(p\)的阶,记作\(δ_p(a)\)
设\(p\)是正整数,\(a\)是整数,若\(δ_p(a)=?(p)\),则称\(a\)为模\(p\)的一个原根
原根个数不唯一
若\(P\)为素数,设\(G\)为\(P\)的原根,那么\(G^i\ mod\ P\ (i<G<P,\ 0<i<P)\)的结果两两不同
模数有\(998244353,1004535809,469762049\),原根都为\(3\)
\(ω_n≡g^{\frac{p-1}{n}}\ mod\ p\)
\(code :\)
#define P 998244353
#define G 3
#define Gi (P+1)/G
ll qp(ll x,ll y)
{
ll ans=1;
while(y)
{
if(y&1) ans=(ans*x)%P;
x=(x*x)%P;
y>>=1;
}
return ans;
}
void NTT(ll *a,int type)
{
for(int i=0;i<lim;++i)
if(i<r[i])
swap(a[i],a[r[i]]);
for(int mid=1;mid<lim;mid<<=1)
{
ll wn=qp(type==1?G:Gi,(P-1)/(mid<<1));
for(int j=0;j<lim;j+=(mid<<1))
{
ll w=1;
for(int k=0;k<mid;++k,w=(w*wn)%P)
{
int x=a[j+k],y=(w*a[j+k+mid])%P;
a[j+k]=(x+y)%P,a[j+k+mid]=(x-y+P)%P;
}
}
}
}
......
while(lim<=n+m) lim<<=1,l++;
for(int i=0;i<lim;++i)
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
NTT(A,1),NTT(B,1);
for(int i=0;i<=lim;++i) A[i]=(A[i]*B[i])%P;
NTT(A,-1);
inv=qp(lim,P-2);
for(int i=0;i<=n+m;++i)
printf("%lld ",(A[i]*inv)%P);
\(f(x)g(x)≡1\ (mod\ x^n)\)
称\(g(x)\)为\(f(x)\)的逆元,模\(x^n\)的意义是将\(\geqslant\ n\)的项都忽略掉
设
\(f(x)g^\prime (x)≡1\ (mod\ x^{\lceil \frac{n}{2}\rceil})\)
由原式得
\(f(x)g (x)≡1\ (mod\ x^{\lceil \frac{n}{2}\rceil})\)
得
\(f(x)(g^\prime (x)-g(x))≡0\ (mod\ x^{\lceil \frac{n}{2}\rceil})\)
同除\(f(x)\)
\(g^\prime (x)-g(x)≡0\ (mod\ x^{\lceil \frac{n}{2}\rceil})\)
平方
\(g^{\prime2} (x)+g^2(x)-2g^\prime (x)g(x)≡0\ (mod\ x^n)\)
同乘\(f(x)\)
\(g^{\prime2} (x)f(x)+g(x)-2g^\prime (x)≡0\ (mod\ x^n)\)
得
\(g(x)≡g^\prime (x)(2-g^\prime (x)f(x))\ (mod\ x^n)\)
\(code :\)
ll qp(ll x,ll y)
{
ll ans=1;
while(y)
{
if(y&1) ans=(ans*x)%P;
x=(x*x)%P;
y>>=1;
}
return ans;
}
void NTT(ll *a,int type)
{
for(int i=0;i<lim;++i)
if(i<r[i])
swap(a[i],a[r[i]]);
for(int mid=1;mid<lim;mid<<=1)
{
ll wn=qp(type==1?G:Gi,(P-1)/(mid<<1));
for(int j=0;j<lim;j+=(mid<<1))
{
ll w=1;
for(int k=0;k<mid;++k,w=(w*wn)%P)
{
ll x=a[j+k],y=w*a[j+k+mid]%P;
a[j+k]=(x+y)%P,a[j+k+mid]=(x-y+P)%P;
}
}
}
if(type==1) return;
ll inv=qp(lim,P-2);
for(int i=0;i<lim;++i)
a[i]=a[i]*inv%P;
}
void work(int deg,ll *a,ll *b)
{
if(deg==1)
{
b[0]=qp(a[0],P-2);
return;
}
work((deg+1)>>1,a,b);
lim=1,l=0;
while(lim<(deg<<1)) lim<<=1,l++;
for(int i=0;i<lim;++i)
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
for(int i=0;i<deg;++i) C[i]=a[i];
for(int i=deg;i<lim;++i) C[i]=0;
NTT(C,1),NTT(b,1);
for(int i=0;i<lim;++i)
b[i]=((ll)2-C[i]*b[i]%P+P)%P*b[i]%P;
NTT(b,-1);
for(int i=deg;i<lim;++i) b[i]=0;
}
......
work(n,A,B);
for(int i=0;i<n;++i) printf("%lld ",B[i]);
\(code :\)
ll qp(ll x,ll y)
{
ll ans=1;
while(y)
{
if(y&1) ans=(ans*x)%P;
x=(x*x)%P;
y>>=1;
}
return ans;
}
void NTT(ll *a,int type)
{
for(int i=0;i<lim;++i)
if(i<r[i])
swap(a[i],a[r[i]]);
for(int mid=1;mid<lim;mid<<=1)
{
ll wn=qp(type==1?G:Gi,(P-1)/(mid<<1));
for(int j=0;j<lim;j+=(mid<<1))
{
ll w=1;
for(int k=0;k<mid;++k,w=(w*wn)%P)
{
ll x=a[j+k],y=w*a[j+k+mid]%P;
a[j+k]=(x+y)%P,a[j+k+mid]=(x-y+P)%P;
}
}
}
if(type==1) return;
ll inv=qp(lim,P-2);
for(int i=0;i<lim;++i)
a[i]=a[i]*inv%P;
}
void work(int deg,ll *a,ll *b)
{
if(deg==1)
{
b[0]=qp(a[0],P-2);
return;
}
work((deg+1)>>1,a,b);
lim=1,l=0;
while(lim<(deg<<1)) lim<<=1,l++;
for(int i=0;i<lim;++i)
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
for(int i=0;i<deg;++i) C[i]=a[i];
for(int i=deg;i<lim;++i) C[i]=0;
NTT(C,1),NTT(b,1);
for(int i=0;i<lim;++i)
b[i]=((ll)2-C[i]*b[i]%P+P)%P*b[i]%P;
NTT(b,-1);
for(int i=deg;i<lim;++i) b[i]=0;
}
void direv()
{
for(int i=1;i<n;++i) A[i-1]=A[i]*i%P;
A[n-1]=0;
}
void inter()
{
for(int i=n-1;i;--i) A[i]=A[i-1]*qp(i,P-2)%P;
A[0]=0;
}
......
work(n,A,B);
direv();
lim=1,l=0;
while(lim<=n*2) lim<<=1,l++;
NTT(A,1),NTT(B,1);
for(int i=0;i<lim;++i) A[i]=(A[i]*B[i])%P;
NTT(A,-1);
inter();
for(int i=0;i<n;++i) printf("%lld ",A[i]);
标签:operator type pre line com inter while 向量 反转
原文地址:https://www.cnblogs.com/lhm-/p/12229547.html