多项式（Ⅱ）：进阶工业

bloodstalk · 2023-06-26 11:16:22 · 个人记录

书接上回 多项式（Ⅰ）：基础工业。这部分主要写一下进阶的一些模板。

这部分你可能需要一点点导数和积分的知识。

多项式牛顿迭代

参考资料：

• 【数学】牛顿切线法在幂级数科技上的应用-nyy

• 多项式乘法逆题解-KAMIYA_KINA

牛顿迭代法，就是用来找函数零点的一个理论复杂度同二分但是实际上快而 nb 多的东西-nyy

其实这个东西在高数上册微分在求方程近似解上的应用里面出现过，其中的切线法，就是牛顿迭代。

下面我们讨论一下多项式的牛顿迭代。

给定多项式 G(x)，已知有 F(x) 满足：

G(F(x)) \equiv 0(\bmod \ x^n)

求出模 x^n 意义下的 F(x)

其实它跟我们下面要讲的多项式求逆的倍增法本质上是一样的，都是倍增的思想。

假设我们已知 F_0(x) \equiv F(x)(\bmod \ x^{\lceil \frac{n}{2} \rceil})，我们要求 F_1(x) \equiv F(x)(\bmod \ x^n)

下面我们先给出牛顿迭代的公式：

接下来我们简单证明一下：

我们在 F_0(x) 处用泰勒公式展开 G(F(x))

根据给出的条件，我们可以知道 \displaystyle F(x) - F_0(x) \equiv 0(\bmod \ x^{\lceil \frac{n}{2} \rceil}) 是成立的，所以说 \displaystyle (F(x)-F_0(x))^k \equiv 0(\bmod x^n) 在 k \ge 2 的时候成立。

所以说，我们泰勒展开后从二次项开始都是 0，有意义的只有前两项，因此：

又因为 G(F_1(x)) \equiv 0(\bmod \ x^n)，我们把方程整理一下，便有：

自此，公式得证。

可以看出，从 \lceil\frac{n}{2}\rceil 到 n，是一个 \log(n) 的过程，所以它的复杂度与二分和倍增是一样的，接下来的一些题目中，我们会运用到它来解题。

多项式求逆（Polyinv）

给定一个多项式 F(x)，求出一个多项式 G(x)，满足 F(x) * G(x) \equiv 1(\bmod \ x^n)。

多项式乘法逆

我们先讨论比较简单的模数是 998244353 这样类型的。

首先有个结论：若一个多项式有乘法逆，当且仅当它的常数项有逆元。

方法1：分治 FFT

目前还没学，咕掉先。时间复杂度是 O(n\log^2 n) 的。

方法2：倍增法

参考资料：多项式求逆-litble

如果多项式 F(x) 只有一项，那么显然 G_0 这个常数项就是 F_0 的逆元。

若有 n 项，考虑递归求解，像类似于归纳法一样往后推。

假设我们已经知道

因为 F(x) * G(x) \equiv 1(\bmod \ x^n)，这就说明对于第 1 \sim n-1 次方项，两个多项式乘起来后系数为 0，这就说明

两个式子相减

即

我们对上面这个式子进行平方。由于 G(x) - H(x) 在模 x^{\lceil \frac{n}{2}\rceil} 意义下为 0，说明从 0\sim \lceil \frac{n}{2}\rceil-1 次项都为 0。设平方后的这个多项式为 P，则

那么对于 i \leq n 的项，j 和 i-j 至少有一项的次数 < \lceil \frac{n}{2 }\rceil，那么就是 0，所以

乘上一个 F(x)，根据 F(x)G(x) \equiv 1(\bmod x^n)，再移一下项，就得到

转化到了这个地方便可以用 NTT 求解了。

复杂度用主定理分析一下：

此外还有一些常用的递归复杂度：

可以用主定理分析一下，感觉基本的这么多就可以了。

代码实现

有一些小细节需要注意：

• 代码中 (n+1)>>1 就是向上取整，至于为什么是向上取整，这样就能保证 \lceil \frac{n}{2} \rceil 的平方一定是能包含到 n 的。

• 为什么 lim < (n<<1) ， 因为这样相乘以后就能保证两个 n 项的多项式相乘后的多项式能存的开。

