FFT与NTT专题

先不管旋转操作，考虑化简这个差异值

$$
begin{aligned}
sum_{i=1}^n(x_i-y_i-c)^2
&=sum_{i=1}^n(x_i-y_i)^2+nc^2-2csum_{i=1}^n(x_i-y_i)
&=sum_{i=1}^nx_i^2+sum_{i=1}^ny_i^2+nc^2-2csum_{i=1}^n(x_i-y_i)-2sum_{i=1}^nx_iy_i
end{aligned}
$$

注意到$sum x^2+sum y^2$是常数，先不管

可以发现，这是一个关于$c$的二次函数

那么我们知道，此时$c$的极值点就在$-frac{b}{2a}$处

所以，我们可以得出$c$的最优值是

$$
frac{sum_{i=1}^n x_i-sum_{i=1}^n y_i}{n}
$$

而分子的两个数均与旋转无关

但是$c$只能是整数

所以判一下$c, c-1, c+1$哪个与上面的式子更接近

注意到旋转唯一能改变的是$sum xy$

而我们要让这个值尽量小

设

$$
F(m)=sum_{i=1}^nx_iy_{i+m}
$$

我们可以看出，这是一个类似卷积的东西

但是一般的卷积是后两式下标的和不变

而这个是差不变

所以把这个式子变一下

设

$$
x_{n-i+1}=x_i
$$

就是将x倒序一下

可以得到

$$
F(m)=sum_{i=1}^nx_{n-i+1}y_{i+m}
$$

不妨设后面$xy$的卷积是$A$，也就是

$$
A(n+m+1)=sum_{i=1}^nx_{n-i+1}y_{i+m}
$$

可以发现，这个$A$就是将$F$整体向右平移了$n+1$

所以

$$
F(m)=A(n+m+1)
$$

为了不丢精度，NTT即可（保证答案不会超过mod）

代码如下 

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

using namespace std;
#define N 100010
#define LL long long
int r[N << 2]; const int mod = 998244353;
inline int (int x, int y) {
    int res = 1;
    for (;y;y >>= 1, x = (LL)x * x % mod) if (y & 1) res = (LL)res * x % mod;
    return res;
}
inline void NTT(int len, int type, int a[]) {
    for (int i = 0;i < len;i++) if (i < r[i]) swap(a[i], a[r[i]]);
    for (int mid = 2;mid <= len;mid <<= 1) {
        int Wn = Pow(3, type ? (mod - 1) / mid : mod - 1 - (mod - 1) / mid);
        for (int i = 0;i < len;i += mid)
            for (int j = i, w = 1, t;j < i + (mid >> 1);j++, w = (LL)w * Wn % mod)
                t = (LL)w * a[j + (mid >> 1)] % mod, a[j + (mid >> 1)] = (a[j] - t + mod) % mod, a[j] = (a[j] + t) % mod;
    }
}
int x[N], y[N], A[N << 2], B[N << 2], res[N];
template<class T> inline T Abs(const T x) {return x > 0 ? x : -x;}
int main() {
    int n, m, sumx = 0, sumy = 0, sumx2 = 0, sumy2 = 0; scanf("%d%d", &n, &m);
    for (int i = 1;i <= n;i++) scanf("%d", &x[i]), A[i] = x[i], sumx += x[i], sumx2 += x[i] * x[i];
    for (int i = 1;i <= n;i++) scanf("%d", &y[i]), B[2 * n - i] = B[n - i] = y[i], sumy += y[i], sumy2 += y[i] * y[i];
    int len = 1, l = 0;
    while (len <= 3 * n) len <<= 1, l++;
    for (int i = 1;i <= len;i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << l - 1);
    NTT(len, 1, A), NTT(len, 1, B);
    for (int i = 0;i < len;i++) A[i] = (LL)A[i] * B[i] % mod;
    NTT(len, 0, A); int Inv = Pow(len, mod - 2);
    for (int i = 0;i < n;i++) res[i] = (LL)A[2 * n - i] * Inv % mod;
    int c = (sumx - sumy) / n; LL ans = 1e18;
    if (Abs((LL)c * n - sumx + sumy) > Abs(((LL)c * n + n - sumx + sumy))) c++;
    if (Abs((LL)c * n - sumx + sumy) > Abs(((LL)c * n - n - sumx + sumy))) c--;
    for (int i = 0;i < n;i++) {
        LL tmp = (LL)sumx2 + sumy2 - 2 * res[i] - (LL)2 * c * (sumx - sumy) + (LL)n * c * c;
        if (tmp < ans) ans = tmp;
    }
    printf("%lldn", ans);
    return 0;
}

