首页 > 编程语言> > 狂恋多项式算法(FFT,NTT,生成函数,插值,求逆...)

狂恋多项式算法(FFT,NTT,生成函数,插值,求逆...)

2021-07-25 12:03:17 作者：互联网

咕咕咕~~

P3803 【模板】多项式乘法（FFT）(MODULE + FFT)

模板

#include <bits/stdc++.h>
using namespace std;
static constexpr int N = 4000010;
static constexpr double eps = 1e-6;
struct cp {
	double r, i;
	cp(double _r = 0, double _i = 0) : r(_r), i(_i) { }
	cp operator+(const cp x) { return cp(r + x.r, i + x.i); }
	cp operator-(const cp x) { return cp(r - x.r, i - x.i); }
	cp operator*(const cp x) { return cp(r * x.r - i * x.i, r * x.i + i * x.r); }
};
static const double pi = acos(-1.0);
void FFT(cp a[], int n, int t = 1) {
    for (int i = 1, j = 0; i < n - 1; ++i) {
        for (int s = n; j ^= s >>= 1, ~j & s; );
        if (i < j) swap(a[i], a[j]);
    }
    for (int d = 0; (1 << d) < n; ++d) {
        int m = 1 << d, m2 = m << 1;
        double o = pi / m * t;
        cp _w(cos(o), sin(o));
        for (int i = 0; i < n; i += m2) {
            cp w(1, 0);
            for (int j = 0; j < m; ++j) {
                cp&A = a[i + j + m], &B = a[i + j], t = w * A;
                A = B - t; B = B + t; w = w * _w;
            }
        }
    }
    if (t == -1) for (int i = 0; i < n; ++i) a[i].r /= n;
}
cp a[N << 1], b[N];
int main() {
    int n, m;
    cin >> n >> m;
    for (int i = 0; i <= n; ++i) cin >> a[i].r;
    for (int i = 0; i <= m; ++i) cin >> b[i].r;

    int siz = 1;
    while ((siz <<= 1) <= n + m);

    FFT(a, siz), FFT(b, siz);
    for (int i = 0; i <= siz; ++i) a[i] = a[i] * b[i];
    FFT(a, siz, -1);

    for (int i = 0; i <= n + m; ++i) {
        printf("%.0f ", a[i].r + eps);
    }
    return 0;
}

P1919 【模板】A*B Problem升级版（FFT快速傅里叶）(MODULE + FFT + 高精度乘法)

把这两个大数的每一位看做多项式的系数，两多项式的乘法就相当于没有进位的数位乘法，可以用 F F T \rm FFT FFT在 O ( ( n + m ) log ⁡ ( n + m ) ) \Omicron((n+m)\log(n+m)) O((n+m)log(n+m))的时间处理。

#include <bits/stdc++.h>
using namespace std;
static constexpr int N = 4000010;
static constexpr double eps = 0.5;
// （模板略）
int ans[N], stk[N], top = 1;
cp _a[N << 1], _b[N];
int main() {
    ios::sync_with_stdio(false), cin.tie(nullptr);
    string a, b;
    cin >> a >> b;
    int n1 = a.length(), n2 = b.length();
    for (int i = 0; i < n1; ++i) _a[i].r = a[n1 - i - 1] - '0';
    for (int i = 0; i < n2; ++i) _b[i].r = b[n2 - i - 1] - '0';
    int n = 1;
    while ((n <<= 1) <= n1 + n2);// n <<= 1;
    
    FFT(_a, n), FFT(_b, n);
    for (int i = 0; i <= n; ++i) _a[i] = _a[i] * _b[i];
    FFT(_a, n, -1);
    
    for (int i = 0; i <= n1 + n2; ++i) {
        ans[i] += int(_a[i].r + eps);
        ans[i + 1] += ans[i] / 10;
        ans[i] %= 10;
        stk[++top] = ans[i];
    }
    while (!stk[top]) top--;
    while (top != 1) {
        cout << stk[top];
        top --;
    }
    return 0;
}

P2553 [AHOI2001]多项式乘法(字符串操作 + FFT)

这个题的难度不在与多项式算法，而是字符串模拟…网上看了好多题解大概是只含有*或者是只含有空格的一个串也有可能出现…可以说是非常毒瘤了…

