[UVA 12298 Super Poker II]FFT/NTT

[UVA 12298 Super Poker II]FFT/NTT

分类: Math FFT/NTT

1. 题目链接

[UVA 12298 Super Poker II]

2. 题意描述

有一副扑克,有无数张牌。对于每个正合数p,恰有4张牌,黑桃p, 红桃p, 方块p。
现在告诉你,有些牌已经丢失了。
给定两个数 a,b , 求区间 [a,b] 中的每个整数 xi , 从4种花色中各选一张牌,问有多少种组合可以使得点数之和等于 xi

3. 解题思路

水题。很显然的fft。每个花色的牌对应一个多项式。那么最终的答案就是四个多项式相乘。
不过,要知道,结果最大可能达到1e14。 普通的fft会有精度误差。这时候可以用long double来提高精度。 或者用ntt。ntt中间运算过程中可能出现爆long long, 需要用到 O(1) 或者 O(logn) 快速乘啊。

还学到一个奇怪的姿势。可以一次对多个多项式一起进行fft变换,最后就只需要一次fft反变换就好了。23333…..

4. 实现代码

/** FFT**/
#include <queue>
#include <stack>
#include <ctime>
#include <cmath>
#include <cctype>
#include <cstdio>
#include <string>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <algorithm>

using namespace std;

typedef long long LL;
typedef long double LB;
typedef unsigned long long ULL;
typedef pair<int, int> PII;
typedef pair<LL, LL> PLL;
typedef vector<int> VI;

const int INF = 0x3f3f3f3f;
const LL INFL = 0x3f3f3f3f3f3f3f3fLL;
const long double PI = acos(-1.0);
void debug() { cout << endl; }
template<typename T, typename ...R> void debug (T f, R ...r) { cout << "[" << f << "]"; debug (r...); }

const int MAXN = 262144 + 5; /// 数组大小应为 2^k
//typedef complex<double> CP;
struct CP {
    long double x, y;
    CP() {}
    CP(long double x, long double y) : x(x), y(y) {}
    inline long double real() { return x; }
    inline CP operator * (const CP& r) const {
        return CP(x * r.x - y * r.y, x * r.y + y * r.x);
    }
    inline CP operator - (const CP& r) const {
        return CP(x - r.x, y - r.y);
    }
    inline CP operator + (const CP& r) const {
        return CP(x + r.x, y + r.y);
    }
    inline CP conj(const CP &r) { return CP(r.x, -r.y); }
};
CP a[MAXN], b[MAXN];
int r[MAXN];
void fft_init(int nm, int k) {
    for(int i = 0; i < nm; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (k - 1));
}
void fft(CP ax[], int nm, int op) {
    for(int i = 0; i < nm; ++i) if(i < r[i]) swap(ax[i], ax[r[i]]);
    for(int h = 2, m = 1; h <= nm; h <<= 1, m <<= 1) { // 枚举长度
        CP wn = CP(cos(op * 2 * PI / h), sin(op * 2 * PI / h));
        for(int i = 0; i < nm; i += h) { // 枚举所有长度为 h 的区间
            CP w(1, 0); // 旋转因子
            for(int j = i; j < i + m; ++j, w = w * wn) { // 枚举角度
                CP t = w * ax[j + m]; // 蝴蝶操作
                ax[j + m] = ax[j] - t;
                ax[j] = ax[j] + t;
            }
        }
    }
    if(op == -1) for(int i = 0; i < nm; ++i) ax[i].x /= nm;
}
void trans(LL ax[], LL bx[], int n, int m) {
    int nm = 1, k = 0;
    while(nm < 2 * n || nm < 2 * m) nm <<= 1, ++k;
    for(int i = 0; i < n; ++i) a[i] = CP(ax[i], 0);
    for(int i = 0; i < m; ++i) b[i] = CP(bx[i], 0);
    for(int i = n; i < nm; ++i) a[i] = CP(0, 0);
    for(int i = m; i < nm; ++i) b[i] = CP(0, 0);
    fft_init(nm, k);
    fft(a, nm, 1); fft(b, nm, 1);
    for(int i = 0; i < nm; ++i) a[i] = a[i] * b[i];
    fft(a, nm, -1);

    nm = n + m - 1;
    for(int i = 0; i < nm; ++i) ax[i] = (LL)(a[i].real() + 0.5);
}