• 为什么要把 n \sim \text{lim}-1 设成 0，因为模 x^n 意义下这些地方是同余于 0 的，为了避免干扰计算，就把系数也设成 0 了。后面的那个 b_i 设 0 也是同理的。

il void Polyinv(int n,ll a[],ll b[])
    {
        if(n == 1) { b[0] = inv(a[0]); return ; }
        Polyinv((n+1)>>1,a,b);
        static ll c[N];
        lim = 1;
        while(lim < (n<<1)) lim <<= 1;
        for(re int i=0;i<lim;i++) R[i] = (R[i>>1]>>1) + ((i&1) ? lim/2 : 0);
        for(re int i=0;i<n;i++) c[i] = a[i];
        for(re int i=n;i<lim;i++) c[i] = 0;
        NTT(c,lim,1) , NTT(b,lim,1);
        for(re int i=0;i<lim;i++) b[i] = 1ll * ((2 * b[i] % mod) - (c[i] * b[i] % mod * b[i] % mod) + mod) % mod;
        NTT(b,lim,-1);
        for(re int i=n;i<lim;i++) b[i] = 0;
        return ;
    }

方法3：牛顿迭代法

其实牛顿迭代法推出来的式子是和倍增法一样的，不过牛顿迭代法用途更为广泛。

我们先把 F(x) * G(x) \equiv 1(\bmod \ x^n) 来转化一下，把它们拆成相加的形式，那么就有

你会发现这符合牛顿迭代法的样子，我们像牛顿迭代法那样的构造一个外层函数 H(t)，使它 能把上面这个式子表达出来，因为我们是要求 G(x) 的，所以 G(x) 应该放在 H(t) 里，那么就构造

令 t = G(x)，我们便得到了我们想要的式子：H(G(x)) \equiv 0(\bmod \ x^n)

类似的，假设我们已知 G_0(x) \equiv G(x)(\bmod \ x^{\lceil \frac{n}{2} \rceil})，我们要求 G_1(x) \equiv G(x)(\bmod \ x^n)

我们将它套入牛顿迭代公式，则有

把 H 的解析式套进去，然后求导的时候要注意，我们是对 G(x) 求导，所以这里的 F(x) 是一个常数，常数求导后等于 0，那么有

化简一下，则有

你会发现这和我们倍增得到的结果是一样的，这也说明它们本质上其实是相同的。

任意模数多项式乘法逆

把 NTT 换成 MTT ，然后再注意亿点细节即可。

感谢万能的 UOJ 群友解答。

三模 NTT 版本

三模 NTT 版本常数巨大无比，由于相乘会爆 long long，所以在一次递归中，三模 NTT 需要调用 18 次 DFT。常数可想而知，又加上题解区三模 NTT 版本的少之又少，导致我调了半上午 + 一下午才有了一个具体的理解。下面进入正文。

首先，因为 CRT 的合并会导致数很大，所以我们分成两次计算。

首先，我们计算出三个模数下的 F(x) \times H(x)，然后我们直接把它用 exCRT 合并起来，目的是直接把这个转化成 \bmod \ 10^9+7 的意义下(模板题给定的模数是 10^9+7)进行，这样就不会出锅了。

之后，我们把它取反，变成 -F(x)H(x)，然后再将其变成 2 - F(x)H(x)，也就是在这个多项式的 x^0 项加上一个 2。这么简单的问题我竟然当时没绕过弯来，其实就是这个 2 其实是 x^0 次方的系数，你不能在所有的 x^k 的系数前都加上 2，这样显然不对。

然后我们在让这个多项式与 H(x) 相乘，这样便得到了 2H(x) - F(x)H(x)^2，也就是我们的 G(x)。

还是，注意取模。

代码

只截取了主要部分，其中 inv 函数求逆元，是用快速幂实现的。

#include<bits/stdc++.h>
#define ll long long
#define next nxt
#define re register
#define il inline
const int N = 4e5 + 5;
const ll mod1 = 469762049;
const ll mod2 = 998244353;
const ll mod3 = 1004535809;
const ll Mod = 1e9 + 7;
const ll g = 3;
using namespace std;

int n,m,lim;
ll invg,invn,inv1,inv12;
ll a[N],b1[N],b2[N],b3[N],R[N],B[N],aa[N],bb[N];

