多项式运算

I_am_WAQ · 2025-07-30 10:08:19 · 算法·理论

前置知识：数学分析系列。

以下文章中，我们记 [x^i]f(x) 表示取多项式 f (x) 中 i 次项的系数。举个例子，f (x) = 3x^2+2x+1，那么 [x^2]f(x) = 3。在计算机科学中，我们最为关注的正是多项式的系数，它们代表了各种丰富的意义。相比之下，对于不定元 x，我们关注的就比较少了，因此这又称为形式幂级数。下文中，若涉及无穷级数，因为我们不给 x 一个具体值，所以不必关注它的敛散性。

多项式运算，常被称为“多项式重工业”。也正是因此，有不少人认为多项式的每个知识都极其考验代码实现能力。但事实上，其整体代码都不算太难，重点知识只在于 \text{FFT}/\text{NTT} 和多项式牛顿迭代。如果理解得透彻，那么剩下的都是简单的套用。

由于本人代码习惯不好，层层封装，取模频繁，故常数较大，给出的运行时间几乎是最劣时间，鼓励各位读者尝试更优秀的写法。

为了保证其具有一定的速度，这里我们不再将各种操作都封装在结构体中。于是，这里首先给出结构体内的元素和一些基本函数，下文将不加解释的直接使用。

结构体中，n 代表多项式次数，N 代表大于 n 得最小 2 的整数次幂（方便后文操作），数组 p 用于存多项式系数，p_i = [x^i]f(x)，保证 p_n \ne 0。

clear 函数用于清空多项式，set 函数用于给多项式赋次数，cover 函数用于为多项式赋常数，tight 函数用于去除高位零。 :::info[参考实现]

struct poly
{
    int n,N,p[270000];
    void clear()
    {
        n=N=0;
        memset(p,0,sizeof p);
    }
    void set(int k)
    {
        n=k;
        N=1;
        while(N<=k)
            N<<=1;
    }
    void cover(int x)
    {
        set(0);
        memset(p,0,sizeof p);
        p[0]=x;
    }
    void tight()
    {
        while(n>=0&&!p[n])
            n--;
        if(n<0)
            n=0;
    }
};

:::

四则运算

首先介绍 poly 结构体一些最基本的运算（注意不是数学定义）。

程序实现中为了方便，我们将多项式对常数 k 的取模的运算重载为多项式对 x^k 取模，具体实现为截断。多项式与常数 k 的乘法，表示多项式数乘。times_x 函数代表多项式乘 x^k，用于优化多项式与单项式和简单多项式的乘法，具体实现为系数平移。在与确定的简单多项式相乘时，我们就可以调用该函数来 \text{O} (n) 完成。

:::info[参考实现]

poly times_x(poly f,int k)
{
    int i;
    poly P;
    P.clear();
    P.set(f.n+k);
    for(i=k;i<=P.n;i++)
        P.p[i]=f.p[i-k];
    return P;
}
poly operator %(poly F,int k)
{
    int i;
    poly P;
    P.clear();
    P.set(k-1);
    for(i=0;i<k;i++)
        P.p[i]=F.p[i];
    P.tight();
    return P;
}
poly operator *(int k,poly F)
{
    int i;
    poly P;
    P.clear();
    P.set(F.n);
    for(i=0;i<=F.n;i++)
        P.p[i]=(long long)k*F.p[i]%mod;
    P.tight();
    return P;
}

:::

加减

多项式加减的定义：[x^k](f(x) \pm g(x)) = [x^k]f(x) \pm [x^k]g(x)。

加减的实现是简单的，直接上代码。 :::info[参考实现]

poly operator +(poly F,poly G)
{
    int i;
    poly P;
    P.clear();
    P.set(max(F.n,G.n));
    for(i=0;i<=P.n;i++)
        P.p[i]=(F.p[i]+G.p[i])%mod;
    P.tight();
    return P;
}
poly operator -(poly F,poly G)
{
    int i;
    poly P;
    P.clear();
    P.set(max(F.n,G.n));
    for(i=0;i<=P.n;i++)
        P.p[i]=(F.p[i]-G.p[i]+mod)%mod;
    P.tight();
    return P;
}

:::

乘法

多项式乘法的定义：[x^k](f(x)g(x)) = \sum\limits_{i + j = k}[x^i]f(x)[x^j]g(x)（如果数学式子感觉复杂，就理解为直接按照分配律去乘，本人习惯写成卷积形式）。直接做是 \text{O} (n ^ 2) 的。

