题目大意:
就是现在给出一个数X, 1 < X <= 80000, 现在X是由至多三个质数通过加法和乘法组成的, 问有多少种这样的质数表达式
可用的表达式形式:
p1
p1 + p2
p1 + p2 + p3
p1*p2
p1*p2*p3
p1*p2 + p3
例如8可以是 3 + 5, 2 + 3 + 3, 2*3 + 2, 2*2*2共四种
大致思路:
感觉还是一个比较麻烦的计数题
对于p1型直接质数筛找出80000以内的质数即可, 大概有7000~8000个质数
对于p1 + p2型可以判断X是否是一个质数的两倍, 也是预处理之后O(1)判断, 然后预处理出两个不同质数的乘积, 这个用FFT预处理出来, 也很容易
对于p1*p2型, 枚举筛出来的质数乘积, 控制不要大于80000找出不超过80000的所有两个质数乘积的数, 这个表达也是唯一的, 存在就+1结束
对于p1*p2*p3型, 枚举质数p3是否是X的约数, 是就判断X/p3是否是p1*p2型即可, 当然这个表达式最多一个, 找到就结果+1结束
对于p1*p2 + p3型首先根据处理出来的80000以内p1*p2的数构造多项式, 与质数那个的多项式相乘, 也是FFT预处理即可
最麻烦的应该就是p1 + p2 + p3型, 因为要去重, 需要分类讨论一下
首先p1 = p2 = p3的情况最简单, 判断X % 3 == 0 && X/3是质数即可
如果p1, p2, p3中有两个相同呢, 这个情况可以枚举不同的那个质数p3, 当(X - p3) / 2也是质数的时候这种情形+1
对于p1, p2, p3都不相同的情形, 首先利用之前预处理出来的p1 + p2的种数, (记住这里去掉p1 = p2的情况, 并且考虑的是p1, p2的组合而不是排列(我们假设p1 < p2), 将这个计数的多项式和质数的多项式相乘, 得到的结果当中, 计算了一次p1, p2, p3中有两个相等的情况, 那么要么p1 = p3, 要么p2 = p3也就是说有两个数相同的算了一次, 三个都不同的算了三次, 分别是(p1, p2, p3), (p1, p3, p2), (p2, p3, p1), (注意前面我们考虑的是组合, 所以(p1, p2, p3)和(p2, p1, p3)一样, 只算一个), 那么我们只需要把这个结果减去之前统计的p1, p2, p3中有一个相同的结果之后除以3即可
讨论比较复杂, 需要细心一点, 还是不太懂的话可以看代码注释的解释
代码如下:
Result : Accepted Memory : 25084 KB Time : 250 ms
/*
* Author: Gatevin
* Created Time: 2015/7/16 16:32:02
* File Name: ZOJ3856.cpp
*/
#include<iostream>
#include<sstream>
#include<fstream>
#include<vector>
#include<list>
#include<deque>
#include<queue>
#include<stack>
#include<map>
#include<set>
#include<bitset>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cctype>
#include<cmath>
#include<ctime>
#include<iomanip>
using namespace std;
const double eps(1e-8);
typedef long long lint;
const lint mod = 1e9 + 7;
bool isPrime[80010];
bool check[80010];
vector<int> prime;
void getPrime()
{
memset(isPrime, 0, sizeof(isPrime));
memset(check, 0, sizeof(check));
isPrime[0] = isPrime[1] = 0;
check[0] = check[1] = 1;
for(int i = 2; i <= 80000; i++)
{
if(check[i]) continue;
isPrime[i] = 1; prime.push_back(i);
for(int j = i; j <= 80000; j += i)
check[j] = 1;
}
return;
}
const double PI = acos(-1.0);
struct Complex
{
double real, image;
Complex(double _real, double _image)
{
real = _real;
image = _image;
}
Complex(){}
};
Complex operator + (const Complex &c1, const Complex &c2)
{
return Complex(c1.real + c2.real, c1.image + c2.image);
}
Complex operator - (const Complex &c1, const Complex &c2)
{
return Complex(c1.real - c2.real, c1.image - c2.image);
}
Complex operator * (const Complex &c1, const Complex &c2)
{
return Complex(c1.real*c2.real - c1.image*c2.image, c1.real*c2.image + c1.image*c2.real);
}
int rev(int id, int len)
{
int ret = 0;
for(int i = 0; (1 << i) < len; i++)
{
ret <<= 1;
if(id & (1 << i)) ret |= 1;
}
return ret;
}
Complex A[1 << 18];
void FFT(Complex *a, int len, int DFT)
{
for(int i = 0; i < len; i++)
A[rev(i, len)] = a[i];
for(int s = 1; (1 << s) <= len; s++)
{
int m = (1 << s);
Complex wm = Complex(cos(DFT*2*PI/m), sin(DFT*2*PI/m));
for(int k = 0; k < len; k += m)
{
Complex w = Complex(1, 0);
for(int j = 0; j < (m >> 1); j++)
{
Complex t = w*A[k + j + (m >> 1)];
Complex u = A[k + j];
A[k + j] = u + t;
A[k + j + (m >> 1)] = u - t;
w = w*wm;
}
}
}
if(DFT == -1) for(int i = 0; i < len; i++) A[i].real /= len, A[i].image /= len;
for(int i = 0; i < len; i++) a[i] = A[i];
return;
}
Complex P[1 << 18];
Complex P1_P2[1 << 18];//p1+p2 && p1 != p2
Complex P1P2[1 << 18];//p1*p2
Complex P1_P2_P3[1 << 18];//p1 + p2 + p3, p1 != p2
Complex P1P2_P3[1 << 18];//p1*p2 + p3
bool isP1P2[80010];
int main()
{
getPrime();//预处理[1, 80000]之间的质数
//预处理p1 + p2且p1 != p2的方案数
for(int i = 0; i <= 80000; i++)
if(isPrime[i]) P[i] = Complex(1, 0);
else P[i] = Complex(0, 0);
int len = (1 << 18);
for(int i = 80001; i < len; i++) P[i] = Complex(0, 0);
FFT(P, len, 1);
for(int i = 0; i < len; i++)
P1_P2[i] = P[i]*P[i];
FFT(P1_P2, len, -1);
for(int i = 2; i <= 80000; i++)//质数来自同一组合
{
if((i & 1) == 0 && isPrime[i >> 1]) P1_P2[i].real -= 1;
P1_P2[i].real /= 2;//p1 + p2 && p1 != p2
}
for(int i = 80001; i < len; i++)//防止1 << 18不够用...
P1_P2[i].real = 0;
FFT(P1_P2, len, 1);
for(int i = 0; i < len; i++)
P1_P2_P3[i] = P1_P2[i]*P[i];
FFT(P1_P2_P3, len, -1); FFT(P1_P2, len, -1);//p1 + p2 + p3 && p1 != p2
//预处理出p1*p2
memset(isP1P2, 0, sizeof(isP1P2));
for(int i = 0, sz = prime.size(); i < sz; i++)
{
for(int j = i; j < sz; j++)
{
if(prime[i]*1. > 80000*1./prime[j]) break;
isP1P2[prime[i]*prime[j]] = 1;
}
}
//预处理p1*p2 + p3的方案个数
for(int i = 0; i <= 80000; i++)
if(isP1P2[i]) P1P2[i] = Complex(1, 0);
else P1P2[i] = Complex(0, 0);
for(int i = 80001; i < len; i++)
P1P2[i] = Complex(0, 0);
FFT(P1P2, len, 1);
for(int i = 0; i < len; i++)
P1P2_P3[i] = P[i]*P1P2[i];
FFT(P1P2_P3, len, -1);
int X;
while(~scanf("%d", &X))
{
lint ans = 0;
if(isPrime[X]) ans++;//p1
ans = (ans + (lint)(P1_P2[X].real + 0.5)) % mod;//p1 + p2 && p1 != p2
if(!(X & 1) && isPrime[X >> 1]) ans++;//p1 + p1
ans = (ans + (lint)(P1P2_P3[X].real + 0.5)) % mod;//p1*p2 + p3
if(X % 3 == 0 && isPrime[X / 3]) ans++;//p1 + p1 + p1
//计算p1 + p1 + p2形式的
int cnt = 0;
for(int i = 0, sz = prime.size(); i < sz && 2*prime[i] < X; i++)
if(isPrime[X - 2*prime[i]] && X != 3*prime[i])
cnt++;
ans = (ans + cnt) % mod;//p1 + p1 + p2
//得到P1_P2_P3(p1 + p2 + p3, p1 != p2, p1和p2作的组合算一次 * p3)
lint tmp = (lint)(P1_P2_P3[X].real + 0.5);//p1 + p2 + ?, , 这一部分将p1 + p1 + p2型算了一次, p1 + p2 + p3型算了3次
ans = (ans + (tmp - cnt) / 3) % mod;//p1 + p2 + p3
if(isP1P2[X]) ans++;//p1*p2
for(int i = 0, sz = prime.size(); i < sz && prime[i] < X; i++)//p1*p2*p3
if(X % prime[i] == 0 && isP1P2[X / prime[i]])
{
ans++;
break;
}
printf("%lld\n", ans % mod);
}
return 0;
}