投稿
  1. 高性能架构探索首页
  2. C/C++

多项式求逆及其应用

admin • • C/C++ • 阅读 160

$part ~ 1 ~ $多项式求逆（乘法逆元QwQ）

01 题目描述

给定一个多项式 \(F(x)\) ，请求出一个多项式 \(G(x)\)， 满足 \(F(x) * G(x) \equiv 1 \pmod{x^n}\)。系数对 \(998244353\) 取模。

输入格式

首先输入一个整数 \(n\)， 表示输入多项式的项数。
接着输入 \(n\) 个整数，第 \(i\) 个整数 \(a_i\) 代表 \(F(x)\) 次数为 \(i-1\) 的项的系数。

输出格式

输出 \(n\) 个数字，第 \(i\) 个整数 \(b_i\) 代表 \(G(x)\) 次数为 \(i-1\) 的项的系数。保证有解。

样例 #1
样例输入 #1
5
1 6 3 4 9
样例输出 #1
1 998244347 33 998244169 1020
提示

对于 \(100\%\) 的数据，\(1 \leq n \leq 10^5\)，$ 0 \leq a_i \leq 10^9$。

02 问题化归

显然，如果我们的问题是求解满足\(F(x)*G(x)\equiv 1\pmod x\)的\(G(x)\)的话肯定都会解，最简单的想法就是求F(x)的常数项在模998244353意义下的逆元，以此作为G(x)的逆元，G(x)其余的项就全部随意啦或者说直接令其为零即可

但是我们的题目现在是求解满足\(F(x)*G(x)\equiv 1\pmod {x^n}\)的G(x)，那么一个很自然的思路就是把我们的\(x^n\)给化归到\(x\)在放回去解决，而在程序实现中我们最通常用的手段是每次除以二的化归

03 化归关系式推导

首先我们可以有这么一个结论：当\(m\leqslant n\)时\(\pmod {x^n}\)肯定强于\(\pmod {x^m}\)，因此我们假设\(F(x)\)在\(\pmod{x^n}\)意义下的逆元是\(G_n(x)\)，那么\(G_n(x)\)肯定也会是\(\pmod{x^{\left\lceil\frac{n}2\right\rceil}}\)意义下的逆元（方便起见，本部分后文中无特殊标明\(\dfrac{n}{2}\)默认为向上取整），即

\[G_n(x)\equiv G_{\frac{n}{2}}(x)\pmod{x^{\frac{n}{2}}}
\]

移项可得

\[G_n(x)-G_{\frac{n}2}(x)\equiv 0\pmod {x^{\frac{n}2}}
\]

平方

\[G^2(x)-2G(x)G_{\frac{n}2}(x)+G_{\frac{n}2}^2(x)\equiv0\pmod{x^n}
\]

于是在这里面我们为了尽量消去未知的\(G(x)\)，我们在两边同时乘以\(F(x)\)，则

\[G(x)-2G_{\frac{n}2}(x)+F(x)G_{\frac{n}2}^2(x)\equiv0\pmod{x^n}
\]

移项便可得到递推式

\[G(x)\equiv G_{\frac{n}2}(x)(2-F(x)G_{\frac{n}2}(x))\pmod{x^n}
\]

04 程序实现

利用好点值转化的关系，求出\(F(x)\)与\(G_{\frac{n}2}(x)\)的点值表示之后直接用这个递推式得到我们\(G(x)\)的点值表示，之后再将点值表示转化回系数表达可一定程度上减少我们利用\(fft/ntt\)加速多项式乘法的常数

另外关于这道题因为我们已经具体要求了998244353的模数，因此直接用\(ntt\)处理会比\(fft\)优秀便捷许多：

void solve(int len,ll *f,ll *g)
{
	if(len==1)
    {
        g[0]=power(f[0],mod-2);
        return;
    }
	solve((len+1)>>1,f,g,n);
	init(len<<1);
	for(register int i=0;i<sup;i++)
    {
        a[i]=i<len?f[i]:0;
        b[i]=(i<(len+1)>>1?g[i]:0);
	}
    ntt(a,sup,1);ntt(b,sup,1);
	for(register int i=0;i<sup;i++)
		a[i]=b[i]*((2-a[i]*b[i]%mod+mod)%mod)%mod;
	ntt(a,sup,0);
	long long invsup=power(sup,mod-2);
	for(register int i=0;i<len;i++)
		g[i]=a[i]*invup%mod;
}

