多项式

首先，由任意n+1个不同的点，可以求出一个唯一的n次多项式（当然可以高斯消元，但事实上一般都用拉格朗日插值

拉格朗日插值

这不就来了

我们知道了n+1个不同的函数上的点，就得过来求一下函数的值对吧

拉插主要是对于一个多项式知道了n+1个点之后，求一个定点的函数值的方法

我们发现：

这个东西我不会证明

背式子：

$$\huge{f(k) = \sum_{i=1}^{n}y_i\prod_{j!=i}^{n}\dfrac{k-x_j}{x_i-x_j}}$$

这个东西在算的时候就可以在算后边式子的时候把下边的乘积和上面的乘积求出来，然后下面整体求逆元就好了

横坐标连续的拉格朗日插值

为什么把这玩意单独拿出来捏？

因为这个东西可以直接O(nlogn)求

我们就假设是从1-n了横坐标

$$\huge{f(k) = \sum_{i=1}^{n}y_i\prod_{j!=i}^{n}\dfrac{k-j}{i-j}}$$

我们发现后边连乘的东西可以O（1）算出来，我们预处理出来一个num表示$\prod_{i=1}^{n}k-i$,那么后边的分子就变成了$num*inv(k-i)$，下面的东西也好说，我们发现在j<i的时候下面的乘积等于$(i-1)!$，在$j>i$的时候，它变成了$(n-i)!*(-1)^{n-i}$直接算就行了

FFT

前置

FFT的算法流程有两个：

DFT and IDFT

DFT是将多项式的系数表达转化成点值表达的算法

IDFT是反过来

我们发现知道一个x很容易O(n)知道它的点值，然后确定一个多项式需要n+1个点值

然后你成功$n^2$了，还不如直接暴力

但是我们的傅里叶先生既然提出了要DFT那就肯定有人家的道理对吧

然后我们发现，n个普通的点如果求点值确实要这么干，比较麻烦

但是n个特殊点呢……

欸，n个特殊点，什么点特殊呢？连续的正整数？（没用），特殊的有理数？（没用）

那用啥啊？

用复数，准确来说是单位根$\omega_{n+1}^0 - \omega_{n+1}^{n}$，他们有一些性质比较特殊，然后可以分治，然后就$n\log n$了

复数大家都会了我不说了直接过

来我们把单位根的性质给出来：

1. 对于任意的n，$\omega_n^0 = 1$
2. $\omega_n^k = (\omega_n^1)^k$
3. $\omega_n^j*\omega_n^k = \omega_n^{j+k}$
4. $\omega_{2n}^{2k} = \omega_n^k$
5. 如果n是偶数，$\omega_n^{k+n/2} = -\omega_n^k$

然后我们就开始下一步！

DFT的加速版本

我们来讲一下FFT怎么分治嗷

我们先假设n是2的整次幂

对于一个多项式$F(x)$，我们把它的偶数项和奇数项分别拿出来变成$FL(x^2),FR(x^2)$

然后我们发现$FL(x^2) = F[0]+F[2]x^2+...+F[n-2]x^{n-2},FR(x^2) = F[1]x+F[3]x^3+...+F[n-1]x^{n-1}$

然后发现：$FL(x^2) = F[0]+F[2]x^2+...+F[n-2]x^{n-2},FR(x^2) = x(F[1]+F[3]x^2+...+F[n-1]x^{n-2})$

那我们让$FR$除掉x，使得$F(x) = FL(x^2)+xFR(x^2)$

然后我们把$\omega_n^k$代入

• $k<n/2$代入的时候

  $F(\omega_n^k) = FL((\omega_n^k)^2)+w_n^kFR((\omega_n^k)^2)$

由于$\omega_n^k = \omega_{n/2}^k$

所以：$F(\omega_n^k) = FL(\omega_{n/2}^k)+w_n^kFR(\omega_{n/2}^k)$

• 那我们再带入$\omega_n^{k+n/2}$的时候会变成$F(\omega_n^{k+n/2}) = FL(\omega_{n/2}^k)-w_n^kFR(\omega_{n/2}^k)$

很神奇地发现，这俩只差这一个正负号

那这个有什么用嘞？

我们只要算出来$FL,FR$在$\omega_{n/2}^k$的值就可以算出$F$在$\omega_n^k~ and~ \omega_n^{k+n/2}$了

