FFT与拉格朗日插值

2022年7月7日和8日，day04

背景

• 计数问题利用生成函数的方法转化为形式幂级数的代数计算问题。
• 假设 $C(x)=A(x)+B(x)$ 或 $A(x)B(x)$ 或 $\frac{A(x)}{B(x)}$ 或 $\sqrt{A(x)}$ 等形式，$C(x)$ 的前 $n+1$ 项的系数 $c_0,c_1,\dots ,c_n$ 仅仅需要 $A(x),B(x)$ 的前 $n+1$ 项系数就可以计算。
• 多项式 $A(x),B(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n+1}$ 代替形式幂级数 $A(x),B(x)$。
• 朴素多项式乘法、求逆、开根号等的复杂度为 $O(n^2)$。
• $FFT/NTT$ 可以 $O(nlog(n))$ 内实现两个 $n$ 次多项式的乘积。

多项式的表示形式

系数表示与点值表示

假设 $f(x)$ 是一个 $n$ 次多项式，则 $f(x)$ 的系数表示为 $f(x)=a_n{x}^n+a_{n-1}{x}^{n-1}+\cdots +a_0$。$f(x)$ 的点值表示为 $(x_0,f(x_0)),(x_1,f(x_1)),\dots ,(x_n,f(x_n))$。

• $n+1$ 个点值可以表示一个 $n$ 次多项式。
• 在点值表示下 $n$ 次多项式的乘法复杂度为 $O(n)$。
  $$\begin{gathered} h(x)=f(x)\times g(x),h(x)是n次多项式 \\ \begin{array}{rl}& (x_0,f(x_0)),(x_1,f(x_1)),\dots ,(x_n,f(x_n))\\ & (x_0,g(x_0)),(x_1,g(x_1)),\dots ,(x_n,g(x_n))\\ \Rightarrow & (x_0,h(x_0)),(x_1,h(x_1)),\dots ,(x_n,h(x_n))\end{array} \\ h(x_i)=f(x_i)\times g(x_i) \end{gathered}$$

复数与单位根

复数的指数形式：$a+bi=r{e}^{i\theta }$，其中 $r=\sqrt{a^2+b^2},tan(\theta )=\frac{b}{a}$。

单位根：$x^n=1$ 在复数域上的根称为 $n$ 次单位根。$n$ 次单位根有 $n$ 个，形式为 $w_n^k=e^{i\frac{2k\pi }{n}}$。

欧拉公式：$e^{ix}=cosx+isinx$。
$$w_n=e^{i\frac{2\pi }{n}}\Rightarrow (w_n{)}^k=e^{i\frac{2k\pi }{n}}\Rightarrow (w_n^k{)}^n=e^{i\frac{2kn\pi }{n}}=1$$
$w_n^k$ 可以理解为 $x^n=1$ 在复数域上的第 $k+1$ 个根。
单位根的性质：

• $w_n^k=w_{2n}^{2k}$
• $w_{2n}^{k+n}=-w_{2n}^k$

FFT(快速傅里叶变换，Fast Fourier Transform)

DFT(离散傅里叶变换，discrete Fourier Transform)

将多项式 $A(x)=a_0+a_{1}x+\cdots +a_{n-1}{x}^{n-1}$ 转化为其点值表达式 $(w_n^k,A(w_n^k))$，$k=0,1,\dots ,n-1$。

不妨设 $n=2^l$
$$\begin{array}{rl}A(x)=& a_0+a_{1}x+\cdots +a_{n-1}{x}^{n-1}\\ A(x)=& (a_0+a_2{x}^2+\cdots +a_{n-2}{x}^{n-2})\\ & +(a_{1}x+a_3{x}^3+\cdots +a_{n-1}{x}^{n-1})\end{array}$$
令 $B(x),C(x)$ 为
$$\begin{gathered} B(x)=a_0+a_{2}x+\cdots +a_{n-2}{x}^{\frac{n}{2}-1} \\ C(x)=a_1+a_{3}x+\cdots +a_{n-1}{x}^{\frac{n}{2}-1} \end{gathered}$$
那么 $A(x)$ 重写为：
$$A(x)=B(x^2)+xC(x^2)$$
分两种情况：

• $0\le k<\frac{n}{2}$
  $$A(w_n^k)=B(w_n^{2k})+w_n^{k}C(w_n^{2k})=B(w_{\frac{n}{2}}^k)+w_n^{k}C(w_{\frac{n}{2}}^k)$$
