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*kint len;//len1，len2分别为两个多项式的最高次数+1int 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;    }}/*做FFTlen必须为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;    }}