抱歉，您的浏览器无法访问本站

本页面需要浏览器支持（启用）JavaScript


了解详情 >

AutumnKite's Blog

博客

归档

友链

关于

说说

留言

多项式相关

This is a picture without description

AutumnKite
笔记

发布于：2019年12月7日
次浏览

太久没更新博客，开一篇多项式相关的博客记一下多项式的一些基础操作吧。

模域下的运算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
const int P = 998244353, G = 3;

void inc(int &a, int b) {
    (a += b) >= P ? a -= P : 0;
}

void dec(int &a, int b) {
    a < b ? a += P - b : a -= b;
}

int plus(int a, int b) {
    return a + b >= P ? a + b - P : a + b;
}

int minus(int a, int b) {
    return a < b ? a + P - b : a - b;
}

int mul(int a, int b) {
    return 1ll * a * b % P;
}

int qpow(int a, int b) {
    int s = 1;
    for (; b; b >>= 1) {
        if (b & 1) {
            s = mul(s, a);
        }
        a = mul(a, a);
    }
    return s;
}

快速傅里叶变换与快速傅里叶逆变换

快速傅里叶变换

在这一节中，除特殊说明外，\(n\) 是 \(2\) 的幂次。

利用 FFT，我们可以在 \(O(n\log n)\) 的时间内求出 \(F(\omega_n^0),F(\omega_n^1),\cdots,F(\omega_n^{n-1})\) 的值，以及由这些值得到 \(F(x)\) 的系数，其中 \(\omega_n^{i}\) 表示单位根，这里不具体展开。

单位根具有一些优美的性质，如

  • \(\omega_{dn}^{dk}=\omega_n^k\)
  • \(\omega_{n}^{\frac{n}{2}}=\omega_2^1=-1\)
  • \(\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}\)
  • \(\omega_n^a=\omega_n^{a \bmod n}\)

FFT 正是利用了这些性质加快速度。

假设原多项式为 \(F(x)=\sum\limits_{i=0}^{n-1} f_ix^i\)，我们把 \(F(x)\) 按指数奇偶性拆成两部分，即 \[F_0(x)=\sum\limits_{i=0}^{\frac{n}{2}-1} f_{2i}x^i,F_1(x)=\sum\limits_{i=0}^{\frac{n}{2}-1} f_{2i+1}x^i\] 那么有 \[F(x)=F_0(x^2)+xF_1(x^2)\]

我们把 \(x=\omega_n^k\ (0\le k < \frac{n}{2})\) 代入 \(F(x)\)，得到 \[F(\omega_n^k)=F_0(\omega_{\frac{n}{2}}^k)+\omega_n^kF_1(\omega_{\frac{n}{2}}^k)\]

同理，把 \(x=\omega_n^{k+\frac{n}{2}}\ (0\le k < \frac{n}{2})\) 代入 \(F(x)\)，得到 \[F(\omega_n^k)=F_0(\omega_{\frac{n}{2}}^k)-\omega_n^kF_1(\omega_{\frac{n}{2}}^k)\]

发现问题转化成了两个子问题，直接分治求解即可。时间复杂度 \(O(n\log n)\)。

直接做常数太大，我们发现分治到底层时下标是原下标的二进制翻转，于是直接翻转完后从下往上合并即可。

可以用 unsigned long long 减少取模次数优化常数。

快速傅里叶逆变换

记多项式 \(G(x)=\sum\limits_{i=0}^{n-1}F(\omega_n^i)x^i\)，考虑对 \(G(x)\) 做 FFT，得到

\[\begin{aligned}G(\omega_n^k)&=\sum\limits_{i=0}^{n-1}F(\omega_n^i)(\omega_n^k)^i\\&=\sum\limits_{i=0}^{n-1}\left(\sum\limits_{j=0}^{n-1} f_j(\omega_n^i)^j\right)(\omega_n^k)^i\\&=\sum\limits_{j=0}^{n-1}f_j\sum\limits_{i=0}^{n-1}(\omega_n^{j+k})^i\end{aligned}\]

发现只有当 \(j+k\equiv 0\pmod n\) 时 \(\sum\limits_{i=0}^{n-1}(\omega_n^{j+k})^i=n\)，其他情况都为 \(0\)。于是

\[\begin{aligned}G(\omega_n^k)&=\sum\limits_{j=0}^{n-1}f_j\sum\limits_{i=0}^{n-1}(\omega_n^{j+k})^i\\&=f_{(n-k)\bmod n}\cdot n\end{aligned}\]

直接翻转后乘上 \(n\) 的逆元即可，时间复杂度同样为 \(O(n\log n)\)。

实现

\(\omega_{2^k}^{i}\) 只有 \(O(n)\) 个，可以在程序的一开始预处理（同时还预处理了逆元）。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
const int MAX_LEN = 1 << 18;

typedef std::vector<int> poly;