• $\frac{n}{2}\le k+\frac{n}{2}<n$
  $$A(w_n^{k+\frac{n}{2}})=B(w_n^{2k})+w_n^{k+\frac{n}{2}}C(w_n^{2k})=B(w_{\frac{n}{2}}^k)-w_n^{k}C(w_{\frac{n}{2}}^k)$$
  时间复杂度：$T(n)=2T(\frac{n}{2})+O(n)\Rightarrow T(n)=O(nlogn)$。
  IDFT(逆离散傅里叶变换，inverse DFT)

将多项式的点值表示 $(w_n^k,b_k)$，$(k=0,1,\dots ,n-1)$ 转化为其系数表示 $A(x)=a_0+a_{1}x+\cdots +a_{n-1}{x}^{n-1}$。
$$\left[\begin{array}{cccc}{w}_n^{0\cdot 0}& w_n^{0\cdot 1}& \cdots & w_n^{0\cdot (n-1)}\\ w_n^{1\cdot 0}& w_n^{1\cdot 1}& \cdots & w_n^{1\cdot (n-1)}\\ ⋮& ⋮& & ⋮\\ w_n^{(n-1)\cdot 0}& w_n^{(n-1)\cdot 1}& \cdots & w_n^{(n-1)\cdot (n-1)}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}{b}_0\\ b_1\\ ⋮\\ b_{n-1}\end{array}\right]$$

