题目链接

  UVa12298 Super Joker II

题意

  有一副超级扑克，包含无数张牌。对于每个正合数p，恰好有4张牌：黑桃p，红桃p，梅花p和方块p（分别用pS、pH、pC 和pD 表示）。没有其他类型的牌。
  给定一个整数n，从4种花色中各选一张牌，问有多少种组合可以使得点数之和等于n。例如，n=24的时候，有一种组合方法是4S+6H+4C+10D，下图所示。

  不巧的是，有些牌已经丢失（题目会提供已经丢失的牌的列表）。并且为了让题目更有趣，我们还会提供两个正整数a和b，你的任务是按顺序输出n=a, n=a+1, n=a+2, …, n=b时的答案。

输入格式

  输入包含不超过25组数据。每组数据的第一行为3个整数a, b, c，其中c是已丢失的牌的张数。第二行包含c个不同的字符串，即已丢失的牌。这些牌形如pS, pH, pC 或者pD，其中p是一个正合数。输入结束标志为a=b=c=0。最多有一组数据满足a=1, b=50000且c≤10000，其他数据满足1≤a≤b≤100, 0≤c≤10。

输出格式

  对于每组数据，输出p 行，每行一个整数。每组数据后输出一个空行。

分析

  只需要对生成函数理解清楚了，就可以套用快速傅里叶变换（FFT）求解。初次接触FFT的话，推荐看看这篇知乎，GhostLX大佬给出了模板：

void dft(Complex a[])
{
    for (int i = 0; i < tot; i++)
    {
        if (i < rev[i])
            swap(a[i], a[rev[i]]); //只需要交换一次就行了，交换两次等于没有换
    }
    for (int mid = 1; mid < tot; mid <<= 1)
    {
        auto w1 = Complex({cos(PI / mid), sin(PI / mid)});
        for (int i = 0; i < tot; i += mid * 2)
        {
            auto wk = Complex({1, 0}); //初始为w(0,mid)
            for (int j = 0; j < mid; j++, wk = wk * w1) //单位根递推式
            {
                auto x = a[i + j], y = wk * a[i + j + mid];
                a[i + j] = x + y, a[i + j + mid] = x - y;
            }
        }
    }
}

  GhostLX大佬模板的单位根递推式像这样写wk = wk * w1，随着递推加深其实会带来精度误差的。推荐用这个模板：

namespace fft {
    #include <cmath>
    #define N 1<<21
    int res[N], tot;
    struct complex {
        double x, y;
        void operator+= (const complex &t) {
            x += t.x; y += t.y;
        }
        complex operator- (const complex &t) const {
            return {x - t.x, y - t.y};
        }
        complex operator* (const complex &t) const {
            return {x * t.x - y * t.y, x * t.y + y * t.x};
        }

    } a[N], b[N];
    
    /**
     * inv=1时求的是傅里叶变换（DFT），inv=-1时求的是傅里叶逆变换（IDFT）
     */
    void fft(complex (&a)[N], int inv) {
        for (int i=0, j=0; i<tot; ++i) { // 原地快速bit reversal
            if(j > i) {complex t = a[i]; a[i] = a[j]; a[j] = t;}
            int k = tot;
            while(j & (k >>= 1)) j &= ~k;
            j |= k;
        }
        for(int step=1; step<tot; step<<=1) {
            // 把每相邻两个“step点DFT”通过一系列蝴蝶变换合并为一个“2*step点DFT”
            double alpha = inv*M_PI / step;
            // 为求高效，我们并不是依次执行各个完整的DFT合并，而是枚举下标k
            // 对于一个下标k，执行所有DFT合并中该下标对应的蝴蝶变换，即通过E[k]和O[k]计算X[k]
            // 蝴蝶变换参考：http://en.wikipedia.org/wiki/Butterfly_diagram
            for(int k=0; k<step; k++) {
                // 计算omega^k
                complex wk = {cos(alpha*k), sin(alpha*k)};
                for(int Ek=k; Ek<tot; Ek += step<<1) { // Ek是某次DFT合并中E[k]在原始序列中的下标
                    int Ok = Ek + step; // Ok是该DFT合并中O[k]在原始序列中的下标
                    complex t = wk * a[Ok]; // 蝴蝶变换：x1 * omega^k
                    a[Ok] = a[Ek] - t;  // 蝴蝶变换：y1 = x0 - t
                    a[Ek] += t;         // 蝴蝶变换：y0 = x0 + t
                }
            }
        }
    }

