[UVA 12298 Super Poker II]FFT/NTT
分类: Math
FFT/NTT
1. 题目链接
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;
}