当然，对于追求时间优化而言，我们的原则是能不递归就不递归，因此我们可以先循环得到我们所需要的递归序然后进行一次正向递推，这样可以一定程度上优化我们的代码：

int main()
{
    int n;
    scanf("%d",&n);
    for(register int i=0;i<n;++i) scanf("%lld",&f[i]);
    F[0]=power(f[0],mod-2);
    for(;n>1;n=(n+1)>>1)
        ite[++ite[0]]=n;
    ite[ite[0]+1]=1;
    for(register int k=ite[0];k;--k)
    {
        init(ite[k]<<1);
        for(register int i=0;i<sup;++i)
        {
            A[i]=i<ite[k]?f[i]:0;
            G[i]=i<ite[k+1]?F[i]:0;
        }
        ntt(A,sup,1),ntt(G,sup,1);
        for(register int i=0;i<sup;++i)
            A[i]=G[i]*((2-A[i]*G[i]%mod+mod)%mod)%mod;
        ntt(A,sup,0);
        for(register int i=0,j=power(sup,mod-2);i<ite[k];++i)
            F[i]=A[i]*j%mod;
    }
	for(int i=0;i<ite[1];i++) printf("%lld ",F[i]);
    putchar('\n');
    return 0;
}

而时间复杂度满足如下递推式：

\[T(n)=T(\left\lceil\frac{n}{2}\right\rceil)+\Theta(n\log_2n)
\]

而显然我们会有\(n\log_2n\geqslant 2\times\dfrac{n}{2}\log_2{\dfrac{n}{2}}\)，借此我们可以得到

\[1\log_2 1+2\log_2 2+4\log_2 4+\cdots+n\log_2n\leqslant 2n\log_2n
\]

因此，我们的时间复杂度由累加法递推可以得到

\[\Theta(n\log_2n)<T(n)<\Theta(2n\log_2n)
\]

统合一下就是：\(T(n)=O(n\log_2n)\)

05 特殊的优化

前面讲到我们将\(\pmod{x^n}\)递推化归为\(\pmod{x^{\frac{n}2}}\)是因为我们在计算机程序设计中经常涉及到二进制的运算，这样的二分是我们最熟悉的

然而现在我们可以再跳出目前我们惯性的二分想法，能不能得到\(\pmod{x^n}\)化归为\(\pmod{x^{\left\lceil\frac{n}{3}\right\rceil}}\)乃至\(\pmod{x^{\left\lceil\frac{n}{m}\right\rceil}}\)的递推式呢？

答案是可以的，以\(\pmod{x^{\left\lceil\frac{n}{3}\right\rceil}}\)为例，我们可以以同样的思路推一波式子（同样以下的\(\dfrac{n}3\)均理解为\(\left\lceil\dfrac{n}3\right\rceil\)）：

\[G_n(x)-G_{\frac{n}3}(x)\equiv 0\pmod {x^{\frac{n}3}}
\]

立方

\[G^3(x)-3G^2(x)G_{\frac{n}3}(x)+3G(x)G_{\frac{n}3}^2(x)-G^3_{\frac{n}3}(x)\equiv0\pmod{x^n}
\]

于是在这里面我们为了尽量消去未知的\(G(x)\)，我们在两边同时乘以\(F^2(x)\)，之后便可移项得到递推式

\[G(x)\equiv G_{\frac{n}3}(x)(3-3F(x)G_{\frac{n}3}(x)+F^2(x)G^2_{\frac{n}3}(x))\pmod{x^n}
\]

而对于这个递推我们发现复杂度满足如下递推式：

\[T(n)=T(\left\lceil\frac{n}{3}\right\rceil)+\Theta(2n\log_22n)
\]

可是对这个我们就不是很好分析了，具体来说的话其实应该跟二分思路的复杂度差不多，如果同样是递归的话由于递归层数偏少可能会快一些，亲测是发现用递推的话三分比二分略慢

int main()
{
    int n;
    scanf("%d",&n);
    for(register int i=0;i<n;++i) scanf("%lld",&f[i]);
    F[0]=power(f[0],mod-2);
    for(;n>1;n=ceil(n/3.0)) ite[++ite[0]]=n;
    ite[ite[0]+1]=1;
    for(register int k=ite[0];k;--k)
    {
        init(ite[k]<<2);
        for(register int i=0;i<sup;++i)
        {
            A[i]=i<ite[k]?f[i]:0;
            G[i]=i<ite[k+1]?F[i]:0;
        }
        ntt(A,sup,1),ntt(G,sup,1);
        for(register int i=0;i<sup;++i)
            A[i]=G[i]*((3ll-(3ll*A[i]%mod)*G[i]%mod+((((A[i]*A[i])%mod)*G[i])%mod)*G[i]%mod+mod)%mod)%mod;
        ntt(A,sup,0);
        for(register int i=0,j=power(sup,mod-2);i<ite[k];++i)
            F[i]=A[i]*j%mod;
    }
	for(int i=0;i<ite[1];i++) printf("%lld ",F[i]);
    putchar('\n');
    return 0;
}

