多项式整理

时间 2019-11-06

标签多项式整理繁體版

原文原文链接

多项式求逆元

多项式求逆元，即已知多项式$A(x)$，咱们须要找到一个多项式$A^{-1}(x)$算法

使得函数

$$A(x)A^{-1}(x)\equiv 1\pmod {x^n}$$ui

咱们称多项式$A^{-1}(x)$为多项式$A(x)$的逆元spa

在这里${x^n}$是一个数，模${x^n}$的意义就是将$\ge n$的项都忽略掉.net

至于咱们为何要模$x^n$，由于经过计算不难发现：除了仅有常数项的多项式的逆元为一个常数以外，其他多项式的逆元均有无穷多项code

算法

这里介绍一种比较经常使用的$O(nlogn)$倍增算法，实际上许多与多项式有关的操做都须要用的倍增算法blog

假设咱们已经求出了多项式$A(x)$在模$x^{\frac{n}{2}}$意义下(就是一半)的逆元，设其为$B(x)$ip

即get

$$A(x)B(x)\equiv 1\pmod {x^{\lceil\frac n2\rceil}}$$string

根据取模运算的性质，$A^{-1}(x)$前$\frac{n}{2}$项与$B(x)$是相同的

$$A(x)A^{-1}(x)\equiv 1\pmod {x^{\lceil\frac n2\rceil}}$$

合并两个式子能够获得

$$A(x)(B(x)-A^{-1}(x))\equiv 0\pmod{x^{\lceil\frac n2\rceil}}$$

这里的$A(x)$不为零

$$B(x)-A^{-1}(x)\equiv 0\pmod{x^{\lceil\frac n2\rceil}}$$

而后两边平方后拆开

$$B^2(x)+A^{-2}(x)-2B(x)A^{-1}(x)\equiv 0\pmod {x^n}$$

两边同乘$A(x)$

$$B^2(x)A(x)+A^{-1}(x)-2B(x)\equiv 0\pmod {x^n}$$

移项

$$A^{-1}(x)\equiv B(x)(2-B(x)*A(x))\pmod {x^n}$$

这样咱们就获得了$A(x)$和$B(x)$的关系

利用NTT计算复杂度为$O(nlogn)$

代码实现

洛谷P4238 【模板】多项式求逆

一份常数还不错的代码

