NTT(快速数论变换）多项式求逆一千五百字解析

qinsanma and Code

已于 2022-07-24 22:48:40 修改

阅读量313

点赞数 1

分类专栏：数论进阶

于 2022-07-24 22:45:08 首次发布

本文链接：https://blog.csdn.net/jisuanji2606414/article/details/125963430

版权

快速幂运算快速傅里叶变换多项式计算 C++ 算法

关键词由CSDN通过智能技术生成

数论进阶专栏收录该内容

17 篇文章 7 订阅

订阅专栏

首先我们有 $\tiny F(x)G(x)\equiv 1(mod \,x^n)$

我们对x^n进行向上取整开根号操作，即 $\tiny x^n\rightarrow x^{\left \lceil n/2 \right \rceil}$

我们设 $\tiny F(x)H(x)\equiv 1(mod \, x^{\left \lceil \frac{n}{2} \right \rceil})$

显然 $\tiny F(x)G(x)\equiv 1 (mod\, x^{\left \lceil \frac{n}{2} \right \rceil})$

两式相减

$\tiny F(x)(G(x)-H(x)) \equiv 0 (mod \, x^{\left \lceil \frac{n}{2} \right \rceil})$

故 $\tiny G(x)-H(x)\equiv 0 (mod \,x^{\left \lceil \frac{n}{2} \right \rceil })$

所以减式构成的多项式在0- $\tiny \left \lceil \frac{n}{2} \right \rceil-1$ 次幂上都是0，否则他们本身就无法被 $\tiny x^{\left \lceil \frac{n}{2} \right \rceil}$ 整除

根据这一发现

对减式进行平方

$\tiny G(x)^2 +H(x)^2 -2G(x)H(x)\equiv 0 (mod \, x^n)$

$\tiny F(x)G(x)\equiv 1(mod \, x^n)$

进行乘F(x)操作再由

$\tiny G(x)\equiv 2H(x)-F(x)H(x)^2 (mod \, x^n)$

由此G(x)用NTT求解

# include<iostream>


using namespace std;
typedef long long int  ll;

const int mod=998244353,g=3,N=2000000;

int n;
ll a[N],b[N],c[N],rev[N];

ll qp(ll base,ll pow)
{
    ll ans=1;

    while(pow)
    {
        if(pow&1)
            ans=ans*base%mod;

        pow>>=1;

        base=base*base%mod;

    }

    return ans;

}

void init(int k)
{
    int len=(1<<k);

    for(int i=0;i<len;i++)
    {
        rev[i]=0;

    }

    for(int i=0;i<len;i++)
    {

        rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));

    }
    return  ;

}

void NTT(ll *a,int n,int flag)
{

    for(int i=0;i<n;i++)
    {
        if(i<rev[i])
        {
            swap(a[i],a[rev[i]]);

        }
    }

    for(int h=1;h<n;h<<=1)
    {
        ll gn=qp(3ll,(mod-1)/(h*2));


        if(flag==-1)
            gn=qp(gn,mod-2);

        for(int j=0;j<n;j+=2*h)
        {
            ll gg=1;

            for(int k=j;k<j+h;k++)
            {
                ll x=a[k];

                ll y=gg*a[k+h]%mod;

                a[k]=(x+y)%mod;

                a[k+h]=((x-y)%mod+mod)%mod;

                gg=gg*gn%mod;


            }
        }
    }

    if(flag==-1)
    {
        ll inv=qp(n,mod-2);

        for(int i=0;i<n;i++)
        {
            a[i]=a[i]*inv%mod;

        }
    }
}
void work(int pow,ll*a,ll *b)
{
    if(pow==1)
    {
        b[0]=qp(a[0],mod-2);

        return ;

    }

    work((pow+1)>>1,a,b);

    int k=1,s=2;

    while((1<<k)<pow+pow-1)
    {
        k++;

        s<<=1;

    }

    init(k);

    for(int i=0;i<pow;i++)    
    {
        c[i]=a[i];
    }
   
 /*
之所以用c存a，而不对a进行直接操作是因为我们每次递归的
都必须是原数列的a
*/

    NTT(c,s,1);

    NTT(b,s,1);

    for(int i=0;i<s;i++)
    {
        b[i]=(ll)(2ll-c[i]*b[i]%mod+mod)%mod*b[i]%mod;
    }

    NTT(b,s,-1);

    for(int i=pow;i<s;i++)
    {
        b[i]=0;
    }

    return ;


}

int main ()
{

    cin>>n;

    for(int i=0;i<n;i++)
    {
        cin>>a[i];

    }

    work(n,a,b);

    for(int i=0;i<n;i++)
    {
        cout<<b[i]<<" ";

    }

    return 0;

}