namespace Poly
{
    il void NTT(ll A[],int n,int type,ll mod)
    {
        //do something
    }
    il ll crt(ll a,ll b,ll c)//采用 exCRT 的合并方法就不会爆 long long了
    {
        ll k = ((b-a) % mod2 + mod2) % mod2 * inv1 % mod2;
        ll x = k * mod1 + a;
        k = ((c-x) % mod3 + mod3) % mod3 * inv12 % mod3;
        x = (x + k * mod1 % Mod * mod2 % Mod) % Mod;
        return x;
    }
    il void Mul(ll *A,ll *B,ll *C,int n,ll modx)
    {
        for(re int i=0;i<n;i++) aa[i] = A[i] , bb[i] = B[i];
        for(re int i=n;i<lim;i++) aa[i] = bb[i] = 0;
        NTT(aa,lim,1,modx) , NTT(bb,lim,1,modx);
        for(re int i=0;i<lim;i++) C[i] = aa[i] * bb[i] % modx;
        NTT(C,lim,-1,modx);
        for(re int i=n;i<lim;i++) C[i] = 0;
        return ;
    }
    il void Polyinv(int n)
    {
        if(n == 1) { B[0] = inv(a[0],Mod); return ; }
        Polyinv((n+1)>>1);
        lim = 1;
        while(lim < (n<<1)) lim <<= 1;
        for(re int i=0;i<lim;i++) R[i] = (R[i>>1]>>1) + ((i&1) ? lim/2 : 0);
        Mul(a,B,b1,n,mod1) , Mul(a,B,b2,n,mod2) , Mul(a,B,b3,n,mod3);//求 F(x)H(x)   
        for(re int i=0;i<n;i++) b1[i] = b2[i] = b3[i] = (Mod - crt(b1[i],b2[i],b3[i])) % Mod;//取反
        b1[0] += 2 , b2[0] += 2 , b3[0] += 2;//x^0项加2
        Mul(b1,B,b1,n,mod1) , Mul(b2,B,b2,n,mod2) , Mul(b3,B,b3,n,mod3);//求 2H(x) - F(x)H(x)^2
        for(re int i=0;i<n;i++) B[i] = crt(b1[i],b2[i],b3[i]);//合并之后就是答案
        return ;
    }
}
using namespace Poly;

signed main()
{
    n = read();
    for(re int i=0;i<n;i++) a[i] = read() % Mod;
    inv1 = inv(mod1 % mod2,mod2) , inv12 = inv(mod1 * mod2 % mod3,mod3);
    Poly::Polyinv(n);
    for(re int i=0;i<n;i++) cout << B[i] << " ";
    return 0;
}

拆系数 FFT 版本

发现拆系数 FFT 还是比三模 NTT 好写+好理解，相比于正常的拆系数 FFT 没有什么大变动。

仍然只保留主要部分。我采用的是 4 次 FFT 版本，相比于三模 NTT 来说，它每轮只需要 8 次调用 DFT，所以最后的速度比 NTT 快了一倍左右，FFT 是 10s，NTT 接近 20s。

#include<bits/stdc++.h>
#define int long long
#define ll long long
#define complex comple
#define double long double
#define re register
#define il inline
const int N = 4e5 + 5;
const int mod = 1e9 + 7;
const double Pi = acos(-1.0);
const ll bas = (1<<15) , bas2 = bas * bas;
using namespace std;
int max(int x,int y){return x > y ? x : y;}
int min(int x,int y){return x < y ? x : y;}

int n,m,lim;
ll a[N],b[N],c[N],R[N];
ll x,AC,AD,BC,BD;

namespace Poly
{
    struct complex
    {
        double x,y;
        complex (double xx = 0 , double yy = 0) { x = xx; y = yy; }
    }f[N],g[N],h[N];
    //运算符重载