我们考虑，多项式不仅有系数表示法，还有点值表示法。一个 n 次多项式可以由 n + 1 个点唯一确定下来。对于多项式 \{(x_i,f(x_i))\} 和多项式 \{(x_i, g(x_i))\}，它们的乘法显然为 \{(x_i,f(x_i)g(x_i))\}。于是乘法可以 \text{O} (n) 完成。可是问题来了，系数表示更为常用，而二者间的转化均是 \text{O} (n ^ 2) 的，时间还是达不到要求。

这时，我们将引入快速傅里叶变换，它的作用就是优化系数表示和点值表示间的转化。

首先需要一些“复数”（[电子版]()）和“复变函数”（电子版）的知识。下文中将用 i 来代表循环变量，\mathrm{i} 来代表虚数单位，注意区分。

这里给出根本的解决方案：使用单位根。令 \omega_n = \mathrm{e}^{\frac{2\pi \mathrm{i}}{n}} = \cos \dfrac{2\pi}{n} + \mathrm{i} \sin \dfrac{2\pi}{n}，容易验证 x ^ n = 1 的根为 \{\omega_n^k|k \in \Z \cap [0, n)\}。

对于多项式 f (x) = \sum\limits_{k = 0} ^ {n - 1} a_kx^k，我们将 \omega_n^k 依次代入，令 b_k = f (\omega_n^k)，即得点值表示 \{b_k\}（n 个点确定了 n - 1 次多项式）。我们将 \{a_k\} \to \{b_k\} 称为快速傅里叶变换（\text{Fast Fourier Transform},\text{FFT}），\{b_k\} \to \{a_k\} 称为快速傅里叶逆变换（\text{Inverse Fast Fourier Transform},\text{IFFT}）。下面我们一个个讲解。

为了以下实现方便，默认 \log_2 n 为整数，若不是则可以向高位补 0，将次数补到最近的 N 使得 \log_2 N 为整数。以下统一记为 n。

首先是 \text{FFT}。

分类讨论。

于是，我们有：

于是 a_k = \dfrac{1}{n} \sum\limits_{i = 0} ^ {n - 1} b_i \omega_n^{-ki}。

证毕。 ::: 接下来考虑如何快速求解该式。我们发现，\text{FFT} 的本质是将 f_k 变成 f'_k = \sum\limits_{i = 0} ^ {n - 1} f_i\omega_n^{ki}。对 b 做一次 \text{FFT}，即得 b'_k = \sum\limits_{i = 0} ^ {n - 1}b_i\omega_n^{ki}，那么发现 a_k = \dfrac{1}{n}b'_{n - k}（注意到 a_0 = b_0）。于是我们只要对 b 做一次 \text{FFT}，就可以快速求出 a。

于是，\text{IFFT} 的时间复杂度与 \text{FFT} 一致，均为 \text{O} (n \log n)。

这样，我们就利用快速傅里叶变换在 \text{O} (n \log n) 的时间内完成了多项式的乘法。下面给出代码。 :::warning[注意事项]

1. 乘完后的多项式为 n + m 次，因此 N（上文分析中的 n）应当是最小大于 n + m 且 \log_2 N \in \Z 的数。于是，数组要开到 2(n+m)。
2. 由于有精度误差，最后的系数要四舍五入。 ::: :::info[参考实现] 封装了名为 comp 的复数结构体。

   void FFT(int n,comp* a)
   {
   int m=n>>1,i;
   if(n==1)
       return;
   comp rt=(comp){1,0},omega=(comp){cos(2*pi/n),sin(2*pi/n)},a0[m],a1[m];
   for(i=0;i<m;i++)
   {
       a0[i]=a[i<<1];
       a1[i]=a[i<<1|1];
   }
   FFT(m,a0);
   FFT(m,a1);
   for(i=0;i<m;i++)
   {
       a[i]=a0[i]+rt*a1[i];
       a[m+i]=a0[i]-rt*a1[i];
       rt=rt*omega;
   }
   }
   void IFFT(int n,comp* a)
   {
   int i;
   FFT(n,a);
   reverse(a+1,a+n);
   for(i=0;i<n;i++)
       a[i].a/=n;
   }

   poly operator *(poly a,poly b)
   {
   int i;
   comp shda[4100000],shdb[4100000];
   poly c;
   c.clear();
   c.set(a.n+b.n);
   for(i=0;i<c.N;i++)
   {
       shda[i]=(comp){(double)a.p[i],0};
       shdb[i]=(comp){(double)b.p[i],0};
   }
   FFT(c.N,shda);
   FFT(c.N,shdb);
   for(i=0;i<c.N;i++)
       shda[i]=shda[i]*shdb[i];
   IFFT(c.N,shda);
   c.n=N;
   for(i=0;i<c.N;i++)
       c.p[i]=shda[i].a+0.5;
   c.tight();
   return c;
   }

   :::