#include <bits/stdc++.h>
using namespace std;
static constexpr int N = 1000010;
static constexpr double eps = 1e-6;
// (模板略)
char poly[N];
cp a1[N], a2[N], ans[N];
int len, i = 1, ind; // ind 表示 index --> 当前已经输入的多项式的度 deg 之和
void setpoly(cp a[]) {
    while (poly[i] != ')' && i < len) {
        int x = 0, z = 0;
        while (!isdigit(poly[i]))
            i++;
        x = poly[i] - '0';
        if (isdigit(poly[i + 1])) {
            i++;
            x = x * 10 + poly[i] - '0';
        }
        i++;
        if (poly[i] == 'a') {
            i++;
            i++;
            z = poly[i] - '0';
            if (isdigit(poly[i+1])) {
                i++;
                z = z * 10 + poly[i] - '0';
            }
            i++;
        }
        ind ++;
        a[z].r += x;
    }
}
int main() {
    while (gets(poly), poly[0]) {
        len = strlen(poly);
        memset(a1, 0, sizeof(a1));
        memset(a2, 0, sizeof(a2));
        memset(ans, 0, sizeof(ans));
        i = 1, ind = 0;
        setpoly(a1);
        i++;
        setpoly(a2);

        int siz = 1;
        while ((siz <<= 1) <= ind);
        
        FFT(a1, siz), FFT(a2, siz);
        for (int i = 0; i <= siz; ++i) ans[i] = a1[i] * a2[i];
        FFT(ans, siz, -1);

        bool first = false;
        for (int j = siz; j >= 0; j--) {
            if (int(ans[j].r + eps) == 0)
                continue;
            if (first)
                putchar('+');
            else
                first = true;
            if (j > 0)
                printf("%.0fa^%d", ans[j].r + eps, j);
            else
                printf("%.0f", ans[j].r + eps);
        }
        putchar(10);
        memset(poly, 0, sizeof(poly));
    }
    return 0;
}

BZOJ3771. Triple(生成函数+容斥原理+FFT)

这类问题，也即组合计数基本都可以构造普通型生成函数来解决。

排列计数常使用构造指数型生成函数解决。

定义 a a a数列： a i a_i ai表示「价值为 i i i的物品的出现次数」(朴素的想就是选了 a i a_i ai次价值为 i i i的元素)
对 a a a构造普通型生成函数 A ( x ) = a 0 x 0 + a 1 x 1 + a 2 x 2 + ⋯ + a n x n A(x) = a_0x^0+a_1x^1+a_2x^2+\cdots+a_nx^n A(x)=a0x0+a1x1+a2x2+⋯+anxn，表示 x 元素选了 1 次的方案数 x元素选了1次的方案数 x元素选了1次的方案数
如果不计重复的，选的方案数就是 A ( x ) + A 2 ( x ) + A 3 ( x ) A(x) + A^2(x)+A^3(x) A(x)+A2(x)+A3(x)

n 排列数 n排列数 n排列数的生成函数表示为 n n n个生成函数的积。

为了叙述方便，我们定义 B ( x ) 为 x 元素选了 2 次的方案数 B(x)为x元素选了2次的方案数 B(x)为x元素选了2次的方案数， C ( x ) 为 x 元素选了 3 次的方案数 C(x)为x元素选了3次的方案数 C(x)为x元素选了3次的方案数
然而题目中说需要将重复的减去，这里就是容斥原理的部分了：选一次没有重复；选两次中可能选到同一元素，我们减去 B ( x ) B(x) B(x)，同时考虑到无序 ( a , b ) = ( b , a ) (a,b)=(b,a) (a,b)=(b,a)，还要除以 2 ! = 2 2! = 2 2!=2；选三次也是类似的想法： ( a , a , b ) , ( a , b , a ) , ( b , a , a ) ) c y c \Big(a,a,b),(a,b,a),(b,a,a)\Big)_{cyc} (a,a,b),(a,b,a),(b,a,a))cyc，减去 3 A ( x ) B ( x ) 3A(x)B(x) 3A(x)B(x)，然而 3 A ( x ) B ( x ) 3A(x)B(x) 3A(x)B(x)中包含了 3 ( a , a , a ) c y c 3(a,a,a)_{cyc} 3(a,a,a)cyc,原本的 A 3 ( x ) A^3(x) A3(x)中有 1 1 1份 ( a , a , a ) c y c (a,a,a)_{cyc} (a,a,a)cyc，加上这部分 2 C ( x ) 2C(x) 2C(x)最后除以 3 ! = 6 3! = 6 3!=6即可。
因此答案生成于

A ( x ) + A 2 ( x ) − B ( x ) 2 ! + A 3 ( x ) − 3 A ( x ) B ( x ) + 2 C ( x ) 3 ! A(x) + \dfrac{A^2(x) - B(x)}{2!} +\dfrac{A^3(x) - 3A(x)B(x)+2C(x)}{3!} A(x)+2!A2(x)−B(x)+3!A3(x)−3A(x)B(x)+2C(x)