    il ll ksm(ll a,ll b)
    {
        //do something
    } ll inv(ll x) { return ksm(x,mod-2); }
    il void FFT(complex A[],int n,int type)
    {
        // do something
    }
    il void Mul(ll a[],ll b[],ll c[],int n)
    {
        for(re int i=0;i<lim;i++) f[i] = g[i] = h[i] = {0,0};
        for(re int i=0;i<n;i++)
        {
            a[i] %= mod , b[i] %= mod;
            f[i].x = (a[i]>>15) , f[i].y = a[i] & 32767;
            h[i].x = (b[i]>>15) , h[i].y = b[i] & 32767;
        }//优化 4 次 FFT
        FFT(f,lim,1) , FFT(h,lim,1);
        g[0] = {f[0].x,-f[0].y};
        for(re int i=1;i<lim;i++) g[i] = {f[lim-i].x,-f[lim-i].y};
        for(re int i=0;i<lim;i++) h[i].x /= lim , h[i].y /= lim , f[i] = f[i] * h[i] , g[i] = g[i] * h[i];
        FFT(f,lim,-1) , FFT(g,lim,-1);
        for(re int i=0;i<n;i++)
        {
            AC = ((ll)((f[i].x+g[i].x)/2+0.5)) % mod;
            AD = ((ll)((f[i].y+g[i].y)/2+0.5)) % mod;
            BD = ((ll)(g[i].x-AC+0.5) % mod + mod) % mod;
            BC = ((ll)(f[i].y-AD+0.5) % mod + mod) % mod;
            AC = bas2 % mod * AC % mod , AD = bas % mod * (AD+BC) % mod;
            c[i] = (AC + AD + BD) % mod;
        }
    }
    il void Polyinv(int n)
    {
        if(n == 1) { b[0] = inv(a[0]); return ; }
        Polyinv((n+1)>>1);
        lim = 1;
        while(lim < (n<<1)) lim <<= 1;
        for(re int i=0;i<lim;i++) R[i] = (R[i>>1]>>1) + ((i&1) ? lim/2 : 0) , c[i] = 0;
        Mul(a,b,c,n);//算第一次
        for(re int i=0;i<n;i++) c[i] = (mod - c[i]) % mod;
        c[0] = (c[0] + 2) % mod;//取反求第二次
        Mul(b,c,b,n);
        for(re int i=n;i<lim;i++) b[i] = 0;
        return ;
    }
}
using namespace Poly;

signed main()
{
    n = read();
    for(re int i=0;i<n;i++) a[i] = read();
    Poly::Polyinv(n);
    for(re int i=0;i<n;i++) cout << b[i] << " ";
    return 0;
}

多项式对数函数（Polyln）

给出 n-1 次多项式 A(x)，求一个 \bmod{\:x^n} 下的多项式 B(x)，满足 B(x) \equiv \ln A(x).

在 \text{mod } 998244353 下进行，且 a_i \in [0, 998244353) \cap \mathbb{Z}

首先若 \ln f(x) 存在，则必有 [x^0]f(x) = 1，因为 e 在模意义下没有对应的值，只有 \ln 1 = 0 是确定的。

在这种情况下，我们考虑给 f 数组进行变换。

因为一个多项式先求导再积分就是原函数加个常数，我们令这个常数为 0，然后就有了下面的变换：

这是一个复合函数的形式，我们先给它导一下

到了这步其实就很简单了。

运用导数的知识，我们可以把分子求出来，多项式的导数还是很简单的，就是每项分别求导再相加就行了，运用 (x^a)' = ax^{a-1}，然后相加就行了，然后注意第 n-1 项是 0 就行。

然后下面的这个式子我们可以用多项式求逆整出来，这个我们也会。

然后就是把这个求出来的多项式积分一下， 这个也好整，运用 \displaystyle \int kx^a = \frac{k}{a+1}x^{a+1} + C（C 是常数） 这个式子往里代就行，同样的也是对每项积分再相加，这样就得到了最终的式子。

时间复杂度是多项式求逆和 NTT 的 O(n\log n)

代码

#include<bits/stdc++.h>
//#define int long long
#define ll long long
#define next nxt
#define re register
#define il inline
const int N = 4e5 + 5;
const int mod = 998244353;
const ll g = 3;
using namespace std;

int n,m,lim,invg,invn;
int R[N];
ll a[N],b[N];