int na, nb, nc;
LL s[MAXN], h[MAXN], c[MAXN], d[MAXN];
char op[100];
bool notprime[MAXN]; // fasle--prime true--compose
void init() {
    memset(notprime, false, sizeof(notprime));
    notprime[0] = notprime[1] = false;
    for(int i = 2; i < MAXN; ++i) {
        if(notprime[i]) continue;
        for(int j = (i << 1); j < MAXN; j += i) {
            notprime[j] = true;
        }
    }
}

int main() {
#ifdef ___LOCAL_WONZY___
    freopen ("input.txt", "r", stdin);
#endif // ___LOCAL_WONZY___
    init();
    while(~scanf("%d %d %d", &na, &nb, &nc)) {
        if(na == 0 && nb == 0 && nc == 0) break;
        for(int i = 0; i < nb; ++i) {
            if(notprime[i]) {
                s[i] = 1;
                h[i] = 1;
                c[i] = 1;
                d[i] = 1;
            } else {
                s[i] = 0;
                h[i] = 0;
                c[i] = 0;
                d[i] = 0;
            }
        }

        while(nc --) {
            scanf("%s", op);
            int len = strlen(op), val;
            char tp = op[len - 1]; op[len - 1] = '\0';
            sscanf(op, "%d", &val);
            switch(tp) {
            case 'S': s[val] = 0; break;
            case 'H': h[val] = 0; break;
            case 'C': c[val] = 0; break;
            case 'D': d[val] = 0; break;
            }
        }

        trans(s, h, nb, nb);
        trans(c, d, nb, nb);
        trans(s, c, nb, nb);
        for(int i = na; i <= nb; ++i) {
            printf("%lld\n", s[i]);
        }
        puts("");
    }
#ifdef ___LOCAL_WONZY___
    cout << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC * 1000 << " ms." << endl;
#endif // ___LOCAL_WONZY___
    return 0;
}
/** NTT **/
#include <queue>
#include <stack>
#include <ctime>
#include <cmath>
#include <cctype>
#include <cstdio>
#include <string>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <algorithm>

using namespace std;

typedef long long LL;
typedef long double LB;
typedef unsigned long long ULL;
typedef pair<int, int> PII;
typedef pair<LL, LL> PLL;
typedef vector<int> VI;

const int INF = 0x3f3f3f3f;
const LL INFL = 0x3f3f3f3f3f3f3f3fLL;
const double PI = acos(-1.0);
void debug() { cout << endl; }
template<typename T, typename ...R> void debug (T f, R ...r) { cout << "[" << f << "]"; debug (r...); }

const int MAXN = 262144 + 5;

const LL G = 3, MOD = 40531930642382849LL; //(479<<21)+1; //998244353;
LL a[MAXN], b[MAXN];
int r[MAXN];
/** 将z = x * y % MOD 替换为 z = quick_mul(x, y)**/
LL quick_mul(LL x, LL y) {
    return (x * y - (LL)(x / (long double)MOD * y + 1e-3) * MOD + MOD) % MOD;
}
template<typename T> T quick_pow(T a, T b) {
    T ret = 1;
    while(b) {
        if(b & 1) ret = quick_mul(ret, a);
        a = quick_mul(a, a);
        b >>= 1;
    }
    return ret;
}
void ntt_init(int nm, int k) {
    for(int i = 0; i < nm; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (k - 1));
}
template<typename T> void ntt(T ax[], int nm, int op) {
    for(int i = 0; i < nm; ++i) if(i < r[i]) swap(ax[i], ax[r[i]]);
    for(int h = 2, m = 1; h <= nm; h <<= 1, m <<= 1) {
        T wn = quick_pow(G, (MOD - 1) / h);
        for(int i = 0; i < nm; i += h) {
            T w = 1;
            for(int j = i; j < i + m; ++j, w = quick_mul(w, wn)) {
                T t = quick_mul(w, ax[j + m]);
                ax[j + m] = ax[j] - t + MOD;
                if(ax[j + m] >= MOD) ax[j + m] -= MOD;
                ax[j] = ax[j] + t;
                if(ax[j] >= MOD) ax[j] -= MOD;
            }
        }
    }
    if(op == -1) {
        for(int i = 1; i < nm / 2; i++) swap(ax[i], ax[nm - i]); // Caution Here!
        T inv = quick_pow((LL)nm, MOD - 2);
        for(int i = 0; i < nm; ++i) ax[i] = quick_mul(ax[i], inv);
    }
}
template<typename T> void trans(T ax[], T bx[], int n, int m) {
    int nm = 1, k = 0;
    while(nm < 2 * n || nm < 2 * m) nm <<= 1, ++k;
    for(int i = 0; i < n; ++i) a[i] = ax[i];
    for(int i = 0; i < m; ++i) b[i] = bx[i];
    for(int i = n; i < nm; ++i) a[i] = 0;
    for(int i = m; i < nm; ++i) b[i] = 0;
    ntt_init(nm, k);
    ntt(a, nm, 1); ntt(b, nm, 1);
    for(int i = 0; i < nm; ++i) a[i] = quick_mul(a[i], b[i]);
    ntt(a, nm, -1);

    nm = n + m - 1;
    for(int i = 0; i < nm; ++i) ax[i] = a[i];
}