B – 求和

我们知道

$$
S(n,m)=sum_{i=0}^m(-1)^i{mchoose i}(m-i)^nfrac{1}{m!}
$$

原理很简单，容斥有几个盒子没有放球，有$mchoose i$种选法，再将$n$个球放入$m-i$个盒子。由于盒子是无序的，最后除以$m$的阶乘

那么我们用这个化简原式

注意到第二个$sum$的上界是$i$，非常讨厌

由于斯特林数的性质，把这个$i$换成$n$也没有问题

因为当$m>n$时，$S(n,m)=0$

所以有

$$
begin{aligned}
sum_{i=0}^nsum_{j=0}^nS(i,j)*2^j*j!
&=sum_{j=0}^n2^j*j!sum_{i=0}^nsum_{k=0}^j(-1)^k{jchoose k}(j-k)^ifrac{1}{j!}
&=sum_{j=0}^n2^j*j!sum_{k=0}^jfrac{(-1)^k}{k!}*frac{sum_{i=0}^n(j-k)^i}{(j-k)!}
end{aligned}
$$

注意到后面那个是卷积的形式

第一个多项式很好求，第二个的分子是等比数列

我们设$B$是第二个多项式

显然有

$$
B(0)=0, B(1)=n+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
26
27
28
29
30
31
32
33
34
35
36
37

using namespace std;
#define N 100010
#define LL long long
int r[N << 2]; const int mod = 998244353;
inline int (int x, int y) {
    int res = 1;
    for (;y;y >>= 1, x = (LL)x * x % mod) if (y & 1) res = (LL)res * x % mod;
    return res;
}
inline void NTT(int len, int type, int a[]) {
    for (int i = 0;i < len;i++) if (i < r[i]) swap(a[i], a[r[i]]);
    for (int mid = 2;mid <= len;mid <<= 1) {
        int Wn = Pow(3, type ? (mod - 1) / mid : mod - 1 - (mod - 1) / mid);
        for (int i = 0;i < len;i += mid)
            for (int j = i, w = 1, t;j < i + (mid >> 1);j++, w = (LL)w * Wn % mod)
                t = (LL)w * a[j + (mid >> 1)] % mod, a[j + (mid >> 1)] = (a[j] - t + mod) % mod, a[j] = (a[j] + t) % mod;
    }
}
int A[N << 2], B[N << 2], frac[N];
int main() {
    int n, len = 1, l = 0; scanf("%d", &n);
    frac[0] = 1;
    for (int i = 1;i <= n;i++) frac[i] = (LL)frac[i - 1] * i % mod;
    while (len <= 2 * n) len <<= 1, l++;
    for (int i = 0;i < len;i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << l - 1);
    A[0] = B[0] = 1, B[1] = n + 1;
    for (int i = 1;i <= n;i++) A[i] = (i & 1 ? -1 : 1) * Pow(frac[i], mod - 2), A[i] = (A[i] + mod) % mod;
    for (int i = 2;i <= n;i++) B[i] = ((LL)(Pow(i, n + 1) - 1) * Pow(i - 1, mod - 2) % mod * Pow(frac[i], mod - 2) % mod + mod) % mod;
    NTT(len, 1, A), NTT(len, 1, B);
    for (int i = 0;i < len;i++) A[i] = (LL)A[i] * B[i] % mod;
    NTT(len, 0, A); int Inv = Pow(len, mod - 2);
    int tmp = 1, res = 0;
    for (int i = 0;i <= n;i++) res = (res + (LL)tmp * frac[i] % mod * A[i]) % mod, tmp = tmp * 2 % mod;
    printf("%dn", (LL)res * Inv % mod);
    return 0;
}

C – 序列统计

这题的难点在于转化成原根

注意到要求的是所有数的乘积而非和

如果是和的话直接NTT就好了

那么我们就将乘积转化成和的形式

如果两个数都是某个数的某次方，那么这两个数乘起来就是指数相加

而原根恰好可以表示模$m$剩余系下的每个数

所以把每个数转化成原根的某次方就好了

求原根代码

