ZOJ 3856 Goldbach FFT计数 2013年长沙网络赛

题目大意：

就是现在给出一个数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;
}