首页 > > 程序设计 > C/C++ >

多项式整理

2018-06-17 20:29:57来源：未知阅读 ()

多项式求逆元

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

使得

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

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

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

至于我们为什么要模$x^n$，因为通过计算不难发现：除了仅有常数项的多项式的逆元为一个常数之外，其余多项式的逆元均有无穷多项

算法

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

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

即

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

根据取模运算的性质，$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}$$

由于$$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)$

标签：

版权申明：本站文章部分自网络，如有侵权，请联系：west999com@outlook.com
特别注意：本站所有转载文章言论不代表本站观点，本站所提供的摄影照片，插画，设计作品，如需使用，请与原作者联系，版权归原作者所有