    void workFFT(int n, int m) { // a[0, n], b[0, m]
        int bit = 0;
        while ((1 << bit) < n + m + 1) ++bit;
        tot = 1 << bit;
        fft(a, 1); fft(b, 1);
        for (int i=0; i<tot; ++i) a[i] = a[i] * b[i]; //点表示法直接运算
        fft(a, -1); //逆变换，点表示法转换为多项式表示法
        for (int i=0, j=m+n; i<=j; ++i) res[i]  = a[i].x / tot + 0.5; //向上取整
    }
}

AC 代码

#include <iostream>
#include <cmath>
using namespace std;

#define M 50020
#define N 1<<17
int na, nb, nc, tot; bool f[M] = {false};
struct complex {
    double x, y;
    void operator+= (const complex &t) {
        x += t.x; y += t.y;
    }
    complex operator- (const complex &t) const {
        return {x - t.x, y - t.y};
    }
    complex operator* (const complex &t) const {
        return {x * t.x - y * t.y, x * t.y + y * t.x};
    }

} a[4][N];

void fft(complex (&a)[N], int inv) {
    for (int i=0, j=0, k; i<tot; ++i) {
        if(j > i) {complex t = a[i]; a[i] = a[j]; a[j] = t;}
        for (k = tot; j & (k >>= 1); j &= ~k);
        j |= k;
    }
    for(int step=1; step<tot; step<<=1) {
        double alpha = inv*M_PI / step;
        for(int k=0; k<step; k++) {
            complex wk = {cos(alpha*k), sin(alpha*k)};
            for(int Ek=k; Ek<tot; Ek += step<<1) {
                int Ok = Ek + step; complex t = wk * a[Ok];
                a[Ok] = a[Ek] - t; a[Ek] += t;
            }
        }
    }
}

void solve() {
    int bit = 0; tot = (nb<<1) | 1;
    while ((1 << bit) < tot) ++bit;
    tot = 1 << bit;
    for (int i=0; i<tot; ++i) a[0][i] = a[1][i] = a[2][i] = a[3][i] = {i<nb && f[i] ? 1. : 0., 0.};
    while (nc--) {
        int v; char h; cin >> v >> h;
        a[h == 'S' ? 0 : (h == 'H' ? 1 : (h == 'C' ? 2 : 3))][v].x = 0.;
    }
    for (int i=0; i<4; ++i) fft(a[i], 1);
    for (int i=0; i<tot; ++i) a[0][i] = a[0][i] * a[1][i], a[2][i] = a[2][i] * a[3][i];
    fft(a[0], -1); fft(a[2], -1);
    for (int i=0; i<tot; ++i) a[0][i] = {i<nb ? a[0][i].x / tot : 0., 0.}, a[2][i] = {i<nb ? a[2][i].x / tot : 0., 0.};
    fft(a[0], 1); fft(a[2], 1);
    for (int i=0; i<tot; ++i) a[0][i] = a[0][i] * a[2][i];
    fft(a[0], -1);
    for (int i=na; i<=nb; ++i) cout << (long long)(a[0][i].x / tot + .5) << endl;
    cout << endl;
}

int main() {
    for (int i=2; i*i<M; ++i) if (!f[i]) for (int j=i*i; j<M; j+=i) f[j] = true;
    while (cin >> na >> nb >> nc && (na || nb || nc)) solve();
    return 0;
}