码迷,mamicode.com
首页 > 其他好文 > 详细

多项式

时间:2020-01-22 22:00:21      阅读:62      评论:0      收藏:0      [点我收藏+]

标签: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)\)的得到它们乘积的点值表示法

从多项式的系数表示向点值表示的转化,称为求值,从多项式的点值表示向系数表示的转化,称为插值

\(FFT\)

快速傅里叶变换可以做到\(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));

\(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}\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]);

多项式 \(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]);

多项式

标签:operator   type   pre   line   com   inter   while   向量   反转   

原文地址:https://www.cnblogs.com/lhm-/p/12229547.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!