多校省选模拟33 C

令 \(T = 2*3*5*7*11*47\)

\(f[i][j]\) 表示确定 \(i\) 位数，模 \(T\) 的余数为 \(j\) 的方案数（类似于数位 \(dp\) 直接转移就好了）。

然后发现 \(T\) 可以优化， \(T = 3*7*11*47\) ，我们只要让最终的数不被这 \(4\) 个数整除，然后最后一位强制不填 \(2,5\)，只填 \(1,3,7\)

然后倍增优化

令\(pw = 10^\frac{i}{2}\) , 则有 \(f[i][j*pw+k] = \sum f[\frac{i}{2}][j] * f[\frac{i}{2}][k]\)

本来需要 \(n*T\) 次转移，现在需要 \(T^2logn\) 次，但还能优化那个 \(T^2\)

固定后面的 \(f[\frac{i}{2}][k]\) （可以将其看为 \(h[k]\) ），对于所有合法的的 \(f[\frac{i}{2}][j]\) 求和得到 \(g\) 数组，就可以将 \(T^2\) 的转移用 \(FFT\) 优化到 \(TlogT\) 了

\(g[j] = \sum\limits_{x*pw\equiv j (mod T)} f[\frac{i}{2}][x]\)

\(f[i][j+k] = \sum g[j]*h[k]\)

时间复杂度 \(O(TlogTlogn)\)

注意事项：

• 快速幂算 \(10^{\frac{i}{2}}\) 时千万记得模的是 \(T\)

• \(f[i][j]\) 应该是卷积结果的第 \(j\) 项以及第 \(j+T\) 项的和

• 因为 \(\frac{i}{2}\) 是下取整，如果 \(i\) 是奇数，用 \(f[\frac{i}{2}]\) 转移出来的是 \(f[i-1]\) ，所以再暴力 \(O(T)\) 转移一下得到 \(f[i]\)

代码

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define ll long long
#define db long double
#define rint register int
using namespace std;

const int mod = 9973, T = 10857, N = 1<<15|10;
const double pi = M_PI;
int n, rev[N], nxt[T+10][5], f[T+10], lim, bit;
ll ff[T+10];
struct Com {
    db x, y;
    Com operator + (const Com &B) { return (Com){x + B.x, y + B.y}; }
    Com operator - (const Com &B) { return (Com){x - B.x, y - B.y}; }
    Com operator * (const Com &B) { return (Com){x * B.x - y * B.y, x * B.y + y * B.x}; }
} g[N], h[N];

int Pow(rint a, rint x, rint ans = 1) {
    for(;x;x >>= 1, a = 1ll * a * a % T)
        if(x&1) ans = 1ll * ans * a % T;
    return ans;
}

void FFT(Com *a, rint opt) {
    for(rint i = 0;i < lim; ++i) if(i < rev[i]) swap(a[i], a[rev[i]]);
    for(rint d = 1;d < lim;d <<= 1) {
        rint st = d<<1;
        Com w1 = (Com){cos(pi/d), sin(pi/d) * opt};
        for(rint i = 0;i < lim; i += st) {
            Com w = (Com){1, 0};
            for(rint j = 0;j < d; ++j, w = w * w1) {
                Com x = a[i+j], y = w * a[i+j+d];
                a[i+j] = x + y;
                a[i+j+d] = x - y;
            }
        }
    }
    if(~opt) return;
    for(rint i = 0;i < lim; ++i) a[i].x = a[i].x/lim+0.5;
}

void Solve(rint x) {
    if(!x) return f[0] = 1, void();
    Solve(x/2);
    rint pw = Pow(10, x/2);
    for(rint i = 0;i < lim; ++i) h[i] = g[i] = (Com){0, 0};
    for(rint i = 0;i < T; ++i) h[i].x = f[i], g[1ll*pw*i%T].x += f[i];
    FFT(h, 1), FFT(g, 1);
    for(rint i = 0;i < lim; ++i) h[i] = h[i] * g[i];
    FFT(h, -1);
    for(rint i = 0;i < T; ++i) f[i] = ((ll)h[i].x + (ll)h[i+T].x) % mod;
    if(x&1) {
        memset(ff, 0, sizeof(ff));
        for(rint i = 0;i < T; ++i) {
            ff[nxt[i][0]] += f[i];
            ff[nxt[i][1]] += f[i];
            ff[nxt[i][2]] += f[i];
            ff[nxt[i][3]] += f[i];
            ff[nxt[i][4]] += f[i];
        }
        for(rint i = 0;i < T; ++i) f[i] = ff[i] % mod;
    }
}

int main() {
    scanf("%d", &n);
    lim = 1, bit = 0;
    while(lim <= T*2) lim <<= 1, ++bit; --bit;
    for(rint i = 0;i < lim; ++i) rev[i] = (rev[i>>1]>>1) | ((i&1)<<bit);
    for(rint i = 0;i < T; ++i) {
        nxt[i][0] = (10*i+1)%T;
        nxt[i][1] = (10*i+2)%T;
        nxt[i][2] = (10*i+3)%T;
        nxt[i][3] = (10*i+5)%T;
        nxt[i][4] = (10*i+7)%T;
    }
    Solve(n-1);
    memset(ff, 0, sizeof(ff));
    for(rint i = 0;i < T; ++i) {
        ff[nxt[i][0]] += f[i];
        ff[nxt[i][2]] += f[i];
        ff[nxt[i][4]] += f[i];
    }
    rint ret = 0;
    for(rint i = 0;i < T; ++i) if(i % 3 && i % 7 && i % 11 && i % 47) ret += ff[i] % mod;
    printf("%d\n", ret % mod);
    return 0;
}