zoukankan      html  css  js  c++  java

  • 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;
    }
  • 相关阅读:
    POJ1579Function Run Fun
    C++ 程序员必读书目清单
    zoj2100Seeding(水题)
    接口开发及技术负责
    哪些需求最重要
    地址
    哪些需求最重要
    setTimeOut与 setInterval区别
    项目管理简介
    项目管理简介
  • 原文地址：https://www.cnblogs.com/lijianming180/p/12366543.html
Copyright © 2011-2022 走看看