int na, nb, nc;
LL s[MAXN], h[MAXN], c[MAXN], d[MAXN];
char op[100];
bool notprime[MAXN]; // fasle--prime true--compose
void init() {
    memset(notprime, false, sizeof(notprime));
    notprime[0] = notprime[1] = false;
    for(int i = 2; i < MAXN; ++i) {
        if(notprime[i]) continue;
        for(int j = (i << 1); j < MAXN; j += i) {
            notprime[j] = true;
        }
    }
}

int main() {
#ifdef ___LOCAL_WONZY___
    freopen ("input.txt", "r", stdin);
#endif // ___LOCAL_WONZY___
    init();
    while(~scanf("%d %d %d", &na, &nb, &nc)) {
        if(na == 0 && nb == 0 && nc == 0) break;
        for(int i = 0; i < nb; ++i) {
            if(notprime[i]) {
                s[i] = 1;
                h[i] = 1;
                c[i] = 1;
                d[i] = 1;
            } else {
                s[i] = 0;
                h[i] = 0;
                c[i] = 0;
                d[i] = 0;
            }
        }

        while(nc --) {
            scanf("%s", op);
            int len = strlen(op), val;
            char tp = op[len - 1]; op[len - 1] = '\0';
            sscanf(op, "%d", &val);
            switch(tp) {
            case 'S': s[val] = 0; break;
            case 'H': h[val] = 0; break;
            case 'C': c[val] = 0; break;
            case 'D': d[val] = 0; break;
            }
        }
        trans(s, h, nb, nb);
        trans(c, d, nb, nb);
        trans(s, c, nb, nb);
        for(int i = na; i <= nb; ++i) {
            printf("%lld\n", s[i]);
        }
        puts("");
    }
#ifdef ___LOCAL_WONZY___
    cout << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC * 1000 << " ms." << endl;
#endif // ___LOCAL_WONZY___
    return 0;
}
#include <queue>
#include <stack>
#include <ctime>
#include <cmath>
#include <cctype>
#include <cstdio>
#include <string>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <algorithm>

using namespace std;

typedef long long LL;
typedef long double LB;
typedef unsigned long long ULL;
typedef pair<int, int> PII;
typedef pair<LL, LL> PLL;
typedef vector<int> VI;

const int INF = 0x3f3f3f3f;
const LL INFL = 0x3f3f3f3f3f3f3f3fLL;
const long double PI = acos(-1.0);
void debug() { cout << endl; }
template<typename T, typename ...R> void debug (T f, R ...r) { cout << "[" << f << "]"; debug (r...); }