// luogu-judger-enable-o2
#include<cstdio>
#include<cstring>
#include<algorithm>
#define getchar() (p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 1<<21, stdin), p1 == p2) ? EOF : *p1++)
#define swap(x,y) (x ^= y, y ^= x, x ^= y)
#define mul(a, b) 1ll * a * b % P
#define add(a, b) (a + b >= P ? a + b - P : a + b)
#define dec(a, b) (a - b <  0 ? a - b + P : a - b)
#define rg register 
const int MAXN = (1 << 21) + 10, P = 998244353, G = 3, Gi = 332748118; 
char buf[1<<21], *p1 = buf, *p2 = buf;
inline int read() { 
    char c = getchar(); int x = 0, f = 1;
    while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();}
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x * f;
}
char obuf[1<<24], *O=obuf;
void print(int x) {
    if(x > 9) print(x / 10);
    *O++= x % 10 + '0';
}
inline int fastpow(int a, int k) {
    int  base = 1;
    while(k) {
        if(k & 1) base = mul(a, base);
        a = mul(a, a); k >>= 1;
    }
    return base % P;
}
int N, r[MAXN], X[MAXN], Y[MAXN], A[MAXN], B[MAXN], Og[MAXN];
inline void NTT(int  *A, int type, int len) {
    int limit = 1, L = 0;
    while(limit < len) limit <<= 1, L++;
    for(rg int i = 0; i < limit; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));    
    for(rg int i = 0; i < limit; i++) if(i < r[i]) swap(A[i], A[r[i]]);
    for(rg int mid = 1; mid < limit; mid <<= 1) {    
        int  R = mid << 1;
        int  W = fastpow(G, (P - 1) / R); Og[0] = 1;
        for(rg int j = 1; j < mid; j++) Og[j] = mul(Og[j - 1], W);
        for(rg int j = 0; j < limit; j += R) {
            for(rg int k = 0; k < mid; k++) {
                 const int x = A[j + k], y = mul(Og[k], A[j + k + mid]);
                 A[j + k] = add(x, y), A[j + k + mid] = dec(x, y);
            }
        }
    }
    if(type == -1) {
        std::reverse(&A[1], &A[limit]);
        for(int i = 0, inv = fastpow(len , P - 2); i < limit; i++)
            A[i] = 1ll * A[i] * inv % P;
    }
}
void Inv(int *a, int *b, int len) {// a要求的多项式 b逆元 len：要求的逆元的长度 
    if(len == 1) {
        b[0] = fastpow(a[0], P - 2);
        return ;
    }
    Inv(a, b, len >> 1);
    for(rg int i = 0; i < len; i++) A[i] = a[i], B[i] = b[i];
    NTT(A, 1, len << 1); NTT(B, 1, len << 1);
    for(rg int i = 0; i < (len << 1); i++) A[i] = mul(mul(A[i], B[i]), B[i]) ;
    NTT(A, -1, len << 1);
    for(rg int i = 0; i < len; i++) b[i] = (1ll * (b[i] << 1) % P + P - A[i] ) % P;
}
int main() {
    #ifdef WIN32
    freopen("a.in","r",stdin);
    #endif 
    N = read();
    for(int i = 0; i < N; i++) X[i] = (read() + P) % P;
       int Len; for(Len = 1; Len < N; Len <<= 1);
       Inv(X, Y, Len);
       for(int i = 0; i < N; i++) print(Y[i]), *O++ = ' ';
    fwrite(obuf, O-obuf, 1 , stdout);
    return 0;
}

多项式除法

给定多项式$A(x)$,$B(x)$

咱们须要找到多项式$D(x)$，$R(x)$，使得

$$A(x) = D(x)B(x) + R(x)$$

在这里$A(x)$为$N$次多项式，$B(x)$为$M$次多项式

那么$D(x)$为$N-M+1$次多项式，$R(x)$的次数$<M$，这样咱们能够保证解的惟一性

算法

咱们考虑如何解决上面的问题

首先$R(x)$具体的值是不用考虑的，由于咱们求出$D(x)$后能够把$D(x)$带入从而求得$R(x)$

另外，根据咱们求逆元的经验，没有模的多项式除法是有无穷多项的

这个其实也好解决，咱们强制让全部多项式$\pmod {x^{n - m +1}}$

那么接下来咱们只须要解决余数，也就是$R(x)$的问题了

考虑到咱们的模数为$x^{n - m +1}$，也就是说咱们会丢弃全部多项式的前$n-m+1$项

可是$R(x)$是一个$M-1$次多项式，直接模确定是消不掉的，咱们考虑能不能让它的系数乘上$x^{n-m+1}$还能保证要求的多项式跟原来多项式意义相同

这里，咱们定义翻转操做

$$A^R(x) = x^n A(\frac{1}{x}) $$

也就是将多项式的系数进行翻转

下面是神仙推导

$$A(x)=B(x)D(x)+R(x)$$

将$x$用$\frac{1}{x}$替代，两边同乘$x^N$

$$\begin{eqnarray*} x^n A(\frac{1}{x}) &=& x^{n - m}D(\frac{1}{x}) x^mB(\frac{1}{x}) + x^{n - m + 1}x^{m - 1}R(\frac{1}{x}) \\ A^R(x) &=& D^R(x)B^R(x) + x^{n - m + 1}R^R(x) \end{eqnarray*}$$

而后经过$\pmod {x ^{n - m +1}}$就能够把$R(x)$消掉啦

获得

$$A^R(x)\equiv B^R(x)D^R(x)\pmod{x^{n-m+1}}$$