namespace Poly
{
    il ll ksm(ll a,ll b)
    {
        //快速幂
    } ll inv(ll x) { return ksm(x,mod-2); }
    il void NTT(ll A[],int n,int type)
    {
        // do something
    }
    il void Polymul(int n,ll a[],ll b[])
    {
        //多项式相乘
    }
    il void Polyinv(int n,ll a[],ll b[])
    {
        //多项式求逆
    }
    il void Polyderiv(int n,ll a[],ll b[])
    {
        for(re int i=0;i<n-1;i++)
            b[i] = a[i+1] * (i+1) % mod;
        b[n-1] = 0;
    }
    il void Polyinte(int n,ll a[],ll b[])
    {
        b[0] = 0;
        for(re int i=1;i<=n;i++) b[i] = a[i-1] * inv(i) % mod;
    }
    il void Polyln(int n,ll a[],ll b[])
    {
        static ll F[N],G[N];
        memset(F , 0 , sizeof F);
        memset(G , 0 , sizeof G);
        Polyderiv(n,a,F);
        Polyinv(n,a,G);
        Polymul(n,F,G);
        Polyinte(n,F,b);
    }
}
using namespace Poly;

signed main()
{
    n = read();
    for(re int i=0;i<n;i++) a[i] = read();
    Poly::Polyln(n,a,b);
    for(re int i=0;i<n;i++) cout << b[i] << " ";
    return 0;
}

多项式指数函数（Polyexp）

给出 n-1 次多项式 F(x)，求一个 \bmod{\:x^n} 下的多项式 G(x)，满足 G(x) \equiv \text e^{F(x)}。系数对 998244353 取模。

如果你细心的话，你会发现洛谷上对数函数和指数函数的样例刚好是反过来的，这也恰恰说明了这俩是逆运算（其实这是废话）。

和对数函数相同的，我们要保证 [x^0]A(x)= 0，这样才能使得 B(x) 的常数项在模意义下有意义。

在这里我们会用到牛顿迭代法来解决。

首先，我们可以给左右两边取个 \ln，那么就有：

移项，

我们发现这又变成了牛顿迭代的形式了，我们再令 G(x) 为自变量，给它套一个外层函数 H(t)，那么就是

将 G(x) 代入

还是一样的，假设已知 G_0(x) \equiv F(x)(\bmod \ x^{\lceil \frac{n}{2} \rceil})，我们要求 G_1(x) \equiv F(x)(\bmod \ x^n)

代入牛顿迭代公式，就有

把 H 的解析式套进去

化简一下

其中，G_0(x) 我们是已知的，这也就是说，我们可以采取类似于多项式求逆的递归，从小到大的把 G(x) 求出来。

时间复杂度 T(n) = T(n/2) + O(n\log n) = O(n\log n)

代码

#include<bits/stdc++.h>
//#define int long long
#define ll long long
#define next nxt
#define re register
#define il inline
const int N = 4e5 + 5;
const int mod = 998244353;
const ll g = 3;
using namespace std;

int n,m,lim,invg,invn;
int R[N];
ll a[N],b[N];

il int read()
{
    int f=0,s=0;
    char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) f |= (ch=='-');
    for(; isdigit(ch);ch=getchar()) s = (s<<1) + (s<<3) + (ch^48);
    return f ? -s : s;
}

namespace Poly
{
    il ll ksm(ll a,ll b)
    {
        //do something
    } ll inv(ll x) { return ksm(x,mod-2); }
    il void NTT(ll A[],int n,int type)
    il void Polymul(int n,ll a[],ll b[])
    il void Polyinv(int n,ll a[],ll b[])
    il void Polyderiv(int n,ll a[],ll b[])
    il void Polyinte(int n,ll a[],ll b[])

   il void Polyln(int n,ll a[],ll b[])
    {
        static ll F[N],G[N];
        memset(F , 0 , sizeof F);
        memset(G , 0 , sizeof G);
        Polyderiv(n,a,F);
        Polyinv(n,a,G);
        Polymul(n,F,G);
        Polyinte(n,F,b);
    }
    il void Polyexp(int n,ll a[],ll b[])
    {
        if(n == 1) { b[0] = 1; return ; }
        Polyexp((n+1)>>1,a,b);
        static ll lnb[N];
        Polyln(n,b,lnb);
        lim = 1;
        while(lim < (n<<1)) lim <<= 1;
        for(re int i=0;i<n;i++) lnb[i] = ((a[i] - lnb[i]) % mod + mod) % mod;
        for(re int i=n;i<lim;i++) lnb[i] = b[i] = 0;
        lnb[0]++;
        Polymul(n,b,lnb);
        for(re int i=n;i<lim;i++) b[i] = 0;
    }
}
using namespace Poly;