1
2
3
4
5
6
7
8
inline int G(int x) {
    if (x == 2) return 1;
    for (int i = 2, flg = 1;i;i++, flg = 1) {
        for (int j = 2;j * j < x;j++) 
            if ((x - 1) % j == 0 && Pow(i, (x - 1) / j, x) == 1) {flg = 0; break;}
        if (flg) return i;
    }
}

D – 残缺的字符串

带通配符的字符串匹配问题

首先考虑不带通配符的怎么做

那么拓展KMP， 后缀数组都可以

但是我们有一个更高级的方法：FFT求字符串匹配

首先我们需要定义“匹配”

所以设差异函数$g(i)$表示从$B$串的$i$位置开始，与$A$串的差异程度

有

$$
g(x)=sum_{i=x}^{x+m-1}(B_i-A_{i-x+1})^2
$$

显然，只有当$A$串从$x$位置开始与$B$串完全相同，$g$的值才为0

化简原式

$$
begin{aligned}
g(x)&=sum_{i=x}^{x+m-1}A_{i-x+1}^2+sum_{i=x}^{x+m-1}B_i^2-2sum_{i=x}^{x+m-1}B_iA_{i-x+1}
&=sum_{i=1}^mA_i^2+sum_{i=1}^mB_{i+x-1}^2-2sum_{i=1}^mA_iB_{i+x-1}
end{aligned}
$$

前两项可以通过预处理前缀和得出，后面的是一个下标差相等的卷积

那么模仿之前的套路，我们将$A$序列倒序一下再求卷积就行了

解决了不带通配符的问题，再考虑带通配符

这个通配符是可以匹配任意字符的，所以把差异函数改一下

$$
g(x)=sum_{i=x}^{x+m-1}(B_i-A_{i-x+1})^2A_{i-x+1}B_i
$$

当$i$处的字符是$*$时，我们设那个地方的值为0

化简得

$$
=sum_{i=1}^mA_i^3B_{i-x+1}+sum_{i=1}^mA_iB_{i-x+1}^3-2sum_{i=1}^xA_i^2B_{i-x+1}^2
$$

做3次FFT即可

代码如下

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

using namespace std;
#define N 300010
int r[N << 2]; const double PI = acos(-1);
inline void FFT(int len, int type, complex<double> a[]) {
    for (int i = 1;i < len;i++) if (i < r[i]) swap(a[i], a[r[i]]);
    for (int mid = 2;mid <= len;mid <<= 1) {
        complex<double> Wn(cos(2 * PI / mid), type * sin(2 * PI / mid));
        for (int i = 0;i < len;i += mid) {
            complex<double> w(1), t;
            for (int j = i;j < i + (mid >> 1);j++, w *= Wn) 
                t = w * a[j + (mid >> 1)], a[j + (mid >> 1)] = a[j] - t, a[j] += t;
        }
    }
}
char a[N], b[N]; complex<double> A[N << 2], B[N << 2]; int a1[N], b1[N];
#define LL long long
LL res[N]; vector<int> ans;
#define Clear(x) for (int i = 0;i < len;i++) x[i] = 0;
inline void mul(int len, int n, int k) {
    FFT(len, 1, A), FFT(len, 1, B);
    for (int i = 0;i < len;i++) A[i] *= B[i];
    FFT(len, -1, A);
    for (int i = 1;i <= n;i++) res[i] += k * (LL)(A[i].real() / len + 0.5);
}
int main() {
    int n, m; scanf("%d%d%s%s", &m, &n, a + 1, b + 1);
    for (int i = 1;i <= m;i++) a1[m - i] = a[i] == '*' ? 0 : a[i] - 'a' + 1;
    for (int i = 1;i <= n;i++) b1[i] = b[i] == '*' ? 0 : b[i] - 'a' + 1;
    int len = 1, l = 0;
    while (len <= m + n) len <<= 1, l++;
    for (int i = 0;i < len;i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << l - 1);
    for (int i = 0;i < m;i++) A[i] = a1[i] * a1[i] * a1[i];
    for (int i = 1;i <= n;i++) B[i] = b1[i];
    mul(len, n, 1); Clear(A) Clear(B)
    for (int i = 0;i < m;i++) A[i] = a1[i];
    for (int i = 1;i <= n;i++) B[i] = b1[i] * b1[i] * b1[i];
    mul(len, n, 1); Clear(A) Clear(B)
    for (int i = 0;i < m;i++) A[i] = a1[i] * a1[i];
    for (int i = 1;i <= n;i++) B[i] = b1[i] * b1[i];
    mul(len, n, -2);
    for (int i = m;i <= n;i++)
        if (res[i] == 0) ans.push_back(i - m + 1);
    printf("%dn", ans.size());
    for (auto i : ans) printf("%d ", i);
    return 0;
}