设 $n\times n$ 的矩阵 $\Omega$，其中 ${\Omega }_{i,j}=w_n^{ij}$，$(0\le i,j<n)$，设向量 $\mathbf{a}=(a_0,a_1,\cdots ,a_{n-1})$，$\mathbf{b}=(b_0,b_1,\cdots ,b_{n-1})$，则 IDFT 相当于求解方程
$$\Omega \mathbf{a}=\mathbf{b}$$
记 ${\bar{\Omega }}_{i,j}=w_n^{-ij}$，那么 $M=\bar{\Omega }\cdot \Omega$ 为
$$\left[\begin{array}{cccc}{w}_n^{-0\cdot 0}& w_n^{-0\cdot 1}& \cdots & w_n^{-0\cdot (n-1)}\\ w_n^{-1\cdot 0}& w_n^{-1\cdot 1}& \cdots & w_n^{-1\cdot (n-1)}\\ ⋮& ⋮& & ⋮\\ w_n^{-(n-1)\cdot 0}& w_n^{-(n-1)\cdot 1}& \cdots & w_n^{-(n-1)\cdot (n-1)}\end{array}\right]\left[\begin{array}{cccc}{w}_n^{0\cdot 0}& w_n^{0\cdot 1}& \cdots & w_n^{0\cdot (n-1)}\\ w_n^{1\cdot 0}& w_n^{1\cdot 1}& \cdots & w_n^{1\cdot (n-1)}\\ ⋮& ⋮& & ⋮\\ w_n^{(n-1)\cdot 0}& w_n^{(n-1)\cdot 1}& \cdots & w_n^{(n-1)\cdot (n-1)}\end{array}\right]=\left[\begin{array}{cccc}n& 0& \cdots & 0\\ 0& n& \cdots & 0\\ ⋮& ⋮& & ⋮\\ 0& 0& \cdots & n\end{array}\right]$$
由于
$$\begin{array}{rl}{M}_{ij}& =w_n^{-i\cdot 0}{w}_n^{0\cdot j}+w_n^{-i\cdot 1}{w}_n^{1\cdot j}+\cdots +w_n^{-i\cdot (n-1)}{w}_n^{(n-1)\cdot j}\\ & =w_n^{(j-i)\cdot 0}+w_n^{(j-i)\cdot 1}+\cdots +w_n^{(j-i)\cdot (n-1)}\\ & =\{\begin{array}{llc}n& \text{if}& j=i\\ \frac{1-(w_n^n{)}^{(j-i)}}{1-w_n^{j-i}}=0& \text{if}& j\ne i\end{array}\end{array}$$
因此有 $\bar{\Omega }\Omega =nI$，那么
$$\mathbf{a}=\frac{1}{n}\bar{\Omega }\mathbf{b}$$
相当于给定 $B(x)=b_0+b_{1}x+\cdots +b_{n-1}{x}^{n-1}$，求点值 $B(w_n^{-k})$，$(0\le k<n)$，那么
$$a_k=\frac{1}{n}B(w_n^{-k})$$

例题1：https://ac.nowcoder.com/acm/contest/26012/A

记 $A$ 和 $B$ 分别为
$$\begin{gathered} A=a_n\cdot 10^n+\cdots +a_i\cdot 10^i+\cdots +a_1\cdot 10+a_0 \\ \Rightarrow A(x)=a_n\cdot x^n+\cdots +a_i\cdot x^i+\cdots +a_1\cdot x+a_0 \\ B=b_m\cdot 10^m+\cdots +b_j\cdot 10^j+\cdots +b_1\cdot 10+b_0 \\ \Rightarrow B(x)=b_m\cdot x^m+\cdots +b_j\cdot x^j+\cdots +b_1\cdot x+b_0 \end{gathered}$$
那么 $C$ 为
$$\begin{gathered} C(x)=A(x)\times B(x) \\ \Rightarrow C(x)=c_{n+m}\cdot x^{n+m}+\cdots +c_i\cdot x^i+\cdots +c_1\cdot x+c_0 \end{gathered}$$

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

const double Pi = acos(-1.0);
const int N = 2e5 + 50;

struct com {
    double x, y;
    com operator + (com r) { return {x + r.x, y + r.y}; }
    com operator - (com r) { return {x - r.x, y - r.y}; }
    com operator * (com r) { return {x * r.x - y * r.y, x * r.y + y * r.x}; }
}f[N], g[N];

int n, m, ans[N];
int rev[N], lim, bit;

int read (com* c) {
    static char s[N]; cin >> s;
    int len = strlen (s);
    for (int i = len - 1; i >= 0; i --) 
        c[i].x = s[len - i - 1] - '0';
    return len;
}

void fft (com* h, int inv) {
    for (int i = 0; i < lim; i ++) 
        if (rev[i] < i) swap (h[i], h[rev[i]]);
    for (int k = 1; k < lim; k <<= 1) {
        com wn = {cos(Pi/k), inv * sin(Pi/k)};
        for (int i = 0; i < lim; i += k << 1) {
            com w = {1, 0};
            for (int j = 0; j < k; w = w * wn, j ++) {
                com x = h[i + j], y = w * h[i + j + k];
                h[i + j] = x + y;
                h[i + j + k] = x - y;
            }
        }
    }
    if (inv == 1) return ;
    for (int i = 0; i < lim; i ++) h[i].x /= lim;
}

int main () {
    n = read (f), m = read (g);
    for (lim = 1; lim <= n + m; lim <<= 1) bit ++;
    for (int i = 0; i < lim; i ++) 
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
    fft (f, 1), fft (g, 1);
    for (int i = 0; i < lim; i ++) f[i] = f[i] * g[i];
    fft (f, -1);
    for (int i = 0; i < lim; i ++) {
        ans[i] += f[i].x + 0.5;
        if (ans[i] > 9) {
            ans[i + 1] += ans[i] / 10;
            ans[i] %= 10;
        }
    }
    while (lim > 0 && ! ans[lim]) lim --;
    while (~ lim) cout << ans[lim --];
    return 0;
}

例题2：https://ac.nowcoder.com/acm/contest/26012/B

要注意的是三元组 $(i,j,k)$ 是有序的，因此直接进行FFT会导致方案数错误，考虑容斥
U = { ( i , j , k )   ∣   1 ≤ i , j , k ≤ n } A = { ( i , j , k )   ∣   i = j } B = { ( i , j , k )   ∣   j = k } C = { ( i , j , k )   ∣   i = k } U=\{(i,j,k)\ |\ 1\le i,j,k \le n\}\\ A=\{(i,j,k)\ |\ i=j\}\\ B=\{(i,j,k)\ |\ j=k\}\\ C=\{(i,j,k)\ |\ i=k\} U={(i,j,k) ∣ 1≤i,j,k≤n}A={(i,j,k) ∣ i=j}B={(i,j,k) ∣ j=k}C={(i,j,k) ∣ i=k}
那么答案为
$$\begin{array}{rl}& \frac{1}{3!}f(U-A\cup B\cup C)\\ =& \frac{1}{3!}(f(U)-f(A)-f(B)-f(C)+f(A\cap B)+f(B\cap C)+f(A\cap C)-f(A\cap B\cap C))\\ =& \frac{1}{3!}(f(U)-3f(A)+2f(A\cap B))\end{array}$$
记 $S(x)={\sum }_i{x}^{s_i},S_2(x)={\sum }_i{x}^{2s_i},S_3(x)={\sum }_i{x}^{3s_i}$ ，那么
$$\begin{array}{rl}& \frac{1}{3!}(f(U)-3f(A)+2f(A\cap B))\\ =& \frac{1}{6}(S(x{)}^3-3S_2(x)S(x)+2S_3(x))\end{array}$$

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

const int N = 4e5 + 50;
const double Pi = acos (-1.0);
const int mov = 2e4, lim = 1 << 17;

struct node {
    double x, y;
    node operator + (node r) { return {x + r.x, y + r.y}; }
    node operator - (node r) { return {x - r.x, y - r.y}; }
    node operator * (node r) { return {x * r.x - y * r.y, x * r.y + y * r.x}; }
}s[N], s2[N], s3[N], f1[N], f2[N];

int n, a[N], bit = 17, rev[N];

void fft (node* f, int flag) {
    for (int i = 0; i < lim; i ++)
        if (rev[i] > i) swap (f[i], f[rev[i]]);
    for (int k = 1; k < lim; k <<= 1) {
        node wn = {cos(Pi / k), flag * sin(Pi / k)};
        for (int i = 0; i < lim; i += k << 1) {
            node w = {1, 0};
            for (int j = 0; j < k; j ++, w = w * wn) {
                node x = f[i + j], y = w * f[i + j + k];
                f[i + j] = x + y;
                f[i + j + k] = x - y;
            }
        }
    }
    if (flag == 1) return ;
    for (int i = 0; i < lim; i ++) f[i].x /= lim;
}

int main () {
    cin >> n;
    for (int i = 1; i <= n; i ++) 
        cin >> a[i], a[i] += mov;
    for (int i = 1; i <= n; i ++) {
        s[a[i]] = {1, 0};
        s2[a[i] * 2] = {1, 0};
        s3[a[i] * 3] = {1, 0};
    }
    for (int i = 0; i < lim; i ++) 
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
    fft (s, 1), fft (s2, 1);
    for (int i = 0; i < lim; i ++) {
        f1[i] = s[i] * s[i] * s[i];
        f2[i] = s2[i] * s[i];
    }
    fft (f1, -1), fft (f2, -1);
    for (int i = 0; i < lim; i ++) {
        int cnt = (f1[i].x - 3 * f2[i].x + 2 * s3[i].x) / 6 + 0.5;
        if (cnt > 0) cout << i - mov*3 << " : " << cnt << endl;
    }
    return 0;
}

例题3：https://ac.nowcoder.com/acm/contest/26012/C

$T$ 串在匹配 $S$ 串具有容错范围 $k$，也就是说 $T[j]$ 匹配 $S[i]$，当 $S$ 串周围 $[i-k,i+k]$ 有 $T[j]$，则说明 $T$ 串与 $S[i\dots i+|T|-1]$ 串匹配。
思路： 依次考虑每一个字符 ${}^{\prime }{A}^{\prime }{,}^{\prime }{T}^{\prime }{,}^{\prime }{G}^{\prime }{,}^{\prime }{C}^{\prime }$，下面以 ${}^{\prime }{A}^{\prime }$ 为例
先将 $S$ 串进行扩展，以 $S$ 串中每一个字符 ${}^{\prime }{A}^{\prime }$ 为中心，半径为 $k$ 的范围全部改成 ${}^{\prime }{A}^{\prime }$，构成 $S^{\prime }$ 串。
那么，考虑 $S$ 串中的每一个起点 $i$ ，记 $f(i{,}^{\prime }{A}^{\prime })$ 为
$$\sum _{j=1}^{|T|}[T[j]{=}^{\prime }{A}^{\prime }]\cdot [S^{\prime }[i+j-1]{=}^{\prime }{A}^{\prime }]$$
定义 $T(x)$ 和 $S(x)$ ，这里需要将 $S^{\prime }$ 串进行翻转
$$\begin{gathered} T(x)=\sum _{j=1}^{|T|}[T[j]{=}^{\prime }{A}^{\prime }]x^j \\ S(x)=\sum _{i=1}^{|S|}[S^{\prime }[n-i]{=}^{\prime }{A}^{\prime }]x^i \end{gathered}$$
那么 $[x^u](T(x)\times S(x))$ 为
$$\sum _{j=1}^{|T|}[T[j]{=}^{\prime }{A}^{\prime }]\cdot [S^{\prime }[n-(u-j)]{=}^{\prime }{A}^{\prime }]$$
由于 $i+j-1=n-(u-j)\Rightarrow u=n-i+1$
$$f(i{,}^{\prime }{A}^{\prime })=[x^{n-i+1}](T(x)\times S(x))$$
最后判断满足下面条件的 $i$ 的数量
$$f(i{,}^{\prime }{A}^{\prime })+f(i{,}^{\prime }{T}^{\prime })+f(i{,}^{\prime }{G}^{\prime })+f(i{,}^{\prime }{C}^{\prime })=|T|$$

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

const int N = 6e5 + 50;
const double Pi = acos (-1.0);
const char ch[] = {'A', 'G', 'T', 'C'};

char s[N], t[N];
int n, m, k, a[N];
int lim, bit, rev[N];

struct node {
    double x, y;
    node (double _x=0, double _y=0): x(_x), y(_y) {}
    node operator + (node r) { return node (x + r.x, y + r.y); }
    node operator - (node r) { return node (x - r.x, y - r.y); }
    node operator * (node r) { return node (x*r.x-y*r.y, x*r.y + y*r.x); }
}f[N], g[N];

void fft (node* ff, int inv) {
    for (int i = 0; i < lim; i ++) 
        if (rev[i] > i) swap (ff[i], ff[rev[i]]);
    for (int k = 1; k < lim; k <<= 1) {
        node wn = node (cos (Pi / k), inv * sin (Pi / k));
        for (int i = 0; i < lim; i += k << 1) {
            node w = node (1, 0);
            for (int j = 0; j < k; j ++, w = w * wn) {
                node x = ff[i + j], y = w * ff[i + j + k];
                ff[i + j] = x + y;
                ff[i + j + k] = x - y;
            }
        }
    }
    if (inv == 1) return ;
    for (int i = 0; i < lim; i ++) ff[i].x /= lim;
}

void extend (node* ff) {
    static int st[N];
    for (int i = 0; i < n; i ++) st[i] = ff[i].x + 0.5;
    int cnt = 0, r = 0, l = 0;
    for (int i = 0; i < n; i ++) {
        while (r - i <= k && r < n) cnt += st[r ++];
        if (cnt) ff[i] = node (1, 0);
        cnt -= st[i];
    }
    cnt = 0, l = n - 1;
    for (int i = n - 1; i >= 0; i --) {
        while (i - l <= k && l >= 0) cnt += st[l --];
        if (cnt) ff[i] = node (1, 0);
        cnt -= st[i];
    }
}

void print (node* ff) {
    for (int i = 0; i < lim; i ++) 
        cout << int(ff[i].x + 0.5) << endl;
}

int main () {
    cin >> n >> m >> k >> s >> t;
    for (lim = 1; lim < n + m; lim <<= 1) bit ++;
    for (int i = 0; i < lim; i ++) 
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
    for (int i = 0; i < 4; i ++) {
        for (int j = 0; j < n; j ++) f[j] = node (s[j] == ch[i]);
        for (int j = n; j < lim; j ++) f[j] = node ();
        for (int j = 0; j < m; j ++) g[m - j - 1] = node (t[j] == ch[i]);
        for (int j = m; j < lim; j ++) g[j] = node ();
        extend(f);
        fft (f, 1), fft (g, 1);
        for (int j = 0; j < lim; j ++) f[j] = f[j] * g[j];
        fft (f, -1);
        for (int j = m - 1; j < lim; j ++) a[j] += f[j].x + 0.5;
    }
    int ans = 0;
    for (int i = m - 1; i < lim; i ++) ans += a[i] == m;
    cout << ans << endl;
    return 0;
}

拉格朗日插值

拉格朗日插值定理：$n$ 个点值 $(x_i,y_i)$，满足 $x_i\ne x_j,(i\ne j)$，它们唯一确定一个 $n-1$ 次多项式 $f(x)$，如下
$$f(x)=\sum _{i=1}^n{y}_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}$$