signed main()
{
    n = read();
    for(re int i=0;i<n;i++) a[i] = read();
   invg = inv(g);
    Poly::Polyexp(n,a,b);
    for(re int i=0;i<n;i++) cout << b[i] << " ";
    return 0;
}

多项式幂函数（Polypow）

因为第一遍跑代码的时候 invg 不会被预处理导致调了 3h+。

普通版

给定一个 n-1 次多项式 F(x)，求一个在 \bmod\ x^n 意义下的多项式 G(x)，使得 G(x) \equiv (F(x))^k \ (\bmod\ x^n)。

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

保证 [x^0]F(x)=0。

首先朴素暴力快速幂肯定是 O(n\log n \log k) 的，但实际上我们可以继续优化它，这需要用到对数函数的性质：\ln x^a = a\ln x。

我们先把 F(x)^k 取个 \ln，变成 \ln(F(x)^k)，然后再应用上面的这个性质，就变成了 k\ln(F(x))，然后再对它求个 exp，于是这题就做完了，没错，就这么简单，并且它保证了 [x^0]F(x) = 0，意味着可以取 \ln。

然后就有我上面说的，因为第一遍跑代码的时候 invg 不会被预处理导致调了 3h+，警钟敲烂。

还有一个问题是 k 的大小，在原题中 k \leq 10^{10^5}，很大。然后 EI 太阳神好像证明了一下当 n<p 时，一个多项式的幂次 k 直接对 p 取模是对的，具体可以看题解里的证明。

代码

#include<bits/stdc++.h>
#define int long long
#define ll long long
#define next nxt
#define re register
#define il inline
const int N = 4e5 + 5;
const int mod = 998244353;
const ll g = 3;
using namespace std;
int max(int x,int y){return x > y ? x : y;}
int min(int x,int y){return x < y ? x : y;}

int n,m,lim,invg,invn;
int R[N];
ll k,a[N],b[N];

namespace Poly
{
    il ll ksm(ll a,ll b)
    il ll inv(ll x) { return ksm(x,mod-2); }
    il void NTT(ll A[],int n,int type)
    il void Polymul(int n,ll a[],ll b[])
    il void Polyinv(int n,ll a[],ll b[])
    il void Polyderiv(int n,ll a[],ll b[])
    il void Polyinte(int n,ll a[],ll b[])
    il void Polyln(int n,ll a[],ll b[])
    il void Polyexp(int n,ll a[],ll b[])

    il void Polypow(int n,ll k,ll a[],ll b[])
    {
        static ll lna[N];
        Polyln(n,a,lna);
        for(re int i=0;i<n;i++) lna[i] = lna[i] * k % mod;
        Polyexp(n,lna,b);
        memset(lna , 0 , sizeof lna);
    }
}
using namespace Poly;

signed main()
{
    n = read() , k = read();
    invg = inv(g);//because of you
    for(re int i=0;i<n;i++) a[i] = read();
    Poly::Polypow(n,k,a,b);
    for(re int i=0;i<n;i++) cout << b[i] << " ";
    return 0;
}

加强版

加强版和普通版就差了一个地方，并不保证 a_0 = 1 。

为此，我们应该把这个多项式转化成 a_0=1 的样子。

具体的，我们从 F(x) 中提取一个公因式 b \times x^k，其中 k 满足 a_k \not = 0，且 \forall 0 \leq i <k,a_i =0，b 就是 x^k 这一项的系数，这样就可以保证最低位是 1 了。

举个例子，我们现在是 0 + 2x + 2x^2，我们需要提取出一个 2x，式子就变成了 2x(1+x)，然后对右边的多项式进行快速幂，然后再乘上左边的系数的幂次，再右移即可。

注意一种特殊情况：k \times m \ge n（这里的 k 是上文的 k，m 就是题目中给定的 k，也就是幂次），这也就是说前 k\times m 项都是 0 了，我们也就没必要往后做了，输出 n 个 0 即可。

还有这个题的取模问题，我们提取出来的公因式的次幂，其实是对 \varphi(p) 取模的，因为它是一个单项式，应该用扩展欧拉定理来降幂，也就是模 \varphi(p)，右边的那个多项式还是同理的对 p 取模。

代码

根据上文，你需要存三个 k：