总结个一般的情况：

\[G_n(x)-G_{\frac nm}(x)\equiv 0\pmod {x^{\frac{n}m}}\Rightarrow(G_n(x)-G_{\frac nm}(x))^m\equiv 0\pmod {x^n}
\]

则

\[\sum\limits_{i=0}^m(-1)^iC_m^iG^{m-i}(x)G^i_{\frac nm}(x)\equiv0\pmod{x^n}
\]

于是

\[G(x)+\sum\limits_{i=1}^m(-1)^iC_m^iF^{i-1}(x)G^i_{\frac nm}(x)\equiv0\pmod{x^n}
\]

得到

\[G(x)\equiv-\sum\limits_{i=1}^m(-1)^iC_m^iF^{i-1}(x)G^i_{\frac nm}(x)\pmod{x^n}
\]

最特殊的时候，我们取m=n，那么

\[G(x)\equiv\sum\limits_{i=0}^{n-1}(-1)^iC_n^{i+1}Q^{i+1}F^i(x)\pmod{x^n}
\]

其中Q为F(x)常数项模998244353意义下的逆元，也就相当于\(G_1(x)\)，而这也就相当于我们最终的直接表达递推式，只是对于\(F^i(x)\)，我们不太好直接处理，会导致复杂度达到不尽如人意的\(O(n^2\log_2n)\). 而手玩一下样例便可通过数学归纳法知道，当\(m\geqslant3\)时，随着m的增大时间复杂度是单增的，所以比较以上最好的想法肯定还是m取2或者3

$part ~ 2 ~ $多项式开方

01 题目描述

给定一个 \(n-1\) 次多项式 \(A(x)\)，求一个在 \({} \bmod x^n\) 意义下的多项式 \(B(x)\)，使得 \(B^2(x) \equiv A(x) \pmod{x^n}\)。若有多解，请取零次项系数较小的作为答案。

多项式的系数在 \({}\bmod 998244353\) 的意义下进行运算。

输入格式

第一行一个正整数 \(n\)。

接下来 \(n\) 个整数，依次表示多项式的系数 \(a_0, a_1, \dots, a_{n-1}\)。

保证 \(a_0 = 1\)。

输出格式

输出 \(n\) 个整数，表示答案多项式的系数 \(b_0, b_1, \dots, b_{n-1}\)。

样例 #1
样例输入 #1
3
1 2 1
样例输出 #1
1 1 0
样例 #2
样例输入 #2
7
1 8596489 489489 4894 1564 489 35789489
样例输出 #2
1 503420421 924499237 13354513 217017417 707895465 411020414
提示

对于 \(100 \%\) 的数据：\(1 \le n \leq 10^5\)，\(0 \le a_i < 998244353\)。

02 化归

这是多项式求逆的一个简单应用，首先就是考验基本功的推式子

推式子的思路当然是与上面相同，同样还是先假设我们要求的多项式变为

\[B_n^2(x)\equiv A\pmod{x^n}
\]

那么我们会有

\[B_n(x)-B_{\frac n2}(x)\equiv0\pmod{x^{\frac n2}}
\]

平方

\[B_n^2(x)-2B_n(x)B_{\frac n2}(x)+B_{\frac n2}^2(x)\equiv A(x)-2B_n(x)B_{\frac n2}(x)+B_{\frac n2}^2(x)\equiv0\pmod{x^n}
\]

移项

\[2B_n(x)B_{\frac n2}(x)\equiv A(x)+B^2_{\frac n2}(x)\pmod{x^n}
\]

结合我们$part~ 1 ~ $中的多项式求逆，我们可以得到

\[B_n(x)\equiv2^{-1}B_{\frac n2}^{-1}(x)(A(x)+B^2_{\frac n2}(x))\pmod {x^n}
\]

也就是说，我们每次把之前得到的\(B_{\frac n2}(x)\)在\(\pmod{x^n}\)意义下的逆元得到然后就可以轻松地求解了