构造一个 $f(x)$，$n-1$ 次多项式，满足 $f(x_i)=y_i$。
$$f(x)=f_1(x)+f_2(x)+\cdots +f_i(x)+\cdots +f_n(x)$$
对于每一个 $f_i(x)$，要满足
$$f_i(x_k)=\{\begin{array}{llc}{y}_k& \text{if}& i=k\\ 0& \text{if}& i\ne k\end{array}$$
那么
$$\begin{array}{rl}{f}_i(x)& =A(x-x_1)\cdots (x-x_{i-1})(x-x_{i+1})\cdots (x-x_n)\\ & =A\prod _{j\ne i}(x-x_j)\end{array}$$
由于 $f_i(x_i)=y_i$，有
$$\begin{array}{rl}& f_i(x_i)=A\prod _{j\ne i}(x_i-x_j)=y_i\\ \Rightarrow & A=\frac{y_i}{{\prod }_{j\ne i}(x_i-x_j)}\end{array}$$
因此
$$\begin{array}{rl}f(x)& =f_1(x)+f_2(x)+\cdots +f_n(x)\\ & =\sum _{i=1}^n\frac{y_i}{{\prod }_{j\ne i}(x_i-x_j)}\prod _{j\ne i}(x-x_j)\\ & =\sum _{i=1}^n{y}_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}\end{array}$$
它和 $CRT$ 的关系？