$$A^R(x) * B^{-R}(x)\equiv D^R(x)\pmod{x^{n-m+1}}$$

这样的话咱们求出$B(x)$的逆元，而后用NTT乘起来就行了

多项式取模

实际就是上面的$R(x)$

用多项式除法作就能够

代码实现

// luogu-judger-enable-o2
// luogu-judger-enable-o2
#include<cstdio>
#include<cstring>
#include<algorithm>
#define getchar() (p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 1<<21, stdin), p1 == p2) ? EOF : *p1++)
#define swap(x,y) (x ^= y, y ^= x, x ^= y)
#define mul(a, b) 1ll * a * b % P
#define add(a, b) (a + b >= P ? a + b - P : a + b)
#define dec(a, b) (a - b <  0 ? a - b + P : a - b)
#define rg register 
using namespace std;
const int MAXN = 1e6, P = 998244353, Gi = 3; 
char buf[1<<21], *p1 = buf, *p2 = buf;
inline int read() { 
    char c = getchar(); int x = 0, f = 1;
    while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();}
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x * f;
}
char obuf[1<<24], *O=obuf;
void print(int x) {
    if(x > 9) print(x / 10);
    *O++= x % 10 + '0';
}
inline int fastpow(int a, int k) {
    int  base = 1;
    while(k) {
        if(k & 1) base = mul(a, base);
        a = mul(a, a); k >>= 1;
    }
    return base % P;
}
int N, r[MAXN], A[MAXN], B[MAXN], Og[MAXN], F[MAXN], Q[MAXN], G[MAXN], R[MAXN], Ginv[MAXN], Atmp[MAXN], Btmp[MAXN];
int GetLen(int x) {
    int len;
    for(len = 1; len <= x; len <<= 1);
    return len;
}
inline void NTT(int  *A, int type, int len) {
    int limit = 1, L = 0;
    while(limit < len) limit <<= 1, L++;
    for(rg int i = 0; i < limit; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));    
    for(rg int i = 0; i < limit; i++) if(i < r[i]) swap(A[i], A[r[i]]);
    for(rg int mid = 1; mid < limit; mid <<= 1) {    
        int  R = mid << 1;
        int  W = fastpow(Gi, (P - 1) / R); Og[0] = 1;
        for(rg int j = 1; j < mid; j++) Og[j] = mul(Og[j - 1], W);
        for(rg int j = 0; j < limit; j += R) {
            for(rg int k = 0; k < mid; k++) {
                 const int x = A[j + k], y = mul(Og[k], A[j + k + mid]);
                 A[j + k] = add(x, y), A[j + k + mid] = dec(x, y);
            }
        }
    }
    if(type == -1) {
        std::reverse(&A[1], &A[limit]);
        for(int i = 0, inv = fastpow(len , P - 2); i < limit; i++)
            A[i] = 1ll * A[i] * inv % P;
    }
}
void Inv(int *a, int *b, int len) {// a要求的多项式 b逆元 len：要求的逆元的长度 
    if(len == 1) {
        b[0] = fastpow(a[0], P - 2);
        return ;
    }
    Inv(a, b, len >> 1);
    for(rg int i = 0; i < len; i++) A[i] = a[i], B[i] = b[i];
    NTT(A, 1, len << 1); 
    NTT(B, 1, len << 1);
    for(rg int i = 0; i < (len << 1); i++) A[i] = mul(mul(A[i], B[i]), B[i]) ;
    NTT(A, -1, len << 1);
    for(rg int i = 0; i < len; i++) b[i] = (1ll * (b[i] << 1) % P + P - A[i] ) % P;
    for(rg int i = 0; i < (len << 1); i++) A[i] = B[i] = 0;
}
void Mul(int *a, int *b, int *c, int N, int M) {
    int len = GetLen(max(N, M)) << 1;
    for(int i = 0; i <= N; i++) Atmp[i] = a[i];
    for(int i = 0; i <= M; i++) Btmp[i] = b[i];
    NTT(Atmp, 1, len); NTT(Btmp, 1, len);
    for(int i = 0; i <= len; i++) c[i] = 1ll * Atmp[i] * Btmp[i] % P, Atmp[i] = Btmp[i] = 0;
    NTT(c, -1, len);
}
int main() {
    #ifdef WIN32
    freopen("a.in","r",stdin);
    #endif 
       int N = read(), M = read();
       for(int i = 0; i <= N; i++) F[i] = read();
       for(int i = 0; i <= M; i++) G[i] = read();
       reverse(F, F + N + 1); reverse(G, G + M + 1);
       //for(int i = 0; i <= N - M; i++) Ginv[i] = G[i];//tag
       Inv(G, Ginv, GetLen(N - M));
       Mul(F, Ginv, Q, N - M, N - M);
    reverse(Q, Q + N - M + 1);
       for(int i = 0; i <= N - M; i++) print(Q[i]), *O++ = ' '; *O++ = '\n';
       reverse(F, F + N + 1); reverse(G, G + M + 1);
       Mul(Q, G, R, N - M, M);
       for(int i = 0; i < M; i++) 
       print(dec(F[i], R[i])), *O++ = ' ';
    fwrite(obuf, O-obuf, 1 , stdout);
    return 0;
}