在 【模板】多项式乘法（FFT） 一题中，该代码时间 1.92 \text{s}，空间 195.23 \text{MB}。

注意到，\text{FFT} 中运用到了复数，精度较低，且难以取模。为了解决相关问题，这里我们引入快速数论变换（\text{Fast Number-Theoretic Transform},\text{NTT}）。这需要一点 原根 的知识。

由原根的性质，我们可以得到 g ^ {\varphi(p)} \equiv 1 \pmod p。那么，在模 p 意义下，\omega_n \equiv g ^ {\frac{\varphi(p)}{n}} \pmod p。于是，我们将 \text{FFT} 中的 \omega_n 全都换为 g ^ {\frac{\varphi(p)}{n}} 即可。有兴趣的读者可以证明在 \text{FFT} 中的结论换为原根后仍然成立。

接下来我们考虑选择合适的模数。一般地，\text{NTT} 的模数取质数 998244353，原因有二：

在对每个数进行二进制翻转，即得：

发现它就是有序的数组。事实上，这正是 \text{NTT} 数组变换的规律。

:::info[证明过程] 考虑 \text{NTT} 将奇偶分开的过程，将奇数（末位为 1）放在右边，偶数（末位为 0）放在左边，这正等价于将它们的二进制翻转排序（先排好最高位，再对最高位相同的排次高位，以此类推）。所以最终数组的二进制翻转正是有序的，即原数组本身。

证毕。 :::

将变换完后的数组记作 rev_i，现在我们考虑怎样求出来。一个简单的想法是，在求 i 时，\left\lfloor \dfrac{i}{2} \right\rfloor 的答案必然已经求出，那么我们可以利用它的答案。由于 \left\lfloor \dfrac{i}{2} \right\rfloor 二进制下等价于 i 去掉最低位，我们将 rev_{\lfloor \frac{i}{2} \rfloor} 去除最高位，然后根据 i 的奇偶性填上最高位即可。

这样求出 rev_i 后，我们就可以交换了。

向上回溯，更新数组的过程很简单。先枚举块长，再枚举块，最后对块内的数逐个更新即可。

这样，时间复杂度不变，为 \text{O} (n \log n)，而由于迭代版本不用对每一层开数组，而是全程只用到了 a 数组，空间复杂度降为 \text{O} (n)，同时常数大幅下降。

:::info[参考实现]

void NTT(int n,int* a)
{
    int len,m,h,u,v,rt,omega,i;
    for(i=0;i<n;i++)
    {
        if(rev[i]>i)
            swap(a[rev[i]],a[i]);
    }
    for(len=2;len<=n;len<<=1)
    {
        m=len>>1;
        omega=qpow(gnrt,(mod-1)/len);
        for(h=0;h<n;h+=len)
        {
            rt=1;
            for(i=h;i<h+m;i++)
            {
                u=a[i];
                v=(long long)rt*a[i+m]%mod;
                a[i]=(u+v)%mod;
                a[i+m]=(u-v+mod)%mod;
                rt=(long long)rt*omega%mod;
            }
        }
    }
}
void INTT(int n,int* a)
{
    int invn=qpow(n,mod-2),i;
    NTT(n,a);
    reverse(a+1,a+n);
    for(i=0;i<n;i++)
        a[i]=(long long)a[i]*invn%mod;
}

poly operator *(poly a,poly b)
{
    int i;
    poly c;
    c.clear();
    c.set(a.n+b.n);
    for(i=0;i<c.N;i++)
    {
        rev[i]=rev[i>>1]>>1;
        if(i&1)
            rev[i]|=c.N>>1;
    }
    NTT(c.N,a.p);
    NTT(c.N,b.p);
    for(i=0;i<c.N;i++)
        c.p[i]=(long long)a.p[i]*b.p[i]%mod;
    INTT(c.N,c.p);
    c.tight();
    return c;
}