• 中国剩余定理（$CRT$）：
  $$\{\begin{array}{llcc}x& \equiv & a_1& (\mathrm{m}\mathrm{o}\mathrm{d}{m}_1)\\ x& \equiv & a_2& (\mathrm{m}\mathrm{o}\mathrm{d}{m}_2)\\ & ⋮& & \\ x& \equiv & a_i& (\mathrm{m}\mathrm{o}\mathrm{d}{m}_i)\\ & ⋮& & \\ x& \equiv & a_n& (\mathrm{m}\mathrm{o}\mathrm{d}{m}_n)\end{array}$$
  其中 $(m_i,m_j)=1,(i\ne j)$，如何构造 $x$
  $$x=x_1+x_2+\cdots +x_i+\cdots +x_n$$
  对于每一个 $x_i$，要满足
  $$x_i\mathrm{m}\mathrm{o}\mathrm{d}{m}_k=\{\begin{array}{llc}{a}_i& \text{if}& i=k\\ 0& \text{if}& i\ne k\end{array}$$
  那么
  $$\begin{array}{rl}{x}_i& =A{m}_1\cdots m_{i-1}{m}_{i+1}\cdots m_n\\ & =A\prod _{j\ne i}{m}_j\end{array}$$
  由于 $x_i\mathrm{m}\mathrm{o}\mathrm{d}{m}_i=a_i$，有
  $$\begin{array}{rlr}& (x_i=A\prod _{j\ne i}{x}_j)\equiv a_i& (\mathrm{m}\mathrm{o}\mathrm{d}{m}_i)\\ \Rightarrow & A\equiv a_i(\prod _{j\ne i}{x}_j{)}^{-1}& (\mathrm{m}\mathrm{o}\mathrm{d}{m}_i)\\ \Rightarrow & A\equiv a_i{t}_i& (\mathrm{m}\mathrm{o}\mathrm{d}{m}_i)\end{array}$$
  因此
  $$\begin{array}{rlr}x& \equiv x_1+x_2+\cdots +x_n& \\ & \equiv \sum _{i=1}^n{a}_i{t}_i\prod _{j\ne i}{m}_j& (\mathrm{m}\mathrm{o}\mathrm{d}\prod _{i=1}^n{m}_i)\end{array}$$
  暴力实现：先计算 $g(x)={\prod }_{i=1}^n(x-x_i)$，再求 $n$ 次 $g(x)/(x-x_i)$，按照权值将它们相加。$O(n^2)$
  $$f(x)=\sum _{i=1}^n\frac{y_i}{{\prod }_{j\ne i}(x_i-x_j)}\frac{g(x)}{(x-x_i)}$$
  特殊情况1：只需要求一个 $f(x_0)$。$O(n^2)$
  $$f(x_0)=\sum _{i=1}^n{y}_i\prod _{j\ne i}\frac{x_0-x_j}{x_i-x_j}$$
  特殊情况2：只需要求一个 $f(x_0)$，且 $x_i=i$。$O(n)$
  $$f(x_0)=\sum _{i=1}^n{y}_i\frac{{\prod }_{j\ne i}(x_0-j)}{(-1{)}^{n-i}(n-i)!(i-1)!}$$
  预处理 ${\prod }_{j\ne i}(x_0-j)$ 的前缀积和后缀积，拼在一起。

例题4：https://ac.nowcoder.com/acm/contest/26012/D

通过观察，$f(n)={\sum }_{i=1}^n{i}^k$ 是 $k+1$ 次多项式。例如：
$$\begin{array}{rlrl}& k=0& f(n)& =n\\ & k=1& f(n)& =\frac{n(n+1)}{2}\\ & k=2& f(n)& =\frac{n(n+1)(2n+1)}{6}\\ & & ⋮& \end{array}$$
结合第二类斯特林数的公式，有
$$n^m=\sum _{k=0}^m\left\{\begin{array}{c}m\\ k\end{array}\right\}{n}^{\underline{k}}$$
其中 $n^m$ 的含义为把 $m$ 个互不相同的球放入 $n$ 个不同的盒子的方案数。换一种理解方式为枚举有 $k$ 个盒子一定要放，将 $m$ 个球都放进这 $k$ 个盒子。
$$\sum _{k=0}^m\left(\genfrac{}{}{0ex}{}{n}{k}\right)k!\left\{\begin{array}{c}m\\ k\end{array}\right\}=\sum _{k=0}^m\left\{\begin{array}{c}m\\ k\end{array}\right\}{n}^{\underline{k}}$$
那么
$$f(n)=\sum _{i=1}^n{i}^k=\sum _{i=1}^n\sum _{j=1}^k\left\{\begin{array}{c}k\\ j\end{array}\right\}{i}^{\underline{j}}=\sum _{j=1}^k\left\{\begin{array}{c}k\\ j\end{array}\right\}\sum _{i=1}^n{i}^{\underline{j}}$$
因为
$$\begin{array}{rl}\sum _{i=1}^n{i}^{\underline{j}}& =j(j-1)\cdots 1+(j+1)j\cdots 2+\cdots +n(n-1)\cdots (n-j+1)\\ & =\frac{(j+1)j\cdots 1-j(j-1)\cdots 0}{j+1}+\frac{(j+2)(j+1)\cdots 2-(j+1)j\cdots 1}{j+1}\\ & +\cdots +\frac{(n+1)n\cdots (n-j+1)-n(n-1)\cdots (n-j)}{j+1}\\ & =\frac{(n+1)n\cdots (n-j+1)}{j+1}\\ & =\frac{(n+1{)}^{\underline{j+1}}}{j+1}\end{array}$$
由此可以得知 $f(n)$ 是 $k+1$ 次多项式。那么我们需要 $k+2$ 个点，利用快速幂加前缀和 $O(klogk)$ 得到 $(0,f(0)),(1,f(1)),\cdots (k+1,f(k+1))$。再利用拉格朗日插值 $O(k)$ 得到 $f(n_0)$。

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

const int N = 1e3 + 50;
const int MOD = 9999991;

int n, m, f[N], g[N], x[N];
int pre[N], suf[N], fac, fnv[N];

int fpow (int a, int p) {
    int res = 1;
    for (; p; p >>= 1, a = (ll)a * a % MOD)
        if (p & 1) res = (ll)res * a % MOD;
    return res;
}

int lagrange (int n, int* x, int* y, int k) {
    if (k <= n) return y[k];
    int ans = 0;
    pre[0] = (k - x[0]) % MOD, suf[n + 1] = 1;
    for (int i = 1; i <= n; i ++) pre[i] = pre[i - 1] * (ll)(k - x[i]) % MOD;
    for (int i = n; i >= 0; i --) suf[i] = suf[i + 1] * (ll)(k - x[i]) % MOD;
    for (int i = 0; i <= n; i ++) {
        int up = (ll)(i == 0 ? 1 : pre[i - 1]) * suf[i + 1] % MOD;
        int down = (ll)fnv[i] * (n - i & 1 ? -1 : 1) * fnv[n - i] % MOD;
        ans = (ans + (ll)y[i] * up % MOD * down % MOD) % MOD;
    }
    return (ans % MOD + MOD) % MOD;
}

int main () {
    int T; scanf ("%d", &T);
    fac = 1;
    for (int i = 1; i < N; i ++) fac = (ll)fac * i % MOD;
    fnv[N - 1] = fpow (fac, MOD - 2);
    for (int i = N - 1; i > 0; i --) 
        fnv[i - 1] = fnv[i] * (ll)i % MOD;
    while (T --) {
        scanf ("%d%d", &n, &m);
        for (int i = 0; i <= n; i ++) {
            scanf ("%d", &f[i]);
            x[i] = i;
        }
        f[n + 1] = lagrange(n, x, f, n + 1);
        g[0] = f[0], x[n + 1] = n + 1;
        for (int i = 1; i <= n + 1; i ++) 
            g[i] = (g[i - 1] + f[i]) % MOD;
        while (m --) {
            int lef, rig; scanf ("%d%d", &lef, &rig);
            int ans = lagrange (n + 1, x, g, lef - 1);
            ans = lagrange (n + 1, x, g, rig) - ans;
            printf("%d\n", (ans % MOD + MOD) % MOD);
        }
    }
    return 0;
}

例题5：https://ac.nowcoder.com/acm/contest/26012/E

求 ${\sum }_{i=L}^{R}f(i)$，那么可以联想到 $S(R)-S(L-1)$，其中
$$S(x)=f(1)+f(2)+\cdots +f(x)=\sum _{i=1}^{x}f(i)$$
是 $n+1$ 次多项式，类似于积分，那么我们只需要 $n+2$ 个点值就可以得到 $S(x)$。并且 $O(n)$ 计算 $S(R)$ 和 $S(L-1)$。

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

const int N = 1e3 + 50;
const int MOD = 9999991;

int n, m, f[N], g[N], x[N];
int pre[N], suf[N], fac, fnv[N];

int fpow (int a, int p) {
    int res = 1;
    for (; p; p >>= 1, a = (ll)a * a % MOD)
        if (p & 1) res = (ll)res * a % MOD;
    return res;
}

int lagrange (int n, int* x, int* y, int k) {
    if (k <= n) return y[k];
    int ans = 0;
    pre[0] = (k - x[0]) % MOD, suf[n + 1] = 1;
    for (int i = 1; i <= n; i ++) pre[i] = pre[i - 1] * (ll)(k - x[i]) % MOD;
    for (int i = n; i >= 0; i --) suf[i] = suf[i + 1] * (ll)(k - x[i]) % MOD;
    for (int i = 0; i <= n; i ++) {
        int up = (ll)(i == 0 ? 1 : pre[i - 1]) * suf[i + 1] % MOD;
        int down = (ll)fnv[i] * (n - i & 1 ? -1 : 1) * fnv[n - i] % MOD;
        ans = (ans + (ll)y[i] * up % MOD * down % MOD) % MOD;
    }
    return (ans % MOD + MOD) % MOD;
}

int main () {
    int T; scanf ("%d", &T);
    fac = 1;
    for (int i = 1; i < N; i ++) fac = (ll)fac * i % MOD;
    fnv[N - 1] = fpow (fac, MOD - 2);
    for (int i = N - 1; i > 0; i --) 
        fnv[i - 1] = fnv[i] * (ll)i % MOD;
    while (T --) {
        scanf ("%d%d", &n, &m);
        for (int i = 0; i <= n; i ++) {
            scanf ("%d", &f[i]);
            x[i] = i;
        }
        f[n + 1] = lagrange(n, x, f, n + 1);
        g[0] = f[0], x[n + 1] = n + 1;
        for (int i = 1; i <= n + 1; i ++) 
            g[i] = (g[i - 1] + f[i]) % MOD;
        while (m --) {
            int lef, rig; scanf ("%d%d", &lef, &rig);
            int ans = lagrange (n + 1, x, g, lef - 1);
            ans = lagrange (n + 1, x, g, rig) - ans;
            printf("%d\n", (ans % MOD + MOD) % MOD);
        }
    }
    return 0;
}

例题6：https://ac.nowcoder.com/acm/contest/26012/F

定义 $f(i,x)$ 为 $b_i=x$ 的方案数，$[i,n]$ 已经确立，那么
$$f(i,x)=\sum _{j=a_{i+1}}^{x}f(i+1,j)$$
边界情况 $a_n\le x\le R,f(n,x)=1$，其中 a i : = m a x { a i , a i + 1 , ⋯   , a n } a_i:=max\{a_i,a_{i+1},\cdots,a_n\} ai​:=max{ai​,ai+1​,⋯,an​}。根据经验，可以知道 $f(i,x)$ 是 $n-i$ 次多项式，每次求出 $f(i+1,x)$ ，再通过拉格朗日插值求出 $f(i,x)$。
$$\begin{gathered} S(i+1,x)=\sum _{j=1}^{x}f(i+1,j) \\ \Rightarrow f(i,x)=S(i+1,x)-S(i+1,a_{i+1}-1) \end{gathered}$$

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

const int N = 5e3 + 50;
const int MOD = 1e9 + 7;

int f[N], pre[N], suf[N];
int n, r, a[N], fnv[N], x[N];

int fpow (int a, int p) {
    int res = 1;
    for (; p; p >>= 1, a = (ll)a * a % MOD)
        if (p & 1) res = (ll)res * a % MOD;
    return res;
}

int lagrange (int n, int* x, int* y, int k) {
    int ans = 0;
    pre[0] = (k - x[0]), suf[n + 1] = 1;
    for (int i = 1; i <= n; i ++) pre[i] = (ll)pre[i - 1] * (k - x[i]) % MOD;
    for (int i = n; i >= 0; i --) suf[i] = (ll)suf[i + 1] * (k - x[i]) % MOD;
    for (int i = 0; i <= n; i ++) {
        int up = (ll)(i == 0 ? 1 : pre[i - 1]) * suf[i + 1] % MOD;
        int down = (ll)fnv[i] * fnv[n - i] * (n - i & 1 ? -1 : 1) % MOD;
        ans = (ans + (ll)y[i] * down % MOD * up % MOD) % MOD;
    }
    return (ans % MOD + MOD) % MOD;
}

int main () {
    scanf ("%d%d", &n, &r);
    for (int i = 1; i <= n; i ++) 
        scanf ("%d", &a[i]);
    for (int i = n - 1; i >= 1; i --) 
        a[i] = max (a[i], a[i + 1]);
    for (int i = 0; i <= n; i ++) f[i] = 1;
    for (int i = 0; i < N; i ++) x[i] = i;
    int fac = 1;
    for (int i = 1; i < N; i ++) fac = (ll)fac * i % MOD;
    fnv[N - 1] = fpow (fac, MOD - 2);  
    for (int i = N - 1; i >= 0; i --) fnv[i - 1] = (ll)fnv[i] * i % MOD;
    for (int i = n - 1; i >= 0; i --) {
        for (int j = 1; j <= n; j ++) {
            f[j] = (f[j] + f[j - 1]) % MOD;
        }
        int temp = lagrange (n - i, x, f, a[i + 1] - 1);
        for (int j = 0; j <= n; j ++) {
            f[j] = (f[j] - temp) % MOD;
        }
    }
    printf("%d\n", lagrange(n, x, f, r));
    return 0;
}