多项式
Posted lhm-
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了多项式相关的知识，希望对你有一定的参考价值。
形如 (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))的得到它们乘积的点值表示法
从多项式的系数表示向点值表示的转化，称为求值，从多项式的点值表示向系数表示的转化，称为插值
(FFT)
快速傅里叶变换可以做到(O(nlog n))进行求值和插值
复数相乘，模长相乘，幅角相加
((a+bi)×(c+di)=(ac-bd)+(bc+ad)i)
在复平面上，以原点为圆心，(1)为半径作圆，所得的圆叫单位圆。以圆点为起点，圆的(n)等分点为终点，做(n)个向量，设幅角为正且最小的向量对应的复数为(ω_n)，称为(n)次单位根
(ω_n^k=cosfrac{2πk}{n}+isinfrac{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));(NTT)
若(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} ceil}))
由原式得
(f(x)g (x)≡1 (mod x^{lceil frac{n}{2} ceil}))
得
(f(x)(g^prime (x)-g(x))≡0 (mod x^{lceil frac{n}{2} ceil}))
同除(f(x))
(g^prime (x)-g(x)≡0 (mod x^{lceil frac{n}{2} ceil}))
平方
(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]);多项式 (ln)
(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]);以上是关于多项式的主要内容，如果未能解决你的问题，请参考以下文章
数字信号处理卷积编程实现 ( Matlab 卷积和多项式乘法 conv 函数 | 使用 matlab 代码求卷积并绘图 )