:::

在 【模板】多项式乘法（FFT） 一题中，该代码时间 1.57 \text{s}，空间 44.71 \text{MB}。较递归形式，迭代形式时空都有了肉眼可见的改观。

除法

首先需要定义多项式的逆元。注意到多项式环中次数不为 0 的多项式不存在可逆元，因此一般意义下的多项式乘法逆不可定义（更准确地说，它们都是无穷级数）。我们只能在模意义下定义多项式乘法逆。

定义 f(x) 在模 x^n 意义下的乘法逆为 g(x) 满足 f(x)g(x) \equiv 1 \pmod {x^n}。现在我们考虑如何求解给定 n 下的 g(x)。为方便，我们同样假设 n 为 2 的整数次幂，若不是则高位补零。

注意这里有一个常见误区：将多项式改为点值表示，取逆后还原为系数表示。该做法错在多项式乘法逆为无穷级数，不能因为模 x^n 就只取 n + 1 个点来还原系数。

考虑求解乘法逆就是在模 x ^ n 意义下求解方程 P (g(x)) = \dfrac{1}{g(x)} - f(x) = 0（注意，我们在这里将多项式 g(x) 这一整体当作自变量，因此 f(x) 可看作常数）。

考虑求解方程的 牛顿迭代法，思考可不可以将它类比过来。由于牛顿迭代法具有平方收敛性，我们可以用模 x^{\frac{n}{2}} 的结果（不妨记为 g_0(x)）递推出 x^n 的结果。

将 P(g(x)) 在 g_0(x) 处做泰勒展开，得 P(g(x)) = \sum\limits_{i = 0} ^ \infty \dfrac{P^{(k)}(g_0(x))}{k!}(g(x)-g_0(x))^k。而后我们发现，在模 x^\frac{n}{2} 意义下，g_0(x) 和 g(x) 相等。那么，g (x) - g_0 (x) 的前 \dfrac{n}{2} 位应当相等。于是它们的 2 次方及更高次方都应当为 0。那么，我们就得到了 P(g(x)) = P(g_0(x)) + P'(g_0(x))(g(x)-g_0(x))。g(x) 作为对于模 x^n 意义下 P(g(x)) 的零点，我们可以解得 g (x) = g_0(x) - \dfrac{P(g_0(x))}{P'(g_0(x))}。不得不说，这和牛顿迭代的式子一模一样，而且满足完美的平方收敛（由 \dfrac{n}{2} 推出 n）。

对于 P (g(x)) = \dfrac{1}{g(x)} - f(x)，有 P' (g(x)) = -\dfrac{1}{g(x)^2}。那么 g(x) = g_0(x) - \dfrac{\frac{1}{g_0(x)} - f(x)}{-\frac{1}{g_0(x)^2}} = 2g_0(x) - f(x)g_0(x)^2。

依据此递推即可。运算次数满足 T (n) = T(\dfrac{n}{2}) + \text{O} (n \log n)，运用 主定理，即得时间复杂度 \text{O} (n\log n)。

下面给出实现。

:::warning[注意事项]

1. 多项式乘法逆计算中需计算 f(x)g_0(x)^2，若直接利用上文封装的乘法要做 6 次 \text{NTT}，而其中有很多次是不必要的。由于多项式乘法逆在后文中十分重要，我们考虑减少常数。对原式提取公因式，得 g_0(x)(2 - f(x)g_0(x))。若我们做完里面的 f(x)g_0(x) 后先取模，那么整体的次数即为 2n，我们可以取 2n 个点值。于是先对 f(x)，g_0(x) 做 \text{NTT}，然后按式子计算，最后对 g_0(x) 做 \text{INTT} 即可。这样就只用做 3 次 \text{NTT}。 ::: :::info[参考实现]

   poly poly_inv(poly P,int k)
   {
   int i,j;
   poly ans,f;
   P=P%k;
   ans.clear();
   ans.set(0);
   ans.p[0]=qpow(P.p[0],mod-2);
   for(i=2;i<=P.N;i<<=1)
   {
       f=P%i;
       for(j=0;j<(i<<1);j++)
       {
           rev[j]=rev[j>>1]>>1;
           if(j&1)
               rev[j]|=i;
       }
       NTT(i<<1,f.p);
       NTT(i<<1,ans.p);
       for(j=0;j<(i<<1);j++)
           ans.p[j]=(long long)ans.p[j]*(2ll-(long long)f.p[j]*ans.p[j]%mod+mod)%mod;
       INTT(i<<1,ans.p);
       ans=ans%i;
   }
   ans=ans%k;
   return ans;
   }

   :::

