【ZOJ 3856】Goldbach（FFT）

题意：最多使用三个质数做加法和乘法，问得到X有多少种方案。

这道题的话分几种情况讨论清楚就好了。

1、只有一个数的情况，如果这个数是质数则方案数加1

2、有两个数的情况，可以将两个数做加法和乘法，加法的话将质数得到的数字做一遍FFT，乘法直接做sqrt(n) * sqrt(n)就好了

3、有三个数的情况，分别有三种方法，分别是a + b + c, a + b * c,a * b * c

对于a * b * c，暴力遍历一遍就好，对于a + b * c，首先暴力遍历下b * c，然后再将一个数的两个数的做一遍fft就好，唯一比较麻烦的是a + b + c，这里又可以分为三种情况

（1）三个数都不相同，可以先算出x = a + b的方案数，然后再直接和c做FFT，再将答案除以3就得到方案数

（2）其中三个数相同，那么对于询问n，我们可以得到若n / 3 == 0 && n / 3为质数，则方案数+1

（3）对于只有两个数相同的方案数，那么我们可以先枚举一下。若所求数为9，则两个数相同的方案数有2 + 2 + 5，但是若直接按照（1）中的计数方法进行计算的话，就会出现(2 + 2) + 5,(2 + 5) + 2这两种方案（事实上只算一种），而且对于情况（2），又变得无法区分了。因此我的解决方案是，将x = a + b 中的a == b的情况先去除，这样就保证了计算结果中不会出现（2）即a == b && b == c，然后最后在统计是否满足（2）。并且因为在去除除了a==
b之后，只会出现三个数互不相等和只有一对数相等的情况，则对于询问n，可以求出n = 2 * a + c的个数tt，然后最后的方案数就是(cnt3[n] + tt * 2) / 3。

