2013长沙网络赛G Goldbach FFT

       先贴题目连接:http://acm.zju.edu.cn/changsha/showProblem.do?problemId=28

       题意是给一个数n（1<=n<=80000），问最多用三个素数，通过+，*表示n共有多少种方法。例如：8=2*2*2=2+2*3=2+3+3=3+5，所以共4种。另外注意一下，相同的数不同的符号也记作不同的表示，例如4=2*2=2+2,这是两种不同的表示...

        首先筛素数，之后a,a*b,a*b*c,a*b+c都可以暴力求出来，那么难处理的主要就是a+b和a+b+c这两种情况，这里用两次fft来计算，第一次求出k=a+b中每个k各有多少中表示方法，这里注意判一下重，删掉a+a的情况以及其他情况/2，第二次在求出k的基础上，求出l=k+c共有多少种表示方法，这里注意补上a+b+a的情况后，在处理掉a+b+c==c+a+b的情况，最后特判一下a+a+a的情况就是最后的答案了。另外8000内的素数一共也就不到8000个，所以这题直接暴力预处理也能做..这里拿来做fft的练习了...

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
using namespace std;
typedef long long ll;
const int maxn=80080;
const int mod=1000000007;
bool ok[maxn];
ll prime[maxn];
int cnt=0;

const double PI = acos(-1.0);
struct comp
{
    double r,i;
    comp(double rt=0,double it=0)
    {
        r=rt;
        i=it;
    }
    comp operator +(const comp& b)
    {
        return comp(r+b.r,i+b.i);
    }
    comp operator -(const comp &b)
    {
        return comp(r-b.r,i-b.i);
    }
    comp operator *(const comp &b)
    {
        return comp(r*b.r-i*b.i,r*b.i+i*b.r);
    }
};

void change(comp y[],int len)//二进制转置--雷德算法
{
    int i,j,k;
    for(i = 1, j = len/2;i < len-1;i++)
    {
        if(i < j)swap(y[i],y[j]);
        k = len/2;
        while( j >= k)
        {
            j -= k;
            k /= 2;
        }
        if(j < k)j += k;
    }
}

void fft(comp y[],int len,int on)
/* on=1 DFT 把一个多项式的系数向量转化为点集表示；
on=-1,IDFT 把一个点集转化成多项式的系数向量*/
{
    change(y,len);
    for(int h = 2;h <= len;h <<= 1)
    {
        comp wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
        for(int j = 0;j < len;j += h)
        {
            comp w(1,0);
            for(int k = j;k < j+h/2;k++)
            {
                comp u = y[k];
                comp t = w*y[k+h/2];
                y[k] = u+t;
                y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if(on == -1)
        for(int i = 0;i < len;i++)
            y[i].r /= len;
}

void conv(comp f[],int len)//求f的卷积
{
    fft(f,len,1);
    for (int i=0; i<len; i++)
    f[i]=f[i]*f[i];
    fft(f,len,-1);
}

void ss()
{
    memset(ok,true,sizeof ok);
    ok[2]=true;
    ok[1]=false;
    ok[0]=false;
    for (int i=2; i<maxn; i++)
    {
        if (ok[i])
        {
            for (int j=i+i; j<maxn; j+=i)
            if (ok[j]) ok[j]=false;
        }
    }
    cnt=0;
    for (int i=2; i<maxn; i++)
    if (ok[i])
    {
        prime[cnt++]=i;
    }
//    for (int i=0; i<cnt; i++)
//    cout<<prime[i]<<"\n";
//    cout<<endl;
//    cout<<cnt<<endl;
}
ll ans=0;
ll n;
ll tm[maxn<<2],num[maxn<<2],sum[maxn<<2];
comp x[maxn<<2],y[maxn<<2],z[maxn<<2];
int main()
{
//    freopen("in.txt","r",stdin);
//    freopen("out.txt","w",stdout);
    cnt=0;
    ss();
    memset(tm,0,sizeof tm);
    memset(num,0,sizeof num);
    memset(sum,0,sizeof sum);

    n=80000;
    int len=1;
    while(len<2*n) len<<=1;

    for (int i=0; i<n; i++)
    if (ok[i]) x[i]=comp(1,0);
    else x[i]=comp(0,0);

    for (int i=n+1; i<len; i++)
    x[i]=comp(0,0);

    fft(x,len,1);
    for (int i=0; i<len; i++)
    x[i]=x[i]*x[i];
    fft(x,len,-1);

    for (int i=0; i<len; i++)
    num[i]=(int)(x[i].r+0.5);

    for (int i=0; i<=n; i++)
    if (ok[i]) num[i+i]--;
    for (int i=0; i<=n; i++)
    num[i]/=2;

    for (int i=0; i<=n; i++)
    x[i]=comp(num[i],0);

    for (int i=n+1; i<len; i++)
    x[i]=comp(0,0);

    for (int i=0; i<=n; i++)
    if (ok[i]) y[i]=comp(1,0);
    else y[i]=comp(0,0);
    for (int i=n+1; i<len; i++)
    y[i]=comp(0,0);

    fft(x,len,1);
    fft(y,len,1);
    for (int i=0; i<len; i++)
    x[i]=x[i]*y[i];
    fft(x,len,-1);
    for (int i=0; i<len; i++)
    sum[i]=(int)(x[i].r+0.5);


    while(~scanf("%d",&n))
    {
        ans=0;
        if (ok[n]) ans++;//a

        for (int i=0; i<cnt && prime[i]<=n; i++)
          for (int j=i; j<cnt&& prime[j]<n && prime[i]*prime[j]<=n; j++)
          if (prime[i]*prime[j]==n) ans++;  //a*b


        for (int i=0; i<cnt && prime[i]<=n; i++)
          for (int j=i; j<cnt&& prime[i]*prime[j]<=n; j++)
          {
              int tp=n-prime[i]*prime[j];
              if (ok[tp]) ans++;//a*b+c
          }

        for (int i=0; i<cnt && prime[i]<=n; i++)
          for (int j=i; j<cnt&& prime[i]*prime[j]<=n; j++)
           for (int k=j; k<cnt&& prime[i]*prime[j]*prime[k]<=n; k++)
           if (prime[i]*prime[j]*prime[k]==n) ans++;//a*b*c

        if (n%2==0 && ok[n/2]) ans++;
        ans+=num[n];//a+b

        ans%=mod;
        int tmp=0;
        for (int i=0; i<=n; i++)
        {
            if (ok[i])
            {
                int tp=n-i;
                if (tp%2==0 && ok[tp/2]) tmp++;
            }
        }
        ans=(ans+(sum[n]+2*tmp)/3)%mod;
        if (n%3==0 && ok[n/3]) ans++;
        ans%=mod;
        //a+b+c
        printf("%d\n",ans);
    }
    return 0;
}