E – 万径人踪灭

假设当前确定了一个对称中心$i$

那么当两个位置$j,k$关于i对称且这两个位置的字母相同时对答案有贡献

对称则意味着$j+k=i*2​$，可以FFT

枚举字符，然后FFT

假设这个中心有x对这样的位置

那么每一对都是独立的，可以选也可以不选，但是不能都不选

所以此时的答案为$2^x-1$

题目要求不能全部连续，那么最后再跑一边manacher，减去全部连续的答案即可

代码如下

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

using 大专栏  FFT与NTT专题s="keyword">namespace std;
#define N 100010
#define LL long long
int r[N << 2]; const int mod = 998244353;
inline int (int x, int y, int p = mod) {
    int res = 1;
    for (;y;y >>= 1, x = (LL)x * x % p) if (y & 1) res = (LL)res * x % p;
    return res;
}
inline void NTT(int len, int type, int a[]) {
    for (int i = 0;i < len;i++) if (i < r[i]) swap(a[i], a[r[i]]);
    for (int mid = 2;mid <= len;mid <<= 1) {
        int Wn = Pow(3, type ? (mod - 1) / mid : mod - 1 - (mod - 1) / mid);
        for (int i = 0;i < len;i += mid)
            for (int j = i, t, w = 1;j < i + (mid >> 1);j++, w = (LL) w * Wn % mod)
                t = (LL)w * a[j + (mid >> 1)] % mod, a[j + (mid >> 1)] = (a[j] - t + mod) % mod, a[j] = (a[j] + t) % mod;
    }
}
int A[N << 2], B[N << 2], pal[N << 2], Pow2[N]; char s[N], t[N * 2]; const int p = 1e9 + 7;
inline int manacher(int n) {
    t[0] = '#', t[n * 2 + 1] = '$', t[n * 2 + 2] = '@';
    for (int i = 1;i <= n;i++) t[i * 2 - 1] = '$', t[i * 2] = s[i]; 
    int pos = 1, mx = 1, res = 0; pal[1] = 1; 
    for (int i = 2;i <= n * 2;i++) {
        if (i <= mx) pal[i] = min(mx - i + 1, pal[2 * pos - i]);
        else pal[i] = 1;
        while (t[i - pal[i]] == t[i + pal[i]]) pal[i]++;
        if (i + pal[i] - 1 > mx) mx = i + pal[i] - 1, pos = i;
        res = (res + pal[i] / 2) % p;
    }
    return res;
}
int main() {
    scanf("%s", s + 1); int n = strlen(s + 1);
    for (int i = 1;i <= n;i++) if (s[i] == 'a') A[i] = 1; else B[i] = 1;
    int len = 1, l = 0, ans = 0; Pow2[0] = 1;
    for (int i = 1;i <= n;i++) Pow2[i] = Pow2[i - 1] * 2 % p;
    while (len <= n * 2) len <<= 1, l++;
    for (int i = 0;i < len;i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << l - 1);
    NTT(len, 1, A), NTT(len, 1, B);
    for (int i = 0;i < len;i++) A[i] = (LL)A[i] * A[i] % mod, B[i] = (LL)B[i] * B[i] % mod;
    NTT(len, 0, A), NTT(len, 0, B); int Inv = Pow(len, mod - 2);
    for (int i = 2, t;i <= n * 2;i++) {
        t = (LL)(A[i] + B[i]) * Inv % mod + (i & 1 ^ 1);
        ans = (ans + Pow2[t / 2] - 1) % p;
    }
    printf("%dn", (ans - manacher(n) + p) % p);
    return 0;
}

F – 性能优化

这道题利用到了FFT的原理

如果模数是质数，那么非常好办

但是这题不仅模数不是质数，而且求的是循环卷积，直接FFT会爆炸

贴一篇我觉得很好的题解

img

这个rev数组可以模拟FFT的过程，递归地求出来

单位根满足消去律，上面的$F(omega_n^i)$指的是$ileq frac{n}{p}$的情况

对于剩余的情况，有$omega_{frac{n}{p}}^i=omega_{frac{n}{p}}^{i-frac{n}{p}}$

也就是说，代入的$F^{[0]},F^{[1]},cdots,F^{[p-1]}$都相同，但是系数不同