int omega[MAX_LEN], inv[MAX_LEN];

void Init() { // 在程序开始调用
    for (int m = 1; m < MAX_LEN; m <<= 1) {
        int w = qpow(G, (P - 1) / (m << 1));
        omega[m] = 1;
        for (int i = 1; i < m; ++i) {
            omega[m + i] = mul(omega[m + i - 1], w);
        }
    }
    inv[1] = 1;
    for (int i = 2; i < MAX_LEN; ++i) {
        inv[i] = mul(P - P / i, inv[P % i]);
    }
}

每次需要预处理 \(rev[i]\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
poly rev;

int get(int _n) {
    int n = 1;
    while (n < _n) {
        n <<= 1;
    }
    return n;
}

void init(int n) {
    int k = 0;
    while ((1 << k) < n) {
        ++k;
    }
    rev.resize(n), rev[0] = 0, --k;
    for (int i = 1; i < n; ++i) {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << k);
    }
}

void DFT(poly &f, int n) {
    f.resize(n);
    static unsigned long long F[MAX_LEN];
    for (int i = 0; i < n; ++i) {
        F[rev[i]] = f[i];
    }
    for (int m = 1; m < n; m <<= 1) {
        // 注意当 MAX_LEN > 2^18 时这里需要把 F[i] 全部模 P
        for (int p = 0, l = m << 1; p < n; p += l) {
            int *W = omega + m;
            unsigned long long *F0 = F + p, *F1 = F0 + m;
            for (int i = 0; i < m; ++i, ++F0, ++F1, ++W) {
                int t = (*F1) * (*W) % P;
                *F1 = *F0 + P - t, *F0 += t;
            }
        }
    }
    for (int i = 0; i < n; ++i) {
        f[i] = (F[i] % P + P) % P;
    }
}

void IDFT(poly &f, int n) {
    DFT(f, n), std::reverse(f.begin() + 1, f.end());
    int t = qpow(n, P - 2);
    for (int i = 0; i < n; ++i) {
        f[i] = mul(f[i], t);
    }
}

多项式加法与减法

对于每个 \(x_i\) 的系数直接相加/相减即可。