然后我们考虑分治一下，先将$FL,FR$在这些位置的点值求出来再求$F$的点值

这样就可以$n\log n$了

还有就是刚开始说n是2的整次幂的时候，如果不是呢？

（那补成2的整次幂就行了啊，后边全是0

IDFT

结论：把DFT中的$\omega_n^1$换成$\omega_n^{-1}$就行了，然后再除以n

证明很有价值但是有点难，要用到单位根反演

然后我们就可以写出来FFT的初始版本了

首先我们发现IDFT和DFT区别不大，所以我们直接写进一个函数里面，传进去一个bool函数表示是哪个运算就用哪个单位根就行了

然后呢，我们发现数组拷贝比较慢，所以我们发现把，每一次左右分，分完之后i的最后位置会是它二进制反转之后的位置（这个很好证

那么我们先反转之后向上合并就能使常数变小

怎么求出来二进制翻转后在哪呢？

$tr[i]=(tr[i>>1]>>1)|((i\&1)?limit>>1:0);$

前半部分是除掉最后一位反转之后的数字，后半部分是最后一位反转之后是不是1，由于limit>>1之后就是最高位是1所以直接或上

这个优化叫做蝴蝶变换（很优雅的名字

然后就写出来了：

// Code by johnsonloy

#include <bits/stdc++.h>
using namespace std;

#define db double
const db pi = acos(-1);
const int maxn = 4e6 + 10;

int n, m, tr[maxn];
db ans[maxn];

struct node {
    db x = 0, y = 0;
    node operator+(const node a) {
        return {x + a.x, y + a.y};
    }
    node operator-(const node a) {
        return {x - a.x, y - a.y};
    }
    node operator*(const node a) {
        return {x * a.x - y * a.y, y * a.x + x * a.y};
    }
} f[maxn], g[maxn];

void fft(node* f, bool flag) {
    for (int i = 0; i < n; i++)
        if (i < tr[i])
            swap(f[i], f[tr[i]]);
    //首先是蝴蝶变换
    //蝴蝶变换的意思差不多就是最后的序列x位置的数是刚开始的序列的x把二进制反转之后得到的y位置的数字
    for (int p = 2; p <= n; p <<= 1) {
        // p是当前要合并的区间的长度
        int len = p >> 1;
        //左右区间的长度
        node tg = {cos(2 * pi / p), sin(2 * pi / p)};
        //先算出来p次方的复数单位根
        if (flag == 0)
            tg.y *= -1;
        //如果现在是IDFT的话
        for (int k = 0; k < n; k += p) {
            //枚举一下所有的长度是这么多的起始点
            node buf = {1, 0};
            // 0次单位根
            for (int l = k; l < k + len; l++) {
                //每一次进行区间合并的时候的左区间，由于每合并了之后都要让单位根的上指标加1所以在下面有个操作是buf*tg
                node tt = buf * f[len + l];
                f[len + l] = f[l] - tt;
                f[l] = f[l] + tt;
                //这就是推出来的式子了
                // F(w(n,k)) = FL(w(n/2,k))+w(n,k)FR(w(n/2,k))
                // F(w(n,k+n/2)) = FL(w(n/2,k))-w(n,k)FR(w(n/2.k))
                buf = buf * tg;
            }
        }
    }
}

signed main() {
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
    freopen("out.txt", "w", stdout);
#endif
    ios::sync_with_stdio(false);
    cin >> n >> m;
    for (int i = 0; i <= n; i++)
        cin >> f[i].x;
    for (int i = 0; i <= m; i++)
        cin >> g[i].x;
    for (m += n, n = 1; n <= m; n <<= 1)
        ;
    for (int i = 0; i < n; i++)
        tr[i] = (tr[i >> 1] >> 1) | ((i & 1) ? n >> 1 : 0);
    fft(f, 1), fft(g, 1);
    //先把两个函数都转化成点值的形式
    for (int i = 0; i < n; i++)
        f[i] = f[i] * g[i];
    fft(f, 0);
    //然后进行IDFT的合并，从点值转化成系数
    for (int i = 0; i <= m; i++)
        cout << (int)(f[i].x / n + 0.49) << ' ';
    //这个操作……简称四舍五入
    cout << endl;
    return 0;
}