然后分治就可以了

同样地，最后需要除以len，也就是模数$-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
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

using namespace std;
#define LL long long
#define N 500010
inline int (int x, int y, int mod) {
    int res = 1;
    for (;y;y >>= 1, x = (LL)x * x % mod) if (y & 1) res = (LL)res * x % mod;
    return res;
}
int tot, prime[N], r[N];
int GetPos(int x, int dep, int len, int cnt) {
    if (dep == tot + 1) return cnt;
    int tmp = len / prime[dep], s = x % prime[dep];
    return GetPos((x - s) / prime[dep], dep + 1, tmp, cnt + tmp * s);
}
int tmp[N], g;
inline void NTT(int len, int a[], int mod, int type) {
    for (int i = 0;i < len;i++) tmp[r[i]] = a[i];
    for (int i = 0;i < len;i++) a[i] = tmp[i];
    for (int i = tot, block = 1;i >= 1;i--) {
        int mid = block; block *= prime[i];
        int Wn = Pow(g, type ? (mod - 1) / block : mod - 1 - (mod - 1) / block, mod);
        for (int j = 0;j < len;j++) tmp[j] = 0;
        for (int j = 0, Wk = 1;j < len;j += block, Wk = 1)
            for (int k = 0;k < block;k++) {
                for (int l = k % mid, w = 1;l < block;l += mid, w = (LL)w * Wk % mod)
                    tmp[j + k] = (tmp[j + k] + (LL)w * a[j + l]) % mod;
                Wk = (LL)Wk * Wn % mod;
            }
        for (int j = 0;j < len;j++) a[j] = tmp[j];
    } 
}
inline int GetG(int x) {
    if (x == 2) return 1;
    for (int i = 2, flag = 1;i;i++, flag = 1) {
        for (int j = 2;j * j < x;j++)
            if (Pow(i, (x - 1) / j, x) == 1) {flag = 0; break;}
        if (flag == 1) return i;
    }
}
int A[N << 2], B[N << 2];
int main() {
    int n, C; scanf("%d%d", &n, &C);
    for (int i = 0;i < n;i++) scanf("%d", &A[i]);
    for (int i = 0;i < n;i++) scanf("%d", &B[i]);
    int tmp = n; g = GetG(n + 1);
    for (int i = 2;i * i <= tmp;i++)
        while (tmp % i == 0) prime[++tot] = i, tmp /= i;
    if (tmp != 1) prime[++tot] = tmp;
    for (int i = 0;i < n;i++) r[i] = GetPos(i, 1, n, 0);
    NTT(n, A, n + 1, 1), NTT(n, B, n + 1, 1);
    for (int i = 0;i < n;i++) A[i] = (LL)A[i] * Pow(B[i], C, n + 1) % (n + 1);
    NTT(n, A, n + 1, 0); int Inv = Pow(n, n - 1, n + 1);
    for (int i = 0;i < n;i++) printf("%dn", (LL)A[i] * Inv % (n + 1));
    return 0;
}

H – Frightful Formula

算是比较简单的一道题

公式等价于一个表格，往右走有$a$种方法，往下走有$b$种方法，还可以直接从这个格子开始走，有$c$种方法

先不考虑第一行和第一列格子

假设是从$i,j$这个格子开始走的

那么，这个格子需要向右走$n-j$步，向下走$n-i$步

对答案的贡献是

$$
c*a^{n-i}*b^{n-j}*{n-i+n-jchoose n-i}
$$

含义是，从这个格子开始，有$c$种走法，向有走$n-j$次，向下走$n-i$次，在$n-j+n-i$步中，有$n-i​$步是往下走的

那么，把这些空白的格子加起来，我们可以得到

$$
begin{aligned}
csum_{i=2}^nsum_{j=2}^na^{n-i}b^{n-j}{n-i+n-jchoose n-i}
&=csum_{i=0}^{n-2}sum_{j=0}^{n-2}a^ib^j{i+jchoose i}
&=csum_{i=0}^{n-2}frac{a^i}{i!}sum_{j=0}^{n-2}(i+j)!frac{b^j}{j!}
end{aligned}
$$

我们可以枚举$i$，后面的是一个下标差相等的卷积

将多项式逆序一下就可以了

这道题没有给模数，而答案又很大

为了防止丢精度，所以使用MTT