const int MAXN = 262144 + 5; /// 数组大小应为 2^k
//typedef complex<double> CP;
struct CP {
    long double x, y;
    CP() {}
    CP(long double x, long double y) : x(x), y(y) {}
    inline long double real() { return x; }
    inline CP operator * (const CP& r) const {
        return CP(x * r.x - y * r.y, x * r.y + y * r.x);
    }
    inline CP operator - (const CP& r) const {
        return CP(x - r.x, y - r.y);
    }
    inline CP operator + (const CP& r) const {
        return CP(x + r.x, y + r.y);
    }
    inline CP conj(const CP &r) { return CP(r.x, -r.y); }
};
CP a[MAXN], b[MAXN], c[MAXN], d[MAXN];
int r[MAXN];
void fft_init(int nm, int k) {
    for(int i = 0; i < nm; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (k - 1));
}
void fft(CP ax[], int nm, int op) {
    for(int i = 0; i < nm; ++i) if(i < r[i]) swap(ax[i], ax[r[i]]);
    for(int h = 2, m = 1; h <= nm; h <<= 1, m <<= 1) { // 枚举长度
        CP wn = CP(cos(op * 2 * PI / h), sin(op * 2 * PI / h));
        for(int i = 0; i < nm; i += h) { // 枚举所有长度为 h 的区间
            CP w(1, 0); // 旋转因子
            for(int j = i; j < i + m; ++j, w = w * wn) { // 枚举角度
                CP t = w * ax[j + m]; // 蝴蝶操作
                ax[j + m] = ax[j] - t;
                ax[j] = ax[j] + t;
            }
        }
    }
    if(op == -1) for(int i = 0; i < nm; ++i) ax[i].x /= nm;
}
void trans(LL ax[], LL bx[], LL cx[], LL dx[], int n, int m) {
    int nm = 1, k = 0;
    while(nm < 4 * n || nm < 4 * m) nm <<= 1, ++k; // 需要注意最高次幂
    for(int i = 0; i < n; ++i) a[i] = CP(ax[i], 0);
    for(int i = 0; i < m; ++i) b[i] = CP(bx[i], 0);
    for(int i = 0; i < n; ++i) c[i] = CP(cx[i], 0);
    for(int i = 0; i < m; ++i) d[i] = CP(dx[i], 0);
    for(int i = n; i < nm; ++i) a[i] = CP(0, 0);
    for(int i = m; i < nm; ++i) b[i] = CP(0, 0);
    for(int i = n; i < nm; ++i) c[i] = CP(0, 0);
    for(int i = m; i < nm; ++i) d[i] = CP(0, 0);
    fft_init(nm, k);
    fft(a, nm, 1); fft(b, nm, 1);
    fft(c, nm, 1); fft(d, nm, 1);
    for(int i = 0; i < nm; ++i) a[i] = a[i] * b[i] * c[i] * d[i];
    fft(a, nm, -1);

    nm = n + m - 1;
    for(int i = 0; i < nm; ++i) ax[i] = (LL)(a[i].real() + 0.5);
}

int na, nb, nc;
LL S[MAXN], H[MAXN], C[MAXN], D[MAXN];
char op[100];
bool notprime[MAXN]; // fasle--prime true--compose
void init() {
    memset(notprime, false, sizeof(notprime));
    notprime[0] = notprime[1] = false;
    for(int i = 2; i < MAXN; ++i) {
        if(notprime[i]) continue;
        for(int j = (i << 1); j < MAXN; j += i) {
            notprime[j] = true;
        }
    }
}

int main() {
#ifdef ___LOCAL_WONZY___
    freopen ("input.txt", "r", stdin);
#endif // ___LOCAL_WONZY___
    init();
    while(~scanf("%d %d %d", &na, &nb, &nc)) {
        if(na == 0 && nb == 0 && nc == 0) break;
        for(int i = 0; i < nb; ++i) {
            if(notprime[i]) {
                S[i] = 1;
                H[i] = 1;
                C[i] = 1;
                D[i] = 1;
            } else {
                S[i] = 0;
                H[i] = 0;
                C[i] = 0;
                D[i] = 0;
            }
        }

        while(nc --) {
            scanf("%s", op);
            int len = strlen(op), val;
            char tp = op[len - 1]; op[len - 1] = '\0';
            sscanf(op, "%d", &val);
            switch(tp) {
            case 'S': S[val] = 0; break;
            case 'H': H[val] = 0; break;
            case 'C': C[val] = 0; break;
            case 'D': D[val] = 0; break;
            }
        }

        trans(S, H, C, D, nb, nb);
        for(int i = na; i <= nb; ++i) {
            printf("%lld\n", S[i]);
        }
        puts("");
    }
#ifdef ___LOCAL_WONZY___
    cout << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC * 1000 << " ms." << endl;
#endif // ___LOCAL_WONZY___
    return 0;
}
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值