NTT与多项式全家桶

很多时候我们的多项式计数都要有模数，FFT的单位根需要预处理什么的会丢精度

那么我们想找到一个模意义下单位根的替代品，谁呢？原根咯

我们发现单位根的所有性质原根g在模意义下不一定都具备

但是，当有n的时候$g^{(mod-1)/n}$是满足条件的，那么就用它来替代单位根

然后改一下代码就出来了NTT

#include <bits/stdc++.h>
using namespace std;

#define ll long long
#define pii pair<int, int>
#define mp make_pair
const int maxn = 4e6 + 10;
const int mod = 998244353, G = 3;
const ll inf = 0x3f3f3f3f3f3f3f3f;

namespace IO {
void openfile() {
#ifndef ONLINE_JUDGE
    freopen("in.in", "r", stdin);
    freopen("out.out", "w", stdout);
#endif
}

int mul(int x, int y) {
    return 1ll * x * y % mod;
}

int add(ll x, ll y) {
    x += y;
    if (x >= mod)
        x -= mod;
    return x;
}

int sub(int x, int y) {
    return add(x, mod - y);
}

inline int read() {
    int x = 0, f = 0;
    char c = getchar();
    while (!isdigit(c))
        f |= c == '-', c = getchar();
    while (isdigit(c))
        x = x * 10 + c - '0', c = getchar();
    if (f)
        x = -x;
    return x;
}
}  // namespace IO
using namespace IO;

int f[maxn], g[maxn];
int n, m;
int tr[maxn];

int qpow(int x, int y) {
    int ans = 1;
    while (y) {
        if (y & 1)
            ans = mul(ans, x);
        x = mul(x, x);
        y >>= 1;
    }
    return ans;
}

const int invG = qpow(G, mod - 2);

void ntt(int* f, int op) {
    for (int i = 0; i < n; i++)
        if (i < tr[i])
            swap(f[i], f[tr[i]]);
    for (int p = 2; p <= n; p <<= 1) {
        int len = p >> 1;
        int tg = qpow(op ? G : invG, (mod - 1) / p);
        for (int l = 0; l < n; l += p) {
            int buf = 1;
            for (int k = l; k < l + len; k++) {
                int tmp = mul(buf, f[k + len]); 
                f[k + len] = sub(f[k], tmp);
                f[k] = add(f[k], tmp);
                buf = mul(buf, tg);
            }
        }
    }
    if (!op) {
        int invn = qpow(n, mod - 2);
        for (int i = 0; i < n; i++)
            f[i] = mul(f[i], invn);
    }
}

signed main() {
    openfile();
    n = read(), m = read();
    for (int i = 0; i <= n; i++)
        f[i] = read();
    for (int i = 0; i <= m; i++)
        g[i] = read();
    for (m = n + m, n = 1; n <= m; n <<= 1)
        ;
    for (int i = 0; i < n; i++)
        tr[i] = (tr[i >> 1] >> 1) | ((i & 1) ? n >> 1 : 0);
    ntt(f, 1), ntt(g, 1);
    for (int i = 0; i < n; i++)
        f[i] = mul(f[i], g[i]);
    ntt(f, 0);
    for (int i = 0; i <= m; i++)
        printf("%lld ", f[i]);
    cerr << 1.0 * clock() / CLOCKS_PER_SEC << '\n';
    return 0;
}

然后开始全家桶生涯：

先说多项式牛顿迭代吧

多项式求导：

$$F'(x) = \sum_{i = 0} iF[i]x^{i-1}$$

多项式积分：

$$F'(x) = C+\sum_{i=0}\dfrac{F[i]x^{i+1}}{i}$$

定义多项式复合：$F(G(x)) = \sum_{i=0}F[i]G(x)^i$

• 多项式牛顿迭代：

用于解决：已知函数$G$并且$G(F(x)) = 0$求多项式$F(\mod x^n)$