代码如下

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
#include <bits/stdc++.h>
using namespace std;
#define N 200010
#define LL long long
const int mod = 1e6 + 3;
inline int (int x, int y, int p = mod) {
    int res = 1;
    for (;y;y >>= 1, x = (LL)x * x % p) if (y & 1) res = (LL)res * x % p;
    return res;
}
const int ZJK = (1 << 19) + 233;
int frac[N * 2], f[2][N];
inline int C(int n, int r) {return (LL)frac[n] * Pow(frac[r], mod - 2) % mod * Pow(frac[n - r], mod - 2) % mod;}

struct CP {
    double x, y;
    CP(double _x = 0, double _y = 0) : x(_x), y(_y) {}
    CP operator * (const CP &b) {return CP(x * b.x - y * b.y, x * b.y + y * b.x);}
    CP operator + (const CP &b) {return CP(x + b.x, y + b.y);}
    CP operator - (const CP &b) {return CP(x - b.x, y - b.y);}
    CP operator / (const double b) {return CP(x / b, y / b);}
};
int r[ZJK]; const double PI = acos(-1);
inline void FFT(int len, int type, CP a[]) {
    for (int i = 0;i < len;i++) if (i < r[i]) swap(a[i], a[r[i]]);
    for (int mid = 2;mid <= len;mid <<= 1) {
        CP Wn(cos(2 * PI / mid), type * sin(2 * PI / mid));
        for (int i = 0;i < len;i += mid) {
            CP w(1), t;
            for (int j = i;j < i + (mid >> 1);j++, w = w * Wn) 
                t = w * a[j + (mid >> 1)], a[j + (mid >> 1)] = a[j] - t, a[j] = a[j] + t; 
        }
    }
    if (type == -1) for (int i = 0;i < len;i++) a[i] = a[i] / len;
}
#define LL long long
CP a[ZJK], b[ZJK], c[ZJK], d[ZJK];
int A[ZJK], B[ZJK];
inline void MTT(int len, LL m) {
    for (int i = 0;i < len;i++) a[i] = A[i] / m, b[i] = A[i] % m, c[i] = B[i] / m, d[i] = B[i] % m;
    FFT(len, 1, a), FFT(len, 1, b), FFT(len, 1, c), FFT(len, 1, d);
    for (int i = 0;i < len;i++) {
        CP a1 = a[i], b1 = b[i], c1 = c[i], d1 = d[i];
        a[i] = a1 * c1, b[i] = a1 * d1, c[i] = b1 * c1, d[i] = b1 * d1;
    }
    FFT(len, -1, a), FFT(len, -1, b), FFT(len, -1, c), FFT(len, -1, d);
    for (int i = 0;i < len;i++) 
        A[i] = ((LL)(a[i].x + 0.5) * m % mod * m % mod + (LL)(b[i].x + 0.5) * m % mod + (LL)(c[i].x + 0.5) * m % mod
            + (LL)(d[i].x + 0.5)) % mod;
}
int main() {
    int a, b, c, n; scanf("%d%d%d%d", &n, &a, &b, &c), frac[0] = 1;
    for (int i = 1;i <= n;i++) scanf("%d", &f[1][i]);
    for (int i = 1;i <= n;i++) scanf("%d", &f[0][i]);
    for (int i = 1;i <= n * 2;i++) frac[i] = (LL)frac[i - 1] * i % mod;
    int ans = 0, tmp1 = Pow(b, n - 1), tmp2 = Pow(a, n - 1);
    for (int i = 2;i <= n;i++) ans = (ans + (LL)f[0][i] * C(2 * n - 2 - i, n - i) % mod * tmp1 % mod * Pow(a, n - i)) % mod;
    for (int i = 2;i <= n;i++) ans = (ans + (LL)f[1][i] * C(2 * n - 2 - i, n - i) % mod * tmp2 % mod * Pow(b, n - i)) % mod;
    int len = 1, l = 0;
    while (len <= 2 * n) len <<= 1, l++;
    for (int i = 0;i < len;i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << l - 1);
    for (int i = 0, t;i <= n - 2;i++) t = Pow(frac[i], mod - 2), A[i] = (LL)Pow(a, i) * t % mod, B[i] = (LL)Pow(b, i) * t % mod;
    MTT(len, 1000);
    for (int i = 0;i <= 2 * n - 4;i++) ans = (ans + (LL)c * A[i] % mod * frac[i]) % mod;
    printf("%dn", ans); 
    return 0;
}
【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/lijianming180/p/12366543.html