这个式子可以用 F F T FFT FFT加速。

// (模板中 struct cp 里加一句) cp operator/(const double&x) { return cp(r / x, i / x); }
cp A[N], B[N], C[N];
inline void work() {
    int n; in(n);
    int M = INT_MIN; 
    forn(i,1,n) {
        int t; in(t);
        A[t].r ++, B[t << 1].r ++, C[t * 3].r ++;
        M = max(M, t * 3);
    }
    int siz = 1;
    while ((siz <<= 1) <= M << 1);
    FFT(A, siz), FFT(B, siz), FFT(C, siz);
    forn(i, 0, siz) {
        A[i] = A[i] + (A[i] * A[i] - B[i]) / 2 + (A [i] * A[i] * A[i] - cp(3., 0) * A[i] * B[i] + cp(2., 0) * C[i]) / 6.;
    }
    FFT(A, siz, -1);
    forn (i, 1, M << 1) {
        LL ans = LL(A[i].r + .5);
        if (ans) {
            out(i), putchar(32), out(ans), putchar(10);
        }
    }
}

BZOJ3513 [MUTC2013]idiots(反向思考 + 生成函数 + FFT)

要求合法情况，需要判断两个条件，不如找出不合法的情况，再用总数，也即 ( n 3 ) \binom{n}{3} (3n)减去之
两边之和为 x x x的方案数是： A 2 ( x ) A^2(x) A2(x)，也可以理解为一个卷积 ∑ k = i + j a i × a j \sum_{k = i+j}a_i\times a_j ∑k=i+jai×aj，其中有重复选了同一木棍的情况 B ( x ) B(x) B(x)，那么实际的两边等于 x 的方案数就是 C ( x ) = A 2 ( x ) − B ( x ) 2 实际的两边等于x的方案数就是C(x) = \dfrac{A^2(x) - B(x)}{2} 实际的两边等于x的方案数就是C(x)=2A2(x)−B(x)
现在要求两边之和 ≤ x 两边之和\le x 两边之和≤x的情况：实际上就是 ∑ i = 0 x C ( i ) \sum_{i = 0}^xC(i) ∑i=0xC(i)，处理一个前缀和即可。
最终答案就是 1 − A ( x ) C ( x ) ( n 3 ) 1-\dfrac{ A(x)C(x)}{\binom{n}{3}} 1−(3n)A(x)C(x)

#include <bits/stdc++.h>
#define IOS ios::sync_with_stdio(false), cin.tie(nullptr)
#define forn(i, s, e) for (int i = s; i <= int(e); ++i)
#define forr(i, e, s) for (int i = e; i >= int(s); --i)
using namespace std;
using LL = long long;
static constexpr int N = 1000010;

// (模板略)

LL f[N << 1], g[N];
cp a[N];
void solve() {
    double ans;
    LL M = LLONG_MIN, siz = 1, tot, n;
    in(n);
    forn(i, 1, n) {
        LL t; in(t);
        a[t].r ++;
        M = max(M, t);
        f[t << 1] --; // 对于 t + t 的情况，是多余的，需要去重
    }
    forr(i, M, 1) g[i] = g[i + 1] + LL(a[i].r);
    
    while ((siz <<= 1) <= M << 1);
    FFT(a, siz);
    forn(i, 0, siz) a[i] = a[i] * a[i];
    FFT(a, siz, -1);

    forn(i, 0, siz) {
        f[i] += (LL)(a[i].r + .5);
        f[i] >>= 1;
        // 考虑到 i = j + k, (j, k)cyc 于是 i = k + j 也成立, 应去重
    }
    ans = tot = (n) * (n - 1) / 2 * (n - 2) / 3;
    // 2 | ( n × (n - 1) ) 并且有 3 | ( n × (n - 1) × (n - 2) ) 避免溢出
    // (虽然 n ^ 3 < LLONG_MAX ), 但是很秀(bushi)
    
    forn (i, 0, M) ans -= f[i] * g[i];
    // 统计不合法的情况
    
    printf("%.7f\n", 1. * ans / tot);
    // .f 而不是 .lf 老生常谈了
    
    fill_n(f, siz << 1, 0), fill_n(g, siz, 0), fill_n(a, siz, 0);
    // 用多少清多少, 而不是 memset
}
inline void work() {
    int tt; in(tt); while (tt --) solve();
}