在 【模板】多项式乘法逆 一题中，该代码时间 723 \text{ms}，空间 8.14 \text{MB}。

然后就可以讨论多项式的除法了。注意，我们这里说的多项式除法，是多项式环中的带余除法，因此不可以直接乘上多项式乘法逆。

多项式除法的算法来源于系数翻转。我们记 f^R(x) 为 f(x) 系数翻转后的多项式。举个例子，若 f (x) = 4x^3 + x^2 + 2x + 3，那么 f^R(x) = 3x^3 + 2x^2 + x + 4。若 \deg f = n，则容易发现 f^R(x) = x^nf(\dfrac{1}{x})。

接下来开始推式子。下面记 \deg f = n， \deg g = m。f(x) = g(x) q(x) + r(x)。代入 \dfrac{1}{x}，得 f(\dfrac{1}{x}) = g(\dfrac{1}{x})q(\dfrac{1}{x})+r(\dfrac{1}{x})。两边同乘 x^n，得 x^nf(\dfrac{1}{x}) = x^mg(\dfrac{1}{x})\cdot x^{n - m} q(\dfrac{1}{x}) + x^{n - m + 1} \cdot x ^ {m - 1}r(\dfrac{1}{x})，即 f^R(x) = g^R(x)q^R(x)+x^{n - m + 1}r^R(x)。

由于 \deg q = n - m，我们不妨给两边同时对 x^{n - m + 1} 取模。那么，f^R(x) \equiv g^R(x)q^R(x) \pmod {x^{n - m + 1}}。于是，我们就得到了关键公式：

只需要求一次逆，做一次乘法，时间复杂度 \text{O} (n \log n)。得到 q(x) 后，r(x) 是容易的：

最后给出实现。 :::warning[注意事项]

1. 如果即要求商，也要求余数，建议求出 q 后直接手动用 f-gq 求 r，而不要调用多项式之间的模，因为该运算会再调用一次多项式整除，极大地增大了常数。 ::: :::info[参考实现]

   friend poly operator /(poly a,poly b)
   {
   poly q;
   reverse(a.p,a.p+a.n+1);
   reverse(b.p,b.p+b.n+1);
   q=a*poly_inv(b,a.n-b.n+1)%(a.n-b.n+1);
   reverse(q.p,q.p+q.n+1);
   return q;
   }
   friend poly operator %(poly a,poly b)
   {
   poly q,r;
   q=a/b;
   r=a-b*q;
   return r;
   }

   :::

在 【模板】多项式除法 一题中，该代码时间 1.14 \text{s}，空间 18.96 \text{MB}。若改为注意事项中提到的小常数写法，则时间 695 \text{ms}，空间 15.83 \text{MB}。

分析学操作

求导和积分的定义略过。

求导和积分的实现都是容易的，直接给出代码。 :::info[参考实现]

poly poly_diff(poly f)
{
    int i;
    poly P;
    P.clear();
    P.set(f.n-1);
    for(i=0;i<=P.n;i++)
        P.p[i]=(long long)(i+1)*f.p[i+1]%mod;
    return P;   
}
poly poly_intg(poly f)
{
    int i;
    poly P;
    P.clear();
    P.set(f.n+1);
    for(i=1;i<=P.n;i++)
        P.p[i]=(long long)f.p[i-1]*qpow(i,mod-2)%mod;
    return P;
}

:::

初等函数

这里的多项式初等函数与一般的复合函数不同。为了让多项式的初等函数均为 F[x] \mapsto F[x]，我们均采用“泰勒公式”（电子版）来定义。举个例子，对于 f (x) = x^2，一般的复合指数函数为 \mathrm{e} ^ {x^2}，而这里的多项式指数函数指的是 1 + x^2 + \dfrac{x^4}{2!} + \dfrac{x^6}{3!} + \cdots，后者一般取模 x ^ n 意义下的结果，此时二者的值不一定相等。

指对函数

先讲对数函数。注意到 \dfrac{\mathrm{d}}{\mathrm{d} x} \ln x = \dfrac{1}{x}，不再是超越函数，因此我们可以考虑利用求导和积分来转化。