03 程序实现

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
const int g=3,_g=332748118,mod=998244353,xx=(1<<21)+10;
const ll _2=499122177;
int rev[xx],sup,all;

inline int read()
{
    int x=0;char ch=getchar();
    for(;!isalnum(ch);ch=getchar());
    for(;isalnum(ch);ch=getchar()) x=(x<<3)+(x<<1)+ch-'0';
    return x;
}
inline ll power(ll x,int y)
{
    ll ans=1;
    for(;y;y>>=1,(x*=x)%=mod)
        if(y&1) (ans*=x)%=mod;
    return ans;
}
inline void init(int n)
{
    for(sup=1,all=-1;sup<=n;sup<<=1,++all);
    for(register int i=0;i<sup;++i)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<all);
}
inline void ntt(ll *a,const bool op)
{
    for(register int i=0;i<sup;++i)
        if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(register int k=1;k<sup;k<<=1)
    {
        const ll g0=power(op?g:_g,(mod-1)/(k<<1));
        for(register int i=0;i<sup;i+=(k<<1))
        {
            ll gj=1;
            for(register int j=0;j<k;++j,(gj*=g0)%=mod)
            {
                const ll l=a[i+j],r=a[i+k+j]*gj%mod;
                a[i+j]=(l+r)%mod;
                a[i+k+j]=(l-r+mod)%mod;
            }
        }
    }
}
ll tmp1[xx],tmp2[xx],invsup;
inline void inv(ll *a,ll *b,int n)
{
    int ite[21];
    for(register int i=0;i<21;++i) ite[i]=0;
    for(;n>1;n=(n+1)>>1)
        ite[++ite[0]]=n;
    ite[ite[0]+1]=1;
    b[0]=power(a[0],mod-2);
    for(register int i=ite[0];i;--i)
    {
        init(ite[i]<<1);
        for(register int j=0;j<sup;++j)
        {
            tmp1[j]=j<ite[i]?a[j]:0;
            tmp2[j]=j<ite[i+1]?b[j]:0;
        }
        ntt(tmp1,1),ntt(tmp2,1);
        for(register int j=0;j<sup;++j)
            tmp1[j]=tmp2[j]*(2ll-tmp1[j]*tmp2[j]%mod+mod)%mod;
        ntt(tmp1,0);
        invsup=power(sup,mod-2);
        for(register int j=0;j<ite[i];++j)
            b[j]=tmp1[j]*invsup%mod;
    }
}

ll a[xx],b[xx],tmp3[xx],tmp4[xx],tmp5[xx],invb[xx];
int main()
{
    int n=read(),ite[21];
    for(register int i=0;i<21;++i) ite[i]=0;
    for(;n>1;n=(n+1)>>1) ite[++ite[0]]=n;
    ite[ite[0]+1]=1;
    for(register int i=0;i<ite[1];++i) a[i]=read();
    b[0]=1;
    for(register int i=ite[0];i;--i)
    {
        inv(b,invb,ite[i]);
        init(ite[i]<<1);
        for(register int j=0;j<sup;++j)
        {
            tmp3[j]=j<ite[i]?a[j]:0;
            tmp4[j]=j<ite[i+1]?b[j]:0;
            tmp5[j]=j<ite[i]?invb[j]:0;
        }
        ntt(tmp3,1),ntt(tmp4,1),ntt(tmp5,1);
        for(register int j=0;j<sup;++j)
            tmp3[j]=(((tmp3[j]+tmp4[j]*tmp4[j]%mod)%mod)*tmp5[j]%mod)*_2%mod;
        ntt(tmp3,0);
        invsup=power(sup,mod-2);
        for(register int j=0;j<ite[i];++j)
            b[j]=tmp3[j]*invsup%mod;
    }
    for(register int i=0;i<ite[1];++i)
        printf("%lld ",b[i]);
    puts("");
    return 0;
}

原文链接: https://www.cnblogs.com/cjsyh/p/17066347.html

欢迎关注

微信关注下方公众号，第一时间获取干货硬货；公众号内回复【pdf】免费获取数百本计算机经典书籍

    多项式求逆及其应用

原创文章受到原创版权保护。转载请注明出处：https://www.ccppcoding.com/archives/313863

非原创文章文中已经注明原地址，如有侵权，联系删除

关注公众号【高性能架构探索】，第一时间获取最新文章

转载文章受原作者版权保护。转载请注明原作者出处！

赞 (0)
adminadmin
0 0
上一篇 2023年2月16日 下午1:04
下一篇 2023年2月16日 下午1:04

相关推荐