时间复杂度 \(O(n)\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
poly Plus(poly a, poly b) {
    int n = std::max(a.size(), b.size());
    a.resize(n), b.resize(n);
    for (int i = 0; i < n; ++i) {
        inc(a[i], b[i]);
    }
    return a;
}

poly Minus(poly a, poly b) {
    int n = std::max(a.size(), b.size());
    a.resize(n), b.resize(n);
    for (int i = 0; i < n; ++i) {
        dec(a[i], b[i]);
    }
    return a;
}

多项式乘法

DFT 后把点值乘起来，在 IDFT 回去即可。

代码中的参数 _n 的作用在下面会提到。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
poly Mul(poly a, poly b, int _n = -1) {
    int na = (int)a.size(), nb = (int)b.size(), n = _n;
    if (!na || !nb || !n) {
        return {};
    }
    if (n == -1) {
        n = na + nb - 1;
    }
    _n = n, n = get(n);
    init(n), DFT(a, n), DFT(b, n);
    for (int i = 0; i < n; ++i) {
        a[i] = mul(a[i], b[i]);
    }
    IDFT(a, n);
    a.resize(_n);
    return a;
}

可以把计算 \(F^2(x)\) 单独写一个函数，常数可以更小。

1
2
3
4
5
6
7
8
9
10
11
poly Pow2(poly a) {
    int _n = 2 * (int)a.size() - 1;
    int n = get(_n);
    init(n), DFT(a, n);
    for (int i = 0; i < n; ++i) {
        a[i] = mul(a[i], a[i]);
    }
    IDFT(a, n);
    a.resize(_n);
    return a;
}

多项式求逆

已知多项式 \(F(x)\)，求一个多项式 \(G(x)\) 满足 \(F(x)G(x)\equiv 1\pmod{x^n}\)。

假设我们已经求出 \(G_0(x)\) 满足 \[F(x)G_0(x)\equiv 1\pmod {x^{\lceil\frac{n}{2}\rceil}}\] 那么有 \[F(x)(G(x)-G_0(x))\equiv 0\pmod{x^{\lceil\frac{n}{2}\rceil}}\]

由于 \(G(x)-G_0(x)\) 是 \(n-1\) 次多项式，两边同时乘 \(G(x)-G_0(x)\) 得到 \[F(x)(G(x)-G_0(x))^2\equiv 0\pmod{x^n}\]

展开后整理得到 \[G(x)\equiv 2G_0(x)-F(x)G_0^2(x)\pmod {x^n}\]

迭代求出即可。

边界情况，即 \(n=1\) 时，求 \(f_0\) 的逆元即可。

时间复杂度 \(O(n\log n)\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
poly Inv(poly f, int _n) {
    int n = get(_n);
    f.resize(n);
    poly g(1);
    g[0] = qpow(f[0], P - 2);
    for (int m = 2; m <= n; m <<= 1) {
        int l = m << 1;
        poly tmp(f.begin(), f.begin() + m);
        init(l), DFT(tmp, l), DFT(g, l);
        for (int i = 0; i < l; ++i) {
            g[i] = (2 + 1ll * (P - tmp[i]) * g[i]) % P * g[i] % P;
        }
        IDFT(g, l), g.resize(m);
    }
    g.resize(_n);
    return g;
}

多项式除法

已知 \(n-1\) 次多项式 \(A(x)\) 和 \(m-1\) 次多项式 \(B(x)\)，求 \(n-m\) 次多项式 \(Q(x)\) 和小于 \(m-1\) 次的多项式 \(R(x)\) 满足 \(A(x)=B(x)Q(x)+R(x)\)。

下面我们强制 \(R(x)\) 为 \(m-2\) 次，不足则高位补 \(0\)。

用记号 \(F_r(x)\) 表示 \(F(x)\) 的系数翻转后的结果。有 \[F_r(x)=x^{\operatorname{deg}F(x)}F(\frac{1}{x})\]

其中 \(\operatorname{deg}F(x)\) 表示 \(F(x)\) 的次数。

把原式中的 \(x\) 用 \(\frac{1}{x}\) 代替，得到 \[A(\frac{1}{x})=B(\frac{1}{x})Q(\frac{1}{x})+R(\frac{1}{x})\] 两边同乘 \(x^{n-1}\) 得到 \[x^{n-1}A(\frac{1}{x})=x^{m-1}B(\frac{1}{x})x^{n-m}Q(\frac{1}{x})+x^{n-m+1}\cdot x^{m-2}R(x)\] 即 \[A_r(x)=B_r(x)Q_r(x)+x^{n-m+1}R(x)\] 那么有 \[A_r(x)\equiv B_r(x)Q_r(x)\pmod{x^{n-m+1}}\]

由于 \(Q(x)\) 是一个 \(n-m\) 次多项式，所以 \(\bmod x^{n-m+1}\) 不会对 \(Q(x)\) 产生影响。

此时由于去掉了 \(R(x)\)，直接做多项式求逆即可求出 \(Q(x)\)。然后根据 \(R(x)=A(x)-B(x)Q(x)\) 求出 \(R(x)\)。

时间复杂度 \(O(n\log n)\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
std::pair<poly, poly> Div(const poly &a, const poly &b) {
    int n = a.size(), m = b.size();
    poly Q, R;
    if (n < m) {
        R = a, R.resize(m - 1);
        return {Q, R};
    }
    poly A(a), B(b);
    std::reverse(A.begin(), A.end()), A.resize(n - m + 1);
    std::reverse(B.begin(), B.end()), B.resize(n - m + 1);
    Q = Mul(A, Inv(B, n - m + 1));
    Q.resize(n - m + 1), std::reverse(Q.begin(), Q.end());
    R = Minus(a, Mul(b, Q)), R.resize(m - 1);
    return {Q, R};
}

多项式求导与积分

求导有运算法则： - \((u(x)+v(x))'=u'(x)+v'(x)\) - \((u(x)\cdot v(x))'=u'(x)\cdot v(x)+u(x)\cdot v'(x)\)

并且可以扩展到 \(n\) 个函数。同时有 \(C'=0,(x^n)'=nx^{n-1}\)。

于是 \[F'(x)=\left(\sum\limits_{i=0}^{n-1} f_ix^i\right)'=\sum\limits_{i=0}^{n-1}(f_ix^i)'=\sum\limits_{i=1}^{n-1} f_i\cdot ix^{i-1}\]

积分只要反着做一遍就好了。

时间复杂度 \(O(n)\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
poly Der(poly f) {
    int n = (int)f.size();
    for (int i = 1; i < n; ++i) {
        f[i - 1] = mul(f[i], i);
    }
    f.pop_back();
    return f;
}

poly Int(poly f) {
    f.push_back(0);
    int n = f.size();
    for (int i = n - 1; i; --i) {
        f[i] = mul(f[i - 1], inv[i]);
    }
    f[0] = 0;
    return f;
}

泰勒展开

泰勒展开是将一个在 \(x=x_0\) 处具有 \(n\) 阶导数的函数 \(f(x)\) 利用关于 \(x-x_0\) 的 \(n\) 次多项式来逼近函数的方法。

当 \(n\) 趋向于无穷时，我们可以直接用下面这个式子表示泰勒公式：

\[f(x)=\sum\limits_{i=0}^{\infty} \frac{f^{(i)}(x_0)}{i!}(x-x_0)^i\]

而麦克劳林级数就是上式 \(x_0=0\) 的特殊形式。

要对泰勒展开有更好地理解可以看怎样更好地理解并记忆泰勒展开式？- 知乎。

牛顿迭代

牛顿迭代法是一种用于解非线性方程 \(G(x)=0\) 的近似方法。

而在我们研究多项式相关时，我们一般用牛顿迭代法解关于无穷幂级数 \(F(x)\) 的非线性方程 \(G(F(x))=0\)，求 \(F(x)\) 在模 \(x^n\) 意义下得到的多项式，其中 \(G(F(x))\) 是一个自变量为多项式，应变量也为多项式的函数。

我们假设已经求出答案模 \(x^{\lceil\frac{n}{2}\rceil}\) 意义下的多项式 \(F_0(x)\)，我们要求出答案模 \(x^n\) 意义下的多项式 \(F_1(x)\)。即我们已经有 \[G(F_0(x))\equiv 0\pmod{x^{\lceil\frac{n}{2}\rceil}}\]

我们把 \(G(F(x))\) 在 \(F(x)=F_0(x)\) 处泰勒展开，得到 \[\begin{aligned}G(F(x))&=G(F_0(x))\\ &+G'(F_0(x))(F(x)-F_0(x))\\ &+\frac{G''(F_0(x))}{2}(F(x)-F_0(x))^2\\&+\cdots\end{aligned}\]

因为有 \(F_1(x)-F_0(x)\equiv 0\pmod{x^{\lceil\frac{n}{2}\rceil}}\)，所以有 \((F_1(x)-F_0(x))^2\equiv 0\pmod{x^n}\)。

由于 \(F_1(x)\) 是答案在模 \(x^n\) 意义下的多项式，所以上述泰勒展开式中除前两项外的部分代入 \(F_1(x)\) 后模 \(x^n\) 为 \(0\)。所以我们只需要考虑前两项。即 \(F_1(x)\) 只需要满足 \[G(F_0(x))+G'(F_0(x))(F_1(x)-F_0(x))\equiv 0\pmod{x^n}\]

整理，得 \[F_1(x)\equiv F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}\pmod{x^n}\]

于是我们可以在 \(O(n\log n)\) 的时间求出答案在模 \(x^n\) 意义下的多项式。

多项式对数函数

已知多项式 \(F(x)\)（常数项为 \(1\)），求多项式 \(G(x)\) 满足 \(G(x)\equiv \ln F(x)\pmod{x^n}\)。

注意多项式的对数函数可以认为是多项式和麦克劳林级数的复合。这意味着多项式对数函数是一个无穷级数。我们只能求他在模 \(x^n\) 意义下得到的多项式。

类似地，多项式求逆、多项式指数函数、多项式开根等我们都只能求出在模意义下的多项式。

注意 \(F(x)\) 常数项不为 \(1\) 时 \(\ln F(x)\) 无意义。

两边同时求导，得到 \[G'(x)\equiv (\ln F(x))'\pmod {x^{n-1}}\] 由复合函数求导公式，得到 \[G'(x)\equiv \frac{F'(x)}{F(x)}\pmod {x^{n-1}}\] 然后两边同时求积分，得到 \[\ln F(x)\equiv \int \frac{F'(x)}{F(x)}\pmod {x^n}\]

时间复杂度 \(O(n\log n)\)，瓶颈在于多项式乘法和多项式求逆。

1
2
3
4
poly Ln(const poly &a, int n) {
    poly res = Int(Mul(Der(a), Inv(a, n)));
    return res.resize(n), res;
}

多项式指数函数

已知多项式 \(F(x)\)（常数项为 \(0\)），求多项式 \(G(x)\) 满足 \(G(x)\equiv \exp {F(x)}\pmod{x^n}\)。

多项式 \(\exp\) 有两种做法，分别是 \(O(n\log n)\) 的牛顿迭代法与 \(O(n\log^2 n)\) 的分治法，并且分治法可以进一步优化成 \(O(\frac{n\log^2 n}{\log\log n})\)。从实际运行效率上看，\(n\le 10^5\) 时普通分治与牛顿迭代的运行速度不相上下，且分治法的实现更简单。

牛顿迭代法

注意 \(F(x)\) 常数项不为 \(0\) 时 \(\exp F(x)\) 无意义。

两边同时求 \(\ln\)，得到 \[\ln G(x)\equiv F(x)\pmod{x^n}\] 即 \[(\ln G(x))-F(x)\equiv 0\pmod{x^n}\]

定义函数 \[H(G(x))=(\ln G(x))-F(x)\] 那么我们要求 \(G(x)\) 满足 \(H(G(x))\equiv 0\pmod{x^n}\)。

套用牛顿迭代的式子，即假设我们求出了在模 \(x^{\lceil\frac{n}{2}\rceil}\) 意义下的答案 \(G_0(x)\)，那么有 \[G(x)\equiv G_0(x)-\frac{H(G_0(x))}{H'(G_0(x))}\pmod{x^n}\]

又因为有（要注意函数 \(H(G(x))\) 的自变量是 \(G(x)\) 而非 \(x\)） \[\begin{aligned}H'(G(x))&=\ln' G(x)\\ &=\frac{1}{G(x)}\end{aligned}\]

所以，经整理得 \[G(x)\equiv G_0(x)(1-(\ln G_0(x))+F(x))\pmod{x^n}\]

时间复杂度 \(O(n\log n)\)。常数较大。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
poly ExpNewton(poly f, int _n) {
    int n = get(_n);
    f.resize(n);
    poly g(1);
    g[0] = 1;
    for (int m = 2; m <= n; m <<= 1) {
        int l = m << 1;
        poly tmp = Ln(g, m);
        for (int i = 0; i < m; ++i) {
            tmp[i] = minus(f[i], tmp[i]);
        }
        inc(tmp[0], 1);
        init(l), DFT(tmp, l), DFT(g, l);
        for (int i = 0; i < l; ++i) {
            g[i] = mul(g[i], tmp[i]);
        }
        IDFT(g, l), g.resize(m);
    }
    return g.resize(_n), g;
}

分治法

我们在 \(G(x)\equiv \exp F(x)\pmod{x^n}\) 两边同时求导，得到 \[G'(x)\equiv (\exp F(x))F'(x)\pmod{x^n}\]

即 \[G'(x)\equiv G(x)F'(x)\pmod{x^n}\]

两边同时求积分，得 \[G(x)\equiv \int G(x)F'(x)\]

直接分治 NTT 求解即可。时间复杂度 \(O(n\log^2 n)\)。

关于卡常

在处理区间 \([l,r)\) 时我们需要将 \(f[0,r-l)\) 与 \(g[l,\lfloor\frac{l+r+1}{2}\rfloor)\) 作卷积，求积分以后再加到 \(g[\lfloor\frac{l+r+1}{2}\rfloor,r)\) 上（其中 \(f,g\) 分别是 \(F'(x),G(x)\) 的系数数组）。

发现我们真正需要的部分是卷积后中间一部分，前后两部分是没有用的。而通过计算发现后面部分的长度一定不超过前面部分的长度，FFT 是循环卷积，我们只要把长度限定到 \(r-l-1\) 就可以把后面部分加到前面而不影响中间需要的部分。这样做对优化常数很有帮助。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
void ExpSolve(const poly &f, poly &g, int l, int r) {
    if (l + 1 == r) {
        if (l) {
            g[l] = mul(g[l], inv[l]);
        } else {
            g[l] = 1;
        }
        return;
    }
    int md = (l + r + 1) >> 1;
    ExpSolve(f, g, l, md);
    poly tmp(Mul(poly(g.begin() + l, g.begin() + md), 
	  poly(f.begin(), f.begin() + r - l - 1), r - l - 1));
    for (int i = md; i < r; ++i) {
        inc(g[i], tmp[i - l - 1]);
    }
    ExpSolve(f, g, md, r);
}

poly Exp(poly f, int _n) {
    f.resize(_n), f = Der(f);
    poly g(_n);
    ExpSolve(f, g, 0, _n);
    return g;
}

更为优秀的分治做法

咕咕咕

多项式求任意次幂

给定多项式 \(F(x)\) 和正整数 \(k\)，求 \(G(x)\equiv F^k(x)\pmod{x^n}\)。

如果 \(F(x)\) 常数项为 \(1\)，那么有 \(G(x)=e^{k\ln F(x)}\)。

而 \(F(x)\) 常数项不为 \(1\) 时，我们可以把 \(F(x)\) 除以 \(ax^p\) 使得常数项为 \(1\)，最后再乘上 \(a^kx^{kp}\)。

理论上来说这个方法可以求任意实数幂次，但是系数在模意义下运算，\(k\) 为分数、无理数时很难求或者在模域下不存在答案。

时间复杂度 \(O(n\log n)\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
poly Pow(poly a, int k, int n) {
    a.resize(n);
    int t = n;
    for (int i = 0; i < n; ++i) {
        if (a[i]) {
            t = i;
            break;
        }
    }
    if (t == n) {
        if (!k) {
            a[0] = 1;
        }
        return a;
    }
    int vi = qpow(a[t], P - 2), vk = qpow(a[t], k);
    for (int i = 0; i < n - t; ++i) {
        a[i] = mul(a[i + t], vi);
    }
    a.resize(n - t);
    t = std::min(1ll * t * k, 1ll * n);
    if (t == n) {
        a.resize(n);
        for (int i = 0; i < n; ++i) {
            a[i] = 0;
        }
        return a;
    }
    a = Ln(a, n - t);
    for (int i = 0; i < n - t; ++i) {
        a[i] = mul(a[i], k);
    }
    a = Exp(a, n - t), a.resize(n);
    for (int i = n - t - 1; ~i; --i) {
        a[i + t] = mul(a[i], vk);
    }
    for (int i = 0; i < t; ++i) {
        a[i] = 0;
    }
    return a;
}

注意下面部分的代码还未更新。

多项式开根

已知多项式 \(F(x)\)，求多项式 \(G(x)\) 满足 \(G^2(x)\equiv F(x)\pmod{x^n}\)。

定义函数 \(H(G(x))=G^2(x)-F(x)\)，那么我们要求 \(G(x)\) 使得 \(H(G(x))\equiv 0\pmod{x^n}\)。

仍然假设我们已经求出了 \(G_0(x)\) 满足 \(H(G_0(x))\equiv 0\pmod{x^{\lceil\frac{n}{2}\rceil}}\)。根据牛顿迭代的式子，得到 \[\begin{aligned}G(x)&\equiv G_0(x)-\frac{H(G_0(x))}{H'(G_0(x))}\\&=G_0(x)-\frac{G_0^2(x)-F(x)}{2G_0(x)}\\&=\frac{G_0^2(x)+F(x)}{2G_0(x)}\pmod{x^n}\end{aligned}\]

在模意义下运算时，边界我们只要用 Cipolla 算法求二次剩余即可。

时间复杂度 \(O(n\log n)\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
poly Sqrt(poly f, int _n){
    int n = get(_n);
    f.resize(n);
    poly g(1);
    g[0] = Sqrt(f[0]); // 求二次剩余
    for (register int m = 2; m <= n; m <<= 1){
        poly tmp = Pow2(g);
        tmp.resize(m);
        for (register int i = 0; i < m; ++i) inc(tmp[i], f[i]);
        g = Multiply(tmp, Inverse(g, m)), g.resize(m);
        for (register int i = 0; i < m; ++i) g[i] = mul(g[i], inv[2]);
    }
    return g.resize(_n), g;
}

多项式多点求值

给定多项式 \(F(x)\) 和 \(n\) 个点 \(a_i\)，对于每个 \(a_i\)，求 \(F(a_i)\) 的值。

记 \(k=\lfloor\frac{n}{2}\rfloor\)，构造多项式 \(G_0(x)=\prod\limits_{i=1}^{k}(x-a_i)\)，然后用多项式除法求出 \(D_0(x),R_0(x)\) 满足 \(F(x)=D_0(x)G_0(x)+R_0(x)\)。

发现对于 \(i\le k\) 的 \(a_i\)，代入后 \(D_0(x)G_0(x)=0\)，也就是说 \(F(a_i)=R_0(a_i)\)。

对于 \(i > k\) 的部分也可以同理构造多项式 \(G_1(x)\) 求出 \(R_1(x)\) 使得 \(F(a_i)=R_1(a_i)\)。

于是我们把用 \(O(n\log n)\) 的复杂度把原问题分成了两个规模减半的子问题。

每次的 \(G_0(x),G_1(x)\) 可以一开始都预处理出来。

预处理部分复杂度 \(O(n\log^2n)\)，总复杂度 \(O(n\log^2 n)\)。

注意到暴力展开计算常数非常小，在问题规模减到一定数值时可以暴力计算。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
struct tree_node{
    int ls, rs;
    poly g;
};
int _T_cnt;
std::vector<tree_node> _T;
int eval_inter_init(const poly &x, int l, int r){
    int u = _T_cnt++;
    _T[u].ls = _T[u].rs = -1, _T[u].g.resize(r - l + 1);
    if (l + 1 == r){
        _T[u].g[0] = minus(0, x[l]), _T[u].g[1] = 1;
        return u;
    }
    int md = (l + r + 1) >> 1;
    _T[u].ls = eval_inter_init(x, l, md);
    _T[u].rs = eval_inter_init(x, md, r);
    _T[u].g = Multiply(_T[_T[u].ls].g, _T[_T[u].rs].g);
    return u;
}
void Evaluate_solve(const poly &f, const poly &x, poly &y, int u, int l, int r){
    if (r - l <= 256){
        register int n = f.size();
        for (register int k = l; k < r; ++k){
            register unsigned now = 0;
            register unsigned long long b[17], c1, c2, c3, c4;
            b[0] = 1;
            for (register int i = 1; i <= 16; ++i) b[i] = b[i - 1] * x[k] % P;
            for (register int i = n - 1; i - 15 >= 0; i -= 16){
                c1 = now * b[16] + f[i] * b[15] + f[i - 1] * b[14] + f[i - 2] * b[13];
                c2 = f[i - 3] * b[12] + f[i - 4] * b[11] + f[i - 5] * b[10] + f[i - 6] * b[9];
                c3 = f[i - 7] * b[8] + f[i - 8] * b[7] + f[i - 9] * b[6] + f[i - 10] * b[5];
                c4 = f[i - 11] * b[4] + f[i - 12] * b[3] + f[i - 13] * b[2] + f[i - 14] * b[1];
                now = (c1 + c2 + c3 + c4 + f[i - 15]) % P;
            }
            for (register int i = n % 16 - 1; ~i; --i) now = (1ull * now * x[k] + f[i]) % P;
            y[k] = now;
        }
        return;
    }
    int md = (l + r + 1) >> 1;
    poly R;
    R = Divide(f, _T[_T[u].ls].g).second;
    Evaluate_solve(R, x, y, _T[u].ls, l, md);
    R = Divide(f, _T[_T[u].rs].g).second;
    Evaluate_solve(R, x, y, _T[u].rs, md, r);
}
poly Evaluate(const poly &f, const poly &x){
    int m = x.size(), rt = -1;
    _T_cnt = 0, _T.resize(2 * m - 1), rt = eval_inter_init(x, 0, m);
    poly res(m);
    Evaluate_solve(f, x, res, rt, 0, m);
    return res;
}

多项式快速插值

有一个 \(n-1\) 次多项式 \(F(x)\)，给定 \(n\) 个点及对应的点值 \((x_i,y_i)\)，求 \(F(x)\)。

根据拉格朗日插值公式，有 \[\begin{aligned}F(x)&=\sum_{i=1}^{n} y_i\prod_{i\ne j}\frac{x-x_j}{x_i-x_j}\\&=\sum_{i=1}^{n}\frac{y_i}{\prod_{i\ne j}(x_i-x_j)}\prod_{i\ne j}(x-x_j)\end{aligned}\]

令 \(G(x)=\prod_{i=1}^{n}(x-x_i),H_i(x)=\frac{G(x)}{x-x_i}\)，那么 \(\prod_{i\ne j}(x_i-x_j)=H_i(x_i)\)。根据洛必达法则，有 \(H_i(x_i)=G'(x_i)\)。

于是我们用多项式多点求值求出所有 \(G'(x_i)\) 即可。

记 \(a_i=\frac{y_i}{G'(x_i)}\)，那么有 \[F(x)=\sum_{i=1}^{n}a_i\prod_{i\ne j}(x-x_j)\]

然后我们分治求解。设当前区间为 \([l,r)\)，中点 \(m=\lfloor\frac{l+r+1}{2}\rfloor\)，我们要求 \[\begin{aligned}F_{l,r}(x)&=\sum_{i=l}^{r-1}a_i\prod_{l\le j < r,j\ne i}(x-x_j) \\ &=\left(\prod_{i=m}^{r-1}(x-x_i)\right)\sum_{i=l}^{m-1}a_i\prod_{l\le j < m,i\ne j}(x-x_j)+\left(\prod_{i=l}^{m-1}(x-x_i)\right)\sum_{i=m}^{r-1}a_i\prod_{m\le j < r,i\ne j}(x-x_j) \\ &=\left(\prod_{i=m}^{r-1}(x-x_i)\right)F_{l,m}(x)+\left(\prod_{i=l}^{m-1}(x-x_i)\right)F_{m,r}(x)\end{aligned}\]

时间复杂度 \(O(n\log^2 n)\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
poly Divide_2(poly a, int t){
    int n = a.size();
    poly b(n - 1);
    for (register int i = n - 1; i; --i){
        b[i - 1] = a[i];
        dec(a[i - 1], mul(a[i], t));
    }
    return b;
}
poly Inter_solve(const poly &x, const poly &y, int u, int l, int r){
    if (r - l <= 64){
        poly f(r - l), g = _T[u].g;
        for (register int i = l; i < r; ++i){
            poly tmp = Divide_2(g, minus(0, x[i]));
            for (register int j = 0; j < r - l; ++j) f[j] = (f[j] + 1ull * tmp[j] * y[i]) % P;
        }
        return f;
    }
    int md = (l + r + 1) >> 1;
    poly A = _T[_T[u].rs].g, B = Inter_solve(x, y, _T[u].ls, l, md);
    poly C = _T[_T[u].ls].g, D = Inter_solve(x, y, _T[u].rs, md, r);
    register int m = get(r - l);
    init(m), DFT(A, m), DFT(B, m), DFT(C, m), DFT(D, m);
    for (register int i = 0; i < m; ++i) A[i] = (1ull * A[i] * B[i] + 1ull * C[i] * D[i]) % P;
    IDFT(A, m), A.resize(r - l);
    return A;
}
poly Interpolation(poly x, poly y){
    int n = x.size(), rt = -1;
    _T_cnt = 0, _T.resize((n << 1) - 1), rt = eval_inter_init(x, 0, n);
    poly g = Derivative(_T[rt].g);
    poly res(n);
    Evaluate_solve(g, x, res, rt, 0, n);
    for (register int i = 0; i < n; ++i) y[i] = mul(y[i], qpow(res[i], P - 2));
    return Inter_solve(x, y, rt, 0, n);
}

多项式三角函数

不会

多项式反三角函数

不会

MTT

直接 FFT 的关键问题是精度不够。考虑将系数拆成 \(a\times 2^{15}+b\) 的形式，则两个系数相乘会得到 \[ \begin{aligned} & (a_1\times 2^{15}+b_1)(a_2\times 2^{15}+b_2) \\ =& a_1a_2\times 2^{30}+(a_1b_2+a_2b_1)\times 2^{15}+b_1b_2 \end{aligned} \]

做 \(4\) 次 DFT 以及 \(3\) 次 IDFT 即可。

这是 \(7\) 次 DFT 的做法，\(4\) 次甚至 \(3.5\) 次 DFT 的做法暂时还不会。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
struct Complex {
    double x, y;

    Complex operator + (const Complex &rhs) const {
        return {x + rhs.x, y + rhs.y};
    }

    Complex operator - (const Complex &rhs) const {
        return {x - rhs.x, y - rhs.y};
    }

    Complex operator * (const Complex & rhs) const {
        return {x * rhs.x - y * rhs.y, x * rhs.y + y * rhs.x};
    }
};

namespace MTTPoly {
    const int MAX_LEN = 1 << 18;
    const double pi = acos(-1);

    typedef std::vector<int> poly;

    Complex omega[MAX_LEN];

    void Init() {
        for (int m = 1; m < MAX_LEN; m <<= 1) {
            for (int i = 0; i < m; ++i) {
                omega[m + i] = {cos(pi / m * i), sin(pi / m * i)};
            }
        }
    }

    int get(int _n) {
        int n = 1;
        while (n < _n) {
            n <<= 1;
        }
        return n;
    }

    std::vector<int> rev;

    void init(int n) {
        int k = 0;
        while ((1 << k) < n) {
            ++k;
        }
        rev.resize(n), rev[0] = 0, --k;
        for (int i = 1; i < n; ++i) {
            rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << k);
        }
    }

    void DFT(std::vector<Complex> &f, int n) {
        f.resize(n);
        static Complex F[MAX_LEN];
        for (int i = 0; i < n; ++i) {
            F[rev[i]] = f[i];
        }
        for (int m = 1; m < n; m <<= 1) {
            for (int p = 0; p < n; p += m << 1) {
                Complex *W = omega + m;
                Complex *F0 = F + p, *F1 = F0 + m;
                for (int i = 0; i < m; ++i, ++W, ++F0, ++F1) {
                    Complex t = (*F1) * (*W);
                    *F1 = *F0 - t, *F0 = *F0 + t;
                }
            }
        }
        for (int i = 0; i < n; ++i) {
            f[i] = F[i];
        }
    }

    void IDFT(std::vector<Complex> &f, int n) {
        DFT(f, n), std::reverse(f.begin() + 1, f.end());
        for (int i = 0; i < n; ++i) {
            f[i] = f[i] * (Complex){1.0 / n, 0};
        }
    }

    poly Mul(const poly &f, const poly &g, int P) {
        int fn = f.size(), gn = g.size(), _n = fn + gn - 1;
        std::vector<Complex> fa(fn), fb(fn), ga(gn), gb(gn);
        for (int i = 0; i < fn; ++i) {
            fa[i].x = f[i] >> 15, fa[i].y = 0;
            fb[i].x = f[i] & 32767, fb[i].y = 0;
        }
        for (int i = 0; i < gn; ++i) {
            ga[i].x = g[i] >> 15, ga[i].y = 0;
            gb[i].x = g[i] & 32767, gb[i].y = 0;
        }
        int n = get(_n);
        init(n), DFT(fa, n), DFT(fb, n), DFT(ga, n), DFT(gb, n);
        std::vector<Complex> a(n), b(n), c(n);
        for (int i = 0; i < n; ++i) {
            a[i] = fa[i] * ga[i];
            b[i] = fa[i] * gb[i] + fb[i] * ga[i];
            c[i] = fb[i] * gb[i];
        }
        IDFT(a, n), IDFT(b, n), IDFT(c, n);
        poly ans(_n);
        for (int i = 0; i < _n; ++i) {
            long long A = a[i].x + 0.5, B = b[i].x + 0.5, C = c[i].x + 0.5;
            A %= P, B %= P, C %= P;
            ans[i] = ((1 << 30) * A + (1 << 15) * B + C) % P;
            ans[i] = (P + ans[i]) % P;
        }
        return ans;
    }
}

完整代码

完整代码可以见 Codes/Polynomial.cpp at master · AutumnKite/Codes。

博客内容遵循 署名-非商业性使用-相同方式共享 4.0 国际 (CC BY-NC-SA 4.0) 协议

本文永久链接是：https://autumnkite.github.io/polynomial-notes/

更新于：2020年6月14日

FFT

多项式

评论