考虑 \dfrac{\mathrm{d}}{\mathrm{d} x} \ln F(x) = \dfrac{F'(x)}{F(x)}，那么 \begin{aligned} \ln F(x) = \int \dfrac{F'(x)}{F(x)}\mathrm{d} x \end{aligned}。为了保证收敛，一般有 F(x) 常数项为 1，此时令 \ln F(x) 常数项为 0 即可。代入泰勒公式可知，若 F(x) 常数项为其他正整数，则 \ln F(x) 的常数项不收敛，此时多项式对数函数无意义（同时，这也是因为 \mathrm{e} 为超越数，在模意义下无法找到对应的数）。

根据该式子直接写即可。

:::info[参考实现]

poly poly_ln(poly f,int n)
{
    poly P;
    P.clear();
    P.set(n-1);
    P=poly_intg(poly_diff(f)*poly_inv(f,n)%n)%n;
    return P;
}

:::

在 【模板】多项式对数函数（多项式 ln） 一题中，以上代码时间 785 \text{ms}，空间 14.23 \text{MB}。

再讲指数函数。不妨令 n 为 2 的整数次幂。我们需要 g(x) \equiv \mathrm{e}^{f(x)} \pmod {x^n}，两边取对数，即得 \ln g(x) \equiv f(x) \pmod {x^n}，也即 \ln g(x) - f(x) \equiv 0 \pmod {x^n}。不妨令 P(g(x)) = \ln g(x) - f(x)，那么我们要做的就是在模 x^n 意义下求解方程 P(g(x)) = 0。

使用牛顿迭代法。对于 P(g(x)) = \ln g(x) - f(x)，有 P'(g(x)) = \dfrac{1}{g(x)}。于是我们就可以得到 g(x) = g_0(x) - \dfrac{\ln g_0(x) - f(x)}{\frac{1}{g_0(x)}} = g_0(x)(1 - \ln g_0(x) + f(x))。

分析时间复杂度，有 T(n) = T(\dfrac{n}{2}) + \text{O} (n \log n)，运用主定理即得时间复杂度 \text{O} (n \log n)。

最后给出代码。

:::info[参考实现]

poly poly_exp(poly f,int n)
{
    int i;
    poly P;
    P.clear();
    P.set(0);
    P.p[0]=1;
    f.p[0]=1;
    for(i=2;i<=n;i<<=1)
        P=P*(f%i-poly_ln(P,i))%i;
    return P;
}

:::

在 【模板】多项式指数函数（多项式 exp） 一题中，该代码时间 2.32 \text{s}，空间 19.48 \text{MB}。

多项式指数函数还有组合意义，详见 生成函数。

幂函数

这是多项式指对函数最直接的应用。

直接套快速幂即可做到 \text{O} (n \log n \log k) 复杂度。这显然不够，且不适用于非整数次幂。我们继续优化。

考虑 \ln 和 \exp 的性质。由于多项式指对函数由泰勒公式定义，它们与实数的指对函数一样，可以做指对变换。我们惊喜的发现，f (x)^k = \exp (k \ln f(x))。于是套用 \exp 和 \ln 即可做到 \text{O} (n \log n)。

给出代码实现。 :::warning[注意事项]

1. 该算法只适用于常数项为 1 的多项式，否则无法取 \ln。 ::: :::info[参考实现]

   poly poly_qpow(poly f,int k,int n)
   {
   poly P;
   P=poly_exp(k*poly_ln(f,n),n);
   return P;
   }

   ::: 在 【模板】多项式快速幂 一题中，该代码时间 7.92 \text{s}，空间 22.58 \text{MB}。

值得注意的一点是该题不寻常的指数，k 达到了 10 ^ {10 ^ 5} 级别，必须取模。一个误解是根据欧拉定理，k 应当对 \varphi (998244353) 取模。注意到是系数对 998244353 取模，而非多项式本身，所以这种运用是错误的。那该怎样取模呢？注意到公式 \exp (k \ln f(x)) 中，k 直接与 \ln f(x) 的系数相乘，而我们这里所谈的多项式实际上是域 F_{998244353} 上的形式幂级数，所以 k 可以直接对 998244353 取模。

由于后文中有运用，此处特殊的讲一下多项式开根。多项式开根可以看作 \dfrac{1}{2} 次方，然后套用上文的幂函数解法。这里讲解一种基于多项式牛顿迭代的快速做法。