结论：$F(x) = F^*(x) - \dfrac{G(F^*(x))}{G'(F^*(x))}(\mod x^n)$

定义$F^*(x)$代表函数$F$在$\mod x^{n/2}$意义下的函数

一定要好好记住

• 比如多项式求逆：

  根据$B(x)*A(x) - 1 = 0(\mod x^n)$，$A(x)$已知。

  设$G(x)$且$G(B(x)) = A(x)B(x) - 1$这里$A(x)$为系数

  所以$G'(B(x)) = A(x)(\mod x^n)$

  得：$B(x) = B^*(x) - \dfrac{G(B^*(x))}{G'(B^*(x))} = B^*(x) - \dfrac{1}{A(x)}(A(x)B^*(x) - 1)$

  由于$\dfrac{1}{A(x)}$的精度只需要达到$\mod x^{n/2}$所以这玩意就等于$B^*(x)$

  所以：$B(x) = B^*(x)(2-A(x)B^*(x))$

• 还有多项式开根：

  此处$A(x)$已知

  令：$G(B(x)) = B(x)^2-A(x)$

  则$G'(B(x)) = 2B(x)$

  根据牛顿迭代可以得到：$B(x) = B^*(x) - \dfrac{G(B^*(x))}{G'(B^*(x))} = B^*(x) - \dfrac{B^*(x)^2 - A(x)}{2B^*(x)} = \dfrac{A(x)+B^*(x)^2}{2B^*(x)}$

  多项式只有常数的时候这玩意就是二次剩余

多项式带余除法

结论：给定$F,G$用$F$的反转和$G$的反转的逆乘起来模$x^{n+m-1}$就是商，然后余数直接算，别问我为啥有点难

多项式ln,exp

$ln'(x) = \dfrac{1}{x}$

题目有：$ln(A(x)) = B(x)$,两边同时求导：

$\dfrac{d}{dx}ln(a9x)) = \dfrac{d}{dx}B(x)$

$\dfrac{dA(x)}{dx}\dfrac{d}{dA(x)}lnA(x) = B'(x)$

$\dfrac{A'(x)}{A(x)} = B'(x)$

同时积分得到：$B(x) = \int(A'(x)/A(x))$

exp和ln是逆运算，考虑牛顿迭代：

$e^(A(x)) = B(x)$设，$G(B(x)) = lnB(x) - A(x)$

求导得到：$G'(B(x)) = B(x)^{-1}$

$B(x) = B^*(x) - \dfrac{G(B^*(x))}{G'(B^*(x))} = B^*(x) - \dfrac{lnB^*(x) - A(x)}{B*(x)^{-1}}$

$B(x) = B^*(x) - (lnB^*(x) - A(x))B^*(x)$

$B(x) = (1-lnB^*(x)+A(x))B^*(x)$

保证$A[0] = 0$，此时$B[0] = 1$否则会出现无穷求和

多项式快速幂

公式很简单$A^k(x) = e^{ln(A(x))*k}$

但是要保证$A[0] = 1$

如果不保证直接把前面是0的拿出来并且把整个多项式除掉$A[p]$就好了,$p$是第一个非空的$A[p]$

#include <bits/stdc++.h>
using namespace std;

#define cpy(f, g, n) memcpy(f, g, sizeof(int) * (n))
#define clr(f, n) memset(f, 0, sizeof(int) * (n))
#define int long long
#define pii pair<int, int>
#define mp make_pair
const int maxn = 1e6 + 10;
const int mod = 998244353;
const int inf = 0x3f3f3f3f3f3f3f3f;

namespace IO {
void openfile() {
#ifndef ONLINE_JUDGE
    freopen("in.in", "r", stdin);
    freopen("out.out", "w", stdout);
#endif
}

int mul(int x, int y) {
    return 1ll * x * y % mod;
}

int add(int x, int y) {
    x += y;
    if (x >= mod)
        x -= mod;
    return x;
}

int sub(int x, int y) {
    return add(x, mod - y);
}

inline int read() {
    int x = 0, f = 0;
    char c = getchar();
    while (!isdigit(c))
        f |= c == '-', c = getchar();
    while (isdigit(c))
        x = x * 10 + c - '0', c = getchar();
    if (f)
        x = -x;
    return x;
}
}  // namespace IO
using namespace IO;

int qpow(int x, int y = mod - 2) {
    int ans = 1;
    while (y) {
        if (y & 1)
            ans = mul(ans, x);
        x = mul(x, x);
        y >>= 1;
    }
    return ans;
}

namespace duo {
const int _G = 3, invG = qpow(_G);
int n, F[maxn], k, k2, inv[maxn];
int m, G[maxn], tr[maxn];

void print(int* f, int n) {
    for (int i = 0; i < n; ++i)
        printf("%lld ", f[i]);
    puts("");
}

void px(int* f, int* g, int n) {
    for (int i = 0; i < n; i++)
        f[i] = mul(f[i], g[i]);
}

void init_tr(int n) {
    for (int i = 0; i < n; i++)
        tr[i] = (tr[i >> 1] >> 1) | (i & 1 ? n >> 1 : 0);
}

void rev(int* f, int n) {
    for (int i = 0; i < n; i++)
        if (i < n - 1 - i)
            swap(f[i], f[n - 1 - i]);
}

void ntt(int* f, bool op, int n) {
    init_tr(n);
    for (int i = 0; i < n; i++)
        if (i < tr[i])
            swap(f[i], f[tr[i]]);
    for (int p = 2; p <= n; p <<= 1) {
        int len = p >> 1;
        int tg = qpow(op ? _G : invG, (mod - 1) / p);
        for (int l = 0; l < n; l += p) {
            int buf = 1;
            for (int k = l; k < l + len; k++) {
                int tmp = tg * f[k + len];
                f[k + len] = sub(f[k], tmp);
                f[k] = add(f[k], tmp);
                buf = mul(buf, tg);
            }
        }
    }
    if (!op)
        for (int invn = qpow(n), i = 0; i < n; i++)
            f[i] = mul(f[i], invn);
}

void times(int* f, int* g, int n1, int lim) {
    static int s[maxn];
    int n = 1;
    for (; n < n1 + n1; n <<= 1)
        ;
    cpy(s, g, n);
    ntt(f, 1, n), ntt(s, 1, n);
    px(f, s, n);
    ntt(f, 0, n);
    clr(s, n);
    clr(f + lim, n - lim);
}

void jia(int* f, int* g, int n) {
    for (int i = 0; i < n; i++)
        f[i] = add(f[i], g[i]);
}

void jian(int* f, int* g, int n) {
    for (int i = 0; i < n; i++)
        f[i] = sub(f[i], g[i]);
}

void invp(int* f, int m) {
    int n = 1;
    for (; n < m; n <<= 1)
        ;
    static int w[maxn], s[maxn], r[maxn];
    w[0] = qpow(f[0]);
    for (int len = 2; len <= n; len <<= 1) {
        for (int i = 0; i < (len >> 1); i++)
            r[i] = mul(w[i], 2);
        cpy(s, f, len);
        ntt(w, 1, len << 1);
        px(w, w, len << 1);
        ntt(s, 1, len << 1);
        px(w, s, len << 1);
        ntt(w, 0, len << 1);
        clr(w + len, len);
        for (int i = 0; i < len; i++)
            w[i] = sub(r[i], w[i]);
    }
    cpy(f, w, m);
    clr(w, n + n), clr(s, n + n), clr(r, n + n);
}
// R(x) = 2R*(x) - R*(x)^2F(x)

void dao(int* f, int n) {
    for (int i = 0; i < n - 1; i++)
        f[i] = mul(f[i + 1], i + 1);
    f[n - 1] = 0;
}

void init(int n) {
    inv[1] = 1;
    for (int i = 2; i <= n; i++)
        inv[i] = mul(inv[mod % i], mod - mod / i);
}

void jifen(int* f, int n) {
    for (int i = n; i; i--)
        f[i] = mul(f[i - 1], inv[i]);
    f[0] = 0;
}
//牛顿迭代：G(F(x)) = 0 ,F(x) = F*(x)-G(F*(x))/G'(F*(x)) (mod x^n)

void sqrtn(int* f, int m) {
    int n = 1;
    for (; n < m; n <<= 1)
        ;
    static int b[maxn], b2[maxn];
    b[0] = 1;
    for (int len = 2; len <= n; len <<= 1) {
        for (int i = 0; i < (len >> 1); i++)
            b2[i] = mul(b[i], 2);
        invp(b2, len);
        ntt(b, 1, len), px(b, b, len);
        ntt(b, 0, len);
        jia(b, f, len);
        times(b, b2, len, len);
    }
    cpy(f, b, m);
    clr(b, n + n), clr(b2, n + n);
}
// B(x) = (A(x)+B*(x)^2)/2B*(x)

void mof(int* f, int* g, int n, int m) {
    static int q[maxn], t[maxn];
    int l = n - m + 1;
    rev(g, m), cpy(q, g, l), rev(g, m);
    rev(f, n), cpy(t, f, l), rev(f, n);
    invp(q, l), times(q, t, l, l);
    times(g, q, n, n);
    for (int i = 0; i < m; i++)
        g[i] = sub(f[i], g[i]);
    clr(g + m - 1, l);
    cpy(f, q, l);
    clr(f + l, n - l);
}
// QT(x) = FT(x)*GT(x)^{-1}

void lnp(int* f, int n) {
    static int g[maxn];
    cpy(g, f, n);
    dao(f, n), invp(g, n);
    times(f, g, n, n);
    jifen(f, n - 1);
    clr(g, n);
}
// B(x) = ji(A'(x)/A(x))

void exp(int* f, int m) {
    static int b[maxn], b2[maxn];
    int n = 1;
    for (; n < m; n <<= 1)
        ;
    b2[0] = 1;
    for (int len = 2; len <= n; len <<= 1) {
        cpy(b, b2, len >> 1);
        lnp(b, n);
        for (int i = 0; i < len; i++)
            b[i] = sub(f[i], b[i]);
        b[0] = add(b[0], 1);
        times(b2, b, len, len);
    }
    cpy(f, b2, m), clr(b, n), clr(b2, n);
}
// B(x) = (1-lnB*(x)+A(x))B*(x)

void dqpow(int* f, int n, int k) {
    lnp(f, n);
    for (int i = 0; i < n; i++)
        f[i] = mul(f[i], k);
    exp(f, n);
}
// A(x)^k = e^{ln(A(x))*k}

}  // namespace duo

signed main() {
    openfile();
    cerr << 1.0 * clock() / CLOCKS_PER_SEC << '\n';
    return 0;
}

然后就给两道题：

P3338 [ZJOI2014]力

我们要求：$E_j = \sum_{i=1}^{j-1}\dfrac{q_i}{(i-j)^2}-\sum_{i=j+1}^{n}\dfrac{q_i}{(i-j)^2}$

然后我们来转化一下：

设$f[i] = q_i,g[i] = \dfrac{1}{i^2}$原式变成：

$\sum_{i=1}^jf[i]*g[j-i]-\sum_{i=j}^{n}f[i]*g[i-j]$

令$f[0] = g[0] = 0$变成：

$\sum_{i=0}^jf[i]*g[j-i]-\sum_{i=j}^{n}f[i]*g[i-j]$

左边已经是卷积的形式了，右边还不是，但是我们设$f'[i] = f[n-i]$

就变成了：$\sum_{i=1}^jf[i]*g[j-i]-\sum_{i=j}^{n}f'[n-i]*g[i-j]$

$\sum_{i=1}^jf[i]*g[j-i]-\sum_{i=0}^{n-j}f'[n-j-i]*g[i]$

然后就都是卷积了

然后我们设多项式$A(x) = \sum_{i=0}^n f[i],B(x) = sum_{i=0}^n g[i],C(x) = \sum_{i=0}^n f'[i]$

令$L(x) = A(x)*B(x),R(x) = C(x)*B(x)$

所以$ans[i] = L[i]-R[n-i]$

提醒：会被卡精度，所以我们应该将处理$g_i$变成：$1/i/i$这样精度丢的少

CF528D Fuzzy Search

我们考虑一个常见套路嗷，用FFT求字符串匹配：

$(s_i-s_j)^2*s_i*s_j$，通配符的权值为0

然后把模式串倒过来就成了：$(s_i-s_{m-j})^2*s_i*s_{m-j}$

然后就可以卷积了，最后的m+i项如果是0说明全都匹配，否则就不是

这个题可以稍微简单一点，我们对每一个字母考虑，比如本来是AGCAATTCAT和ACAT

对A考虑就是：A##AA###A#和A%A%匹配，$cnt_i$表示第i位开头匹配了多少位

然后使得所有A位置都是1，其他都是0

然后发现没有k的处理，我们发现只要把上面的文本串每一个A的左右两边变成A就行了

然后直接匹配看匹配了多少位，如果匹配了m位就说明可行