【FFT】ZOJ 3856 Goldbach
通道：http://acm.zju.edu.cn/onlinejudge/showProBlem.do?proBlemCode=3856
题意：最多使用三个素数做加法和乘法（无括号），问得到X有多少种方案
思路：
6种情况：
1、A，O(1)判断
2、A+B, FFT处理
3、A*B，lgN*lgN判断
4、A*B*C，lgN*lgN*lgN判断
5、A*B+C，有A*B+C=N，可判断N-C是否在情况3出现过
6、A+B+C，利用之前预处理出来的A+B的种数[记住这里去掉A=B的情况, 并且考虑的是A, B的组合而不是排列(我们假设A < B], 将这个计数的多项式和素数的多项式相乘, 得到的结果当中, 计算了一次A, B, C中有两个相等的情况[A=C或B = C,也就是说有两个数相同的算了一次];三个都不同的算了三次, 分别是(A, B, C), (A, C, B), (B, C, A), [注意前面我们考虑的是组合, 所以(A, B, C)和(B, A, C)一样, 只算一个],所以对于三个都不同的情况要除3。
计算6这种方法，有个特殊的技巧，那就是(s[i]+2*tmp[i])/3,全不相同和全相同可以很好理解，对于其中有2个相同的，可以理解为（1+2*1）/3 也去重了
代码：
1 #include <cstdio> 2 #include <cmath> 3 #include <cstring> 4 #include <vector> 5 #include <algorithm> 6 7 using namespace std; 8 9 typedef long long ll; 10 11 const ll MAX_N = 80007; 12 const ll MAX_M = 100007; 13 const double PI = acos(-1.0); 14 15 struct Complex { 16 double r, i; 17 Complex(double _r, double _i) { 18 r = _r; 19 i = _i; 20 } 21 Complex operator + (const Complex &c) { 22 return Complex(c.r + r, c.i + i); 23 } 24 Complex operator - (const Complex &c) { 25 return Complex(r - c.r, i - c.i); 26 } 27 Complex operator * (const Complex &c) { 28 return Complex(c.r * r - c.i * i, c.r * i + c.i * r); 29 } 30 Complex operator / (const int &c) { 31 return Complex(r / c, i / c); 32 } 33 Complex(){} 34 }; 35 namespace FFT { 36 int rev(int id, int len) { 37 int ret = 0; 38 for(int i = 0; (1 << i) < len; ++i) { 39 ret <<= 1; 40 if(id & (1 << i)) ret |= 1; 41 } 42 return ret; 43 } 44 Complex A[MAX_M << 3]; 45 void FFT(Complex *a, int len, int DFT) { 46 for(int i = 0; i < len; ++i) A[rev(i, len)] = a[i]; 47 for(int s = 1; (1 << s) <= len; ++s) { 48 int m = (1 << s); 49 Complex wm = Complex(cos(PI * DFT * 2 / m), sin(PI * DFT * 2 / m)); 50 for(int k = 0; k < len; k += m) { 51 Complex w = Complex(1, 0); 52 for(int j = 0; j < (m >> 1); j++) { 53 Complex t = w * A[k + j + (m >> 1)]; 54 Complex u = A[k + j]; 55 A[k + j] = u + t; 56 A[k + j + (m >> 1)] = u - t; 57 w = w * wm; 58 } 59 } 60 } 61 if(DFT == -1) for(int i = 0; i < len; ++i) A[i] = A[i] / len; 62 for(int i = 0; i < len; i++) a[i] = A[i]; 63 } 64 }; 65 66 bool p[MAX_N]; 67 ll prime[MAX_N], primeCnt; 68 Complex A[MAX_M << 3], B[MAX_M << 3]; 69 ll s2[MAX_N], s3[MAX_N], s4[MAX_N], s5[MAX_N], tmp[MAX_N]; 70 /* 71 s2 : a * b 72 s3 : a * b * c 73 s4 : a + b 74 s5 : a + b + c 75 s6 : a * b + c 76 tmp: 至少2个数相同 77 */ 78 void init() { 79 memset(p, 1, sizeof p); 80 p[0] = p[1] = false; primeCnt = 0; 81 for (int i = 2; i < MAX_N; ++i) { 82 if (p[i]) { 83 prime[primeCnt++] = i; 84 for (int j = i * 2; j < MAX_N; j += i) 85 p[j] = false; 86 } 87 } 88 for (int i = 0; i < primeCnt && prime[i] * prime[i] < MAX_N; ++i) 89 for (int j = i; j < primeCnt && prime[i] * prime[j] < MAX_N; ++j) { 90 ++s2[prime[i] * prime[j]]; 91 for (int k = j; k < primeCnt && prime[i] * prime[j] * prime[k] < MAX_N; ++k) 92 ++s3[prime[i] * prime[j] * prime[k]]; 93 } 94 int len = 1; 95 while (len <= MAX_N) len <<= 1; len <<= 1; 96 for (int i = 0; i < MAX_N; ++i) if (p[i]) A[i] = Complex(1, 0); 97 else A[i] = Complex(0, 0); 98 for (int i = MAX_N; i < len; ++i) A[i] = Complex(0, 0); 99 FFT::FFT(A, len, 1); 100 for (int i = 0; i < len; ++i) A[i] = A[i] * A[i]; 101 FFT::FFT(A, len, -1); 102 for (int i = 0; i < MAX_N; ++i) s4[i] = (ll) (A[i].r + 0.5); 103 for (int i = 0; i < primeCnt && prime[i] * 2 < MAX_N; ++i) --s4[prime[i] * 2]; 104 for (int i = 0; i < MAX_N; ++i) s4[i] >>= 1; 105 for (int i = 0; i < MAX_N; ++i) A[i] = Complex(s4[i], 0); 106 for (int i = MAX_N; i < len; ++i) A[i] = Complex(0, 0); 107 for (int i = 0; i < MAX_N; ++i) if (p[i]) B[i] = Complex(1, 0); 108 else B[i] = Complex(0, 0); 109 for (int i = MAX_N; i < len; ++i) B[i] = Complex(0, 0); 110 FFT::FFT(A, len, 1); FFT::FFT(B, len, 1); 111 for (int i = 0; i < len; ++i) A[i] = A[i] * B[i]; 112 FFT::FFT(A, len, -1); 113 for (int i = 0; i < primeCnt; ++i) { 114 for (int j = 0; j < primeCnt && prime[i] * 2 + prime[j] < MAX_N; ++j) { 115 ++tmp[prime[i] * 2 + prime[j]]; 116 } 117 } 118 for (int i = 0; i < MAX_N; ++i) { 119 s5[i] = (ll) (A[i].r + 0.5); 120 s5[i] = (s5[i] + tmp[i] * 2) / 3; 121 } 122 for (int i = 0; i < primeCnt && prime[i] * 2 < MAX_N; ++i) { 123 ++s4[prime[i] * 2]; 124 if (prime[i] * 3 < MAX_N) ++s5[prime[i] * 3]; 125 } 126 } 127 128 int n; 129 130 int main() { 131 init(); 132 while (1 == scanf("%d", &n)) { 133 ll ans = 0; if (p[n]) ++ans; 134 ans = ans + s2[n] + s3[n] + s4[n] + s5[n]; 135 // cal s6 136 for (int i = 0; i < primeCnt && prime[i] <= n; ++i) 137 ans += s2[n - prime[i]]; 138 printf("%lld\n", ans); 139 } 140 return 0; 141 }