多项式开根，即寻找 g (x) 满足 g(x)^2 \equiv f (x) \pmod {x^n}，也即求解方程 P(g(x)) = g(x)^2-f(x) = 0。记 g_0 (x) 为模 x ^ {\frac{n}{2}} 的答案，代入公式即得 g(x) = g_0(x) - \dfrac{g_0(x)^2-f(x)}{2g_0(x)} = \dfrac{1}{2} (g_0(x) + \dfrac{f(x)}{g_0(x)})。其中 \dfrac{1}{2} 在模 998244353 意义下取 499122177。

运算次数满足 T (n) = T(\dfrac{n}{2}) + \text{O} (n \log n)，运用主定理即得时间复杂度 \text{O} (n \log n)。

给出代码实现。

:::warning[注意事项]

1. 正如二次剩余有两个解，多项式开根也有两个解。代入二次剩余不同的解作为初值就可以获得两个解。下面的代码默认常数项为 1，取解中常数项较小的那个。 ::: :::info[参考实现]

   poly poly_sqrt(poly f,int n)
   {
   int i;
   poly P;
   P.clear();
   P.set(0);
   P.p[0]=1;
   for(i=2;i<=n;i<<=1)
       P=inv2*(f%i*poly_inv(P,i)%i+P);
   return P;
   }

   ::: 在 【模板】多项式开根 一题中，该代码时间 4.44 \text{s}，空间 14.32 \text{MB}。

三角函数

根据复变函数知识，我们有 \sin f(x) = \dfrac{\exp (\mathrm{i} f(x)) - \exp (-\mathrm{i}f(x))}{2 \mathrm{i}}，\cos f(x) = \dfrac{\exp (\mathrm{i}f(x)) + \exp (-\mathrm{i}f(x))}{2}。直接代入公式即可。其中，我们将 \mathrm{i} 看作 \omega_4，即可以避免复数参与计算。具体的，取 \mathrm{i} 为 911660635，-\mathrm{i} 和 \dfrac{1}{\mathrm{i}} 为 86583718 即可。时间复杂度显然 \text{O} (n \log n)。

:::info[参考实现]

poly poly_sin(poly f,int n)
{
    poly P;
    P=(long long)inv2*invimg%mod*(poly_exp(img*f,n)-poly_exp(invimg*f,n));
    return P;
}
poly poly_cos(poly f,int n)
{
    poly P;
    P=inv2*(poly_exp(img*f,n)+poly_exp(invimg*f,n));
    return P;
}

:::

在 多项式三角函数 一题中，该代码时间 5.28 \text{s}，空间 25.59 \text{MB}。

反三角函数

根据数学分析知识，容易得到 \begin{aligned} \arcsin f(x) = \int \dfrac{f'(x)}{\sqrt {1 - f(x)^2}}\mathrm{d} x \end{aligned}，\begin{aligned} \arctan f(x) = \int \dfrac{f'(x)}{1 + f(x)^2}\mathrm{d} x \end{aligned}。代入公式计算即可。时间复杂度显然 \text{O} (n \log n)。

:::warning[注意事项]

1. 多项式开根有两解。注意到 \begin{aligned} \arccos f(x) = -\int \dfrac{f'(x)}{\sqrt{1 - f(x)^2}}\mathrm{d} x \end{aligned}，所以两个解一个对应 \arcsin，另一个对应 \arccos。至于区分它们，至今没有什么好方法。 ::: :::info[参考实现]

   poly poly_arcsin(poly f,int n)
   {
   poly P,One;
   One.cover(1);
   P=poly_intg(poly_diff(f)*poly_inv(poly_sqrt(One-f*f%n,n),n)%n);
   return P;
   }
   poly poly_arctan(poly f,int n)
   {
   poly P,One;
   One.cover(1);
   P=poly_intg(poly_diff(f)*poly_inv(One+f*f%n,n)%n);
   return P;
   }

   :::

在 多项式反三角函数 一题（看不出来为什么它是黑，感觉和其它多项式题没什么本质区别）中，该代码时间 2.06 \text{s}，空间 25.14 \text{MB}。

以上各题均收录于 多项式模板题单。本文章中给出的代码保证可以通过所有题目，由此可见多项式各模板题没有什么卡常的风险（看都能看出本人的代码常数巨大，写法并不优），可以放心完成。