具体过程可见代码：

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
const int N = 80005;
const double pi = acos(-1.0);
const int mod = 1e9 + 7;
typedef long long LL;
struct Complex{
    double r,i;
    Complex(double r = 0.0,double i = 0.0):r(r),i(i){}
    Complex operator + (const Complex &a)const{
        return Complex(r + a.r,i + a.i);
    }
    Complex operator - (const Complex &a)const{
        return Complex(r - a.r,i - a.i);
    }
    Complex operator * (const Complex &a)const{
        return Complex(r * a.r - i * a.i,r * a.i + i * a.r);
    }
    Complex operator * (const double &a)const{
        return Complex(r * a,i * a);
    }
};
int tn = ceil(log(8 * 1e4) / log(2.0)) + 1;
int prime[N],pcnt;
bool isprime[N];
LL res[N],cnt3[N];
LL cnt[N],cnt2[N],cnt11[N];
Complex x1[N << 2],x2[N << 2],tmp[N << 2];
int rev(int x){
    int res = 0;
    for(int i = 0 ; i < tn ; i ++){
        if(x & 1) res += 1 << tn - 1 - i;
        x >>= 1;
    }
    return res;
}
void fft(Complex A[],int n,int op){
    for(int i = 0 ; i < n ; i ++) tmp[ rev(i) ] = A[i];
    for(int i = 0 ; i < n ; i ++) A[i] = tmp[i];
    for(int i = 1 ; (1 << i) <= n ; i ++){
        int m = 1 << i;
        Complex wn(cos(op * 2 * pi / m),sin(op * 2 * pi / m));
        for(int k = 0 ; k < n ; k += m){
            Complex w(1,0),u,t;
            for(int j = 0 ; j < m / 2 ; j ++){
                u = A[k + j];
                t = w * A[k + j + m / 2];
                A[k + j] = u + t;
                A[k + j + m / 2] = u - t;
                w = w * wn;
            }
        }
    }
    if(op == -1)
        for(int i = 0 ; i < n ; i ++)
            A[i] = A[i] * (1.0 / n);
}
void getprime(){
    isprime[0] = isprime[1] = 1;
    for(int i = 4 ; i < N ; i += 2) isprime[i] = 1;
    for(int i = 3 ; i < N ; i += 2){
        if(isprime[i] == 0)
            for(int j = i + i ; j < N ; j += i)
                isprime[j] = 1;
    }
    pcnt = 0;
    for(int i = 2 ; i < N ; i ++) if(isprime[i] == 0) prime[pcnt ++] = i;
}
int len = 1 << tn;
void init(){
    for(int i = 0 ; i < pcnt ; i ++) cnt[ prime[i] ] ++;
    for(int i = 0 ; prime[i] * prime[i] < N ; i ++)
    for(int j = i ; j < pcnt && prime[i] * prime[j] < N ; j ++)
        cnt2[ prime[i] * prime[j] ] ++;
}
void do1(){
    for(int i = 0 ; i < pcnt ; i ++) res[ prime[i] ] ++;
}
void do2(){
    for(int i = 0 ; i < N ; i ++) res[i] += cnt2[i];//*
    for(int i = 0 ; i < N ; i ++) x1[i] = Complex(cnt[i],0);
    for(int i = N ; i < len ; i ++) x1[i] = Complex(0,0);
    fft(x1,len,1);
    for(int i = 0 ; i < len ; i ++) x1[i] = x1[i] * x1[i];
    fft(x1,len,-1);
    for(int i = 0 ; i < N ; i ++){
        if(i % 2 == 0 && !isprime[i / 2])
            cnt11[i] = (LL(x1[i].r + 0.5) + 1) / 2;
        else
            cnt11[i] = LL(x1[i].r + 0.5) / 2;
        res[i] = (res[i] + cnt11[i]) % mod;
    }//+
}
void do3(){
    for(int i = 0 ; prime[i] * prime[i] * prime[i] < N ; i ++)
    for(int j = i ; prime[i] * prime[j] * prime[j] < N ; j ++)
    for(int k = j ; k < pcnt && prime[i] * prime[j] * prime[k] < N ; k ++){
        res[ prime[i] * prime[j] * prime[k] ] ++;
    }
    for(int i = 0 ; i < N ; i ++) if(res[i] >= mod) res[i] -= mod;
    //* and *
    for(int i = 0 ; i < N ; i ++) x1[i] = Complex(cnt[i],0);
    for(int i = N ; i < len ; i ++) x1[i] = Complex(0,0);
    for(int i = 0 ; i < N ; i ++) x2[i] = Complex(cnt2[i],0);
    for(int i = N ; i < len ; i ++) x2[i] = Complex(0,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 < N ; i ++) res[i] = (res[i] + LL(x1[i].r + 0.5)) % mod;
    // + and *
    for(int i = 0 ; i < N ; i ++) x1[i] = Complex(cnt11[i] - (i % 2 == 0 && !isprime[i / 2]),0);
    for(int i = N ; i < len ; i ++) x1[i] = Complex(0,0);
    for(int i = 0 ; i < N ; i ++) x2[i] = Complex(cnt[i],0);
    for(int i = N ; i < len ; i ++) x2[i] = Complex(0,0);
//    for(int i = 0 ; i < 10 ; i ++) printf("%I64d ",i);printf("\n");
//    for(int i = 0 ; i < 10 ; i ++) printf("%.3f ",x1[i].r);printf("\n");
//    for(int i = 0 ; i < 10 ; i ++) printf("%.3f ",x2[i].r);printf("\n");
    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 < N ; i ++) cnt3[i] = LL(x1[i].r + 0.5);
//    for(int i = 0 ; i < 10 ; i ++) printf("%I64d ",cnt3[i]);printf("\n");
    //+ and +
}
int main()
{
    getprime();
    init();
    do1();
    do2();
    do3();
    int n;
    while(~scanf("%d",&n)){
        LL ans = res[n],tt = 0;
        for(int i = 0 ; i < pcnt && 2 * prime[i] <= n ; i ++){
            if(!isprime[ n - 2 * prime[i] ] && n != prime[i] * 3) tt ++;
        }
        if(n % 3 == 0 && !isprime[n / 3]) ans ++;
        ans = (ans + (cnt3[n] + 2 * tt) / 3) % mod;
        printf("%lld\n",ans);
    }
    return 0;
}

版权声明：本文为博主原创文章，未经博主允许不得转载。