泰勒展开

泰勒公式是一个用函数在某点的信息描述其附近取值的公式。若是函数足够平滑的话，在已知函数在某一点的各阶导数值的状况之下，泰勒公式能够用这些导数值作系数构建一个多项式来近似函数在这一点的邻域中的值。泰勒公式还给出了这个多项式和实际的函数值之间的误差。

一句话应用：构造多项式来逼近函数

$$g(x)=g(0)+\frac{f^{1}(0)}{1!}x+\frac{f^{2}(0)}{2!}x^{2}+……+\frac{f^{n}(0)}{n!}x^{n}$$

$f^n$表示对$f$进行$n$次求导

这里的“多项式”咱们能够直观的理解为一种特殊的“函数”

普通牛顿迭代法

用途：求函数$f(x)$的零点

首先任取一个点$x_0$

而后对$f(x)$泰勒展开

$$f(x)=f(x_0)+f'(x_0)(x-x_0)+\frac {f''(x_0)(x-x_0)^2}2+\cdots$$

咱们只保留线性部分

$$f(x)=f(x_0)+f'(x_0)(x-x0)=0$$

这样的话

$$x=x0-\frac{f(x_0)}{f'(x_0)}$$

而后用$x$替换$x_0$，索然咱们在保留线性部分时丢掉了不少信息，可是咱们进行屡次迭代仍然能够获得最终解

多项式牛顿迭代法

用途：

已知多项式$F(x)$和函数$G(x)$，求

$$G(F(x)) \equiv 0 \pmod {x^n}$$

算法

仍然考虑倍增

当$n=1$时，$F(x)$仅有一个常数项，

上面的式子能够化为

$$G(x) \equiv 0$$

咱们能够直接利用牛顿迭代求解

按照多项式求逆的思路，假设咱们已经求得

$$G(F_0(x)) \equiv 0 \pmod {x^{\lceil \frac{n}{2} \rceil}}$$

而后在$F_0(x)$处进行泰勒展开

$$\begin{eqnarray*} G(F(x)) & = & G(F_0(x)) \\ & + & \frac{G'(F_0(x))}{1!}\left(F(x) - F_0(x)\right) \\ & + & \frac{G''(F_0(x))}{2!}\left(F(x) - F_0(x)\right)^2 \\ & + & \cdots \end{eqnarray*}$$

在这里咱们不难发现一个问题$F_0(x)$的前$\frac{n}{2}$项与$F(x)$是相同的

那么$F(x)-F_0(x)$的前$\frac{n}{2}$项必然是$0$

那么对于$(F(x)-F_0)^2$，前$n$项必然是$0$

同理当指数大于$2$时，前$n$项必然也是$0$

那么咱们若是只保留整数部分，获得的应该是准确解！

$$G(F(x)) \equiv G(F_0(x)) + G'(F_0(x))\left(F(x) - F_0(x)\right) \pmod {x^n}$$

因为$$G(F(x)) \equiv 0 \pmod {x^n}$$

而后就作完了

$$F(x) \equiv F_0(x) - \frac{G(F_0(x))}{G'(F_0(x))} \pmod {x^n}$$

多项式开根

利用牛顿迭代法能够快速的推出多项式开根的作法

多项式开根即已知多项式$A(x)$，求多项式$B(x)$，知足

$B^2(x) \equiv A(x) \pmod{x^n}$

设$F(x)$知足

$F^2(x) - A(x) \equiv 0 \pmod{x^n}$

设$G(F(x)) = F(x)^2-A(x)$

咱们要求的也就是$G(F(x))$的零点

根据求导公式，获得$G'(F(x)) = 2F(x)$(在这里咱们认为$A(x)$为常数)

$$\begin{eqnarray*} F(x) & \equiv & F_0(x) - \frac{F_0^2(x) - A(x)}{2F_0(x)} \\ & \equiv & \frac{F_0^2(x) + A(x)}{2F_0(x)} \pmod {x^n} \end{eqnarray*}$$

$F_0$是在$\pmod {x^{\frac{n}{2}}}$意义下的结果

这里还有一个问题，就是当多项式仅有一个常数项的时候怎么处理

咱们实际是要计算

$$\begin{equation}\label{quadratic}x^2 \equiv a \pmod p~~(p \nmid a)\end{equation}$$

这玩意儿的学名叫作二次剩余

能够看这里https://blog.csdn.net/a_crazy_czy/article/details/51959546

多项式对数与多项式$exp$

麦克劳林级数

在泰勒公式中，取$x=0$，获得的级数

$$\sum ^{\infty }_{n=0}\dfrac {f^{n}\left( 0\right) }{n!}x^{n}$$

称为麦克劳林级数。

对数的计算

其实我并不知道这玩意儿的真正意义是什么，不过有大佬给出了它的定义，那就默认按定义来吧

多项式的对数能够认为是一个多项式和麦克劳林级数的复合

即对于多项式$A(x)$，有

$$\ln (1 - A(x)) = -\sum_{i \geq 1} \frac{A^i(x)}{i}$$

在计算时直接对$lnA(x)$微分后再积分

$$lnA(x)=\int (lnA(x))'=\int\frac{A'(x)}{A(x)}$$

这里用到了复合函数求导的链式法则

$f(g(x))=f’(g(x))g’(x)$

这样在计算的时候能够先求导再积分，这二者的复杂度都是$O(n)$的

加上多项式求逆的复杂度，总复杂度为$O(nlogn)$

多项式exp

多项式的定义为，给定多项式$A(x)$

$$\exp(A(x)) = e^{A(x)} = \sum_{i \geq 0} \frac{A^i(x)}{i!}$$

计算方法：

设$F(x) = e^{A(x)}$

$$lnF(x) = A(x)$$

设$G(F(x)) = lnF(x) - A(x) = 0$

这样就转换成了求多项式零点的问题，直接上牛顿迭代

$$\begin{eqnarray*} F(x) & \equiv & F_0(x) - \frac{G(F_0(x))}{G'(F_0(x))} \\ & \equiv & F_0(x) \left (1 - \ln F_0(x) + A(x) \right ) \pmod {x^n}\end{eqnarray*}$$

复杂度$O(nlogn)$

多项式幂函数

给定多项式$A(x)$和$k$，求

$$A^k(x)\pmod{x^n}$$

用上面的两种方法转换

$$A^k(x)\equiv e^{klnA(x)}\pmod{x^n}$$

时间复杂度$O(nlogn)$，常数巨大！

多项式三角函数

不会

$$e^{iA(x)}=cos(A(x))+isin(A(x))$$

多项式的快速差值与多点求值

此部分参(抄)考(袭)自http://blog.miskcoo.com/2015/05/polynomial-multipoint-eval-and-interpolation

多点求值

给出一个多项式$A(x)$以及$n$个点$x_0, x_1, x_2, \dots, x_{n-1}$

求出$A(x_0), A(x_1), A(x_2),\dots,x_{n-1}$，即每一个点在多项式中的取值

快速差值

给出一个集合$$X = \{(x_i, y_i) : 0 \leq i \leq n\}$$

求一个多项式$A(x)$，知足

$$\forall (x, y) \in X, A(x) = y$$

算法

由于这两个问题的特殊性，所以在计算过程当中可能会用到彼此，你们直接略过就好

多点求值

首先把要求的值分红两半

$$\begin{eqnarray*} X^{[0]} &=& \{x_0, x_1, \cdots, x_{\lfloor \frac{n}{2} \rfloor}\} \\ X^{[1]} &=& \{x_{\lfloor \frac{n}{2} \rfloor+1},x_{\lfloor \frac{n}{2} \rfloor+2}, \cdots, x_{n-1}\} \end{eqnarray*}$$

咱们拿前一半来说，后一半同理

记用$x^[0]$插值获得的$\frac{n}{2}$次多项式为$A^[0]$

构造多项式

$$\begin{eqnarray*} P^{[0]}(x) &=& \prod_{i=0}^{\lfloor \frac{n}{2} \rfloor} (x-x_i) \end{eqnarray*}$$

对于该多项式来讲，当$x \in x^[0]$时，$P^{[0]}(x) = 0$，那么$A(x)$能够表示为

$$A(x) = D(x)P^{[0]}(x) + A^{[0]}(x)$$

这个式子等价于

$$A^{[0]}(x) \equiv A(x) \pmod {P^{[0]}(x)}$$

用多项式取模算就能够

时间复杂度：$O(nlog^n)$

快速插值

$O(N^2)$：

$${A(x)=\sum_{i=1}^n{\frac{\prod_{{j}\neq{i}}{({x}-{x_j})}}{\prod_{{j}\neq{i}}{({x_i}-{x_j})}}{y_i}}}$$

首先按下标分类

$$\begin{eqnarray*} X^{[0]} &=& \{(x_i, y_i) : 0 \leq i \leq \lfloor \frac{n}{2} \rfloor \} \\ X^{[1]} &=& \{(x_i, y_i) : \lfloor \frac{n}{2} \rfloor < i \leq n\} \end{eqnarray*}$$

如今假设已经求出了$X^{[0]}$的插值多项式$A^{[0]}(x)$，
设
$$P(x) = \prod_{i=0}^{\lfloor \frac{n}{2} \rfloor} (x-x_i)$$

$$A(x) = A^{[1]}(x)P(x) + A^{[0]}(x)$$

其中 $A^{[1]}(x)$ 是一个未知的多项式，因为 $\forall (x, y) \in X^{[0]}, P(x) = 0, A^{[0]}(x) = y$
这样的话就知足 $X^{[0]}$的点都在$A(x)$上，问题就变成要将$X^{[1]}$ 内的点插值，使得
$\forall (x_i, y_i) \in X^{[1]}, y_i = A^{[1]}(x_i)P(x_i) + A^{[0]}(x_i)$化简以后获得

$$A^{[1]}(x_i) = \frac{A^{[0]}(x_i)-y_i}{P(x_i)}$$

这样就获得了新的待插值点，利用一样的方法求出插值出$A^{[1]}$而后合并就能够了

因为每一次都要多点求值求出$X^{[1]}$内$P(x)$和$A^{[0]}(x)$的值，最终复杂度是

$T(n) = 2T(\frac{n}{2}) + \mathcal O(n \log^2 n) = \mathcal O(n \log^3 n)$