ZOJ3856 Goldbach

转载请注明出处，谢谢http://blog.csdn.net/bigtiao097?viewmode=contents
参考文章http://blog.csdn.net/gatevin/article/details/46917303

题意：

给一个数X, 1 < X <= 80000, 现在X是由至多三个质数通过加法和乘法组成的, 问有多少种这样的质数表达式
可能的情况有以下几种:
①p1
②p1 + p2
③p1 + p2 + p3
④p1*p2
⑤p1*p2*p3
⑥p1*p2 + p3
例如8 = 2 * 2 * 2 = 2 + 2 * 3 = 2 + 3 + 3 = 3 + 5, so the answer is 4

思路：

①p1：简单素筛以下就好了
②p1 + p2 ：根据①构造两个多项式，相乘（用FFT），要注意2+3和3+2是同一种情况
③p1 + p2 + p3：这个稍微麻烦一点，先拿②得到的多项式与①得到的多项式相乘一下，然后还是要去重的，不过计算要等到读入X根据具体值来算比较好，去掉2+3+3与3+3+2一样的情况，还要特判一下2+2+2的情况（这里只是举一个例子，并不一定只是2、3）
④p1*p2：80000以内的素数大概不超过10000个，枚举一下就好了,这个不用预处理也行，读到X后O(n)即可判断，这个表达是唯一的，但是由于后面⑥的预处理用到了，所以就先预处理一下
⑤p1*p2*p3：这个也是读到X后再枚举，也是唯一的，用到了⑤的结果
⑥p1*p2 + p3：拿④的结果构造多项式，和①的多项式相乘(FFT)即可

这个分类还是有点麻烦的，不小心就会漏掉一些情况，不懂的可以看一下下面的代码，里面有注释

---

Result: Accepted Mem: 16.3MB Time: 410ms
具体代码如下：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 1e9+7;
const double PI = acos(-1);
const int maxn = 8e4+5;//maxn = max(len1,len2)
const int maxm = 4*maxn;//k是>=maxn的最小的2的幂,maxm=2*k
int len;//len1，len2分别为两个多项式的最高次数+1
int num[maxm];
int a[maxn];
int b[maxm],c[maxm];
int d[maxn];
int f[maxm];//复数结构体
int prime[maxn];//第i个素数
bool vis[maxn];
int p;//素数个数
int x;
int t;
ll ans;
void sieve()
{
    p=0;
    int m = sqrt(maxn+0.5);
    memset(vis, 0, sizeof vis);
    for(int i = 2; i <= m; i++) if(!vis[i])
        for(int j = i*i; j < maxn; j+=i) vis[j] = 1;
    for(int i=2;i<maxn;i++) if(!vis[i]) prime[++p]= i;
}
struct Complex
{
    double x, y;//实部为x，虚部为y
    Complex(double x=0, double y=0):x(x),y(y){}
    Complex operator+(const Complex &rhs) const { return Complex(x+rhs.x, y+rhs.y);}
    Complex operator-(const Complex &rhs) const { return Complex(x-rhs.x, y-rhs.y);}
    Complex operator*(const Complex &rhs) const { return Complex(x*rhs.x-y*rhs.y,x*rhs.y+y*rhs.x);}
};
Complex x1[maxm],x2[maxm];
/*
进行FFT和IFFT前的反转变换。
将位置i和(i二进制反转后位置)互换。
len必须取2的幂
*/
void change(Complex *x, int len)
{
    Complex t;
    for(int i = 1, j = len/2; i < len-1; i++)
    {
        //交换下标互反的元素，i<j保证交换一次
        //i做正常的+1，j做反转的+1，始终保持i和j是反转的
        if(i < j)
        {
            t = x[i];
            x[i] = x[j];
            x[j] = t;
        }
        int k = len / 2;
        while(j >= k)
        {
            j -= k;
            k >>= 1;
        }
        if(j < k) j += k;
    }
}
/*
做FFT
len必须为2的幂
on==1时是DFT，on==-1时是IDFT
*/
void fft(Complex *x, int len, int on)
{
    change(x, len);
    for(int h = 2; h <= len; h <<= 1)
    {
        Complex wn(cos(-on*2*PI/h), sin(-on*2*PI/h));
        for(int j = 0; j < len; j += h)
        {
            Complex w(1, 0);
            for(int k = j; k < j+h/2; k++)
            {
                Complex u = x[k];
                Complex t = w*x[k+h/2];
                x[k] = u+t;
                x[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if(on == -1)
        for(int i = 0; i < len; i++) x[i].x /= len;
}
void workFFT(int *xx,int len1,int *yy,int len2,int *num)
{
    len = 1;
    memset(x1,0,sizeof x1);
    memset(x2,0,sizeof x2);
    memset(num,0,sizeof num);
    while(len < len1*2 || len < len2*2) len <<=1;
    for(int i = 0; i < len; i++)
            x1[i] =  x2[i] = Complex(0,0);
    for(int i = 0; i < len1; i++) x1[i] = Complex(xx[i],0);
    for(int i = 0; i < len2; i++) x2[i] = Complex(yy[i],0);
    fft(x1,len,1);
    fft(x2,len,1);
    for(int i = 0; i < len; i++) x1[i] = x1[i]*x2[i];
    fft(x1,len,-1);
    for(int i = 0; i < len; i++) num[i] = (ll)(x1[i].x+0.5)%mod;
}
int main()
{
    sieve();
    //预处理p1的情况
    for(int i=1;i<=p;i++)
        a[prime[i]]=1;
    //预处理p1+p2的情况
    workFFT(a,maxn,a,maxn,b);
    for(int i=1;i<=p;i++)
        b[prime[i]+prime[i]]++;
    for(int i=1;i<maxn;i++)
        b[i]/=2;
    //预处理p1+p2+p3的情况(并没有完全处理完，等读到x再计算
    workFFT(b,maxn,a,maxn,c);
    //预处理p1*p2的情况
    for(int i=1;i<=p;i++)
        for(int j=i;j<=p;j++)
        {
            if(1LL*prime[i]*prime[j]>=maxn)break;
            d[prime[i]*prime[j]]++;
        }
    //预处理p1*p2+p3的情况
    workFFT(d,maxn,a,maxn,f);
    while(~scanf("%d",&x))
    {
        t = 0;
        ans = a[x];
        ans = (ans+b[x])%mod;
        //继续处理p1+p2+p3
        for(int i=1;i<=p;i++)
        {
            if(x<2*prime[i])break;
            if(!vis[x-2*prime[i]]&&x!=3*prime[i])//有两个数相同且不是三个数都相同
                t++;
        }
        ans = (ans+(c[x]+t)/3)%mod;
        if(x%3==0&&!vis[x/3])ans++;
        ans = (ans+d[x])%mod;
        //处理p1*p2*p3的情况
        for(int i=2;i<x;i++)
            if(!vis[i]&&x%i==0&&d[x/i])
        {
            ans++;
            break;
        }
        ans = (ans+f[x])%mod;
        cout<<ans<<endl;
    }
}