[单位根反演] HDU 7013 String Mod

题目大意

对于所有由前 (k) 个小写英文字母组成的长为 (L) 的字符串，设 (A[i][j]) 表示字符串中 a 的数量模 (n) 等于 (i)，b 的数量模 (n) 等于 (j) 这样的字符串个数。输出矩阵 (A)。(2leq kleq 26,1leq Lleq 10^{18},1leq nleq 500)。

题解

设字符串中含有 (x) 个 a 和 (y) 个 b，对 (x,y) 进行枚举，有

[A[i][j]=sum_{x=0}^Lsum_{y=0}^{L-x} [n|x-i][n|y-j]inom{L}{x}inom{L-x}{y}(k-2)^{L-x-y}\ ]

考虑换掉艾弗森约定，可以使用单位根反演：

单位根反演

[forall k,[n|k]=frac{1}{n}sum_{i=0}^{n-1} omega_n^{ik} ]

证明：

1. 当 (n|k) 时，(omega_n^{ik}=omega_n^0=1)，(frac{1}{n}sum_{i=0}^{n-1}omega_n^{ik}=1=[n|k])

2. 当 (n mid k) 时，

[frac{1}{n}sum_{i=0}^{n-1}omega_n^{ik}=frac{1}{n}cdotfrac{omega_n^{nk}-1}{omega_n^k-1}=frac{1}{n}cdotfrac{1-1}{omega_n^k-1}=0=[n|k] ]

所以有

[A[i][j]=sum_{x=0}^Lsum_{y=0}^{L-x} [n|x-i][n|y-j]inom{L}{x}inom{L-x}{y}(k-2)^{L-x-y}\ =sum_{x=0}^Lsum_{y=0}^{L-x} left(frac{1}{n}sum_{p=0}^{n-1}omega_n^{p(x-i)} ight)left(frac{1}{n}sum_{q=0}^{n-1}omega_n^{q(y-j)} ight)inom{L}{x}inom{L-x}{y}(k-2)^{L-x-y}\ =frac{1}{n^2}sum_{p=0}^{n-1}sum_{q=0}^{n-1}sum_{x=0}^Lsum_{y=0}^{L-x} omega_n^{p(x-i)}omega_n^{q(y-j)}inom{L}{x}inom{L-x}{y}(k-2)^{L-x-y}\ =frac{1}{n^2}sum_{p=0}^{n-1}sum_{q=0}^{n-1}sum_{x=0}^Lsum_{y=0}^{L-x} omega_n^{px}omega_n^{qy}omega_n^{-pi}omega_n^{-qj}inom{L}{x}inom{L-x}{y}(k-2)^L(k-2)^{-x}(k-2)^{-y}\ =frac{1}{n^2}sum_{p=0}^{n-1}sum_{q=0}^{n-1}omega_n^{-pi}omega_n^{-qj}sum_{x=0}^Lsum_{y=0}^{L-x} omega_n^{px}omega_n^{qy}inom{L}{x}inom{L-x}{y}(k-2)^{L-x-y}\ =frac{1}{n^2}sum_{p=0}^{n-1}sum_{q=0}^{n-1}omega_n^{-pi}omega_n^{-qj}sum_{x=0}^Linom{L}{x}omega_n^{px}sum_{y=0}^{L-x} inom{L-x}{y}(k-2)^{L-x-y}(omega_n^q)^y\ =frac{1}{n^2}sum_{p=0}^{n-1}sum_{q=0}^{n-1}omega_n^{-pi}omega_n^{-qj}sum_{x=0}^Linom{L}{x}(omega_n^q+k-2)^{L-x}(omega_n^p)^x\ =frac{1}{n^2}sum_{p=0}^{n-1}sum_{q=0}^{n-1}omega_n^{-pi}omega_n^{-qj}(omega_n^p+omega_n^q+k-2)^L ]

令 (P(i,p)=omega_{n}^{-pi},Q(p,q)=(omega_n^p+omega_n^q+k-2)^L,R(q,j)=omega_n^{-qj})

则 (A[i][j]=frac{1}{n^2}sum_{p=0}^{n-1}sum_{q=0}^{n-1}P(i,p)Q(p,q)R(q,j))

可以作两次矩阵乘法得到答案。

Code

#include <bits/stdc++.h>
using namespace std;

#define RG register int
#define LL long long

const LL p = 1e9 + 9;

LL qpow(LL b, LL n) {
    LL x = 1, Power = b % p;
    while (n) {
        if (n & 1) x = x * Power % p;
        Power = Power * Power % p;
        n >>= 1;
    }
    return x;
}

struct Matrix {
    LL a[502][502];
    int n;
    Matrix operator*(const Matrix& r) {
        Matrix res;
        res.n = n;
        memset(res.a, 0, sizeof(res.a));
        for (int k = 0;k < n;++k)
            for (int i = 0;i < n;++i)
                for (int j = 0;j < n;++j)
                    res.a[i][j] = (res.a[i][j] + a[i][k] * r.a[k][j] % p) % p;
        return res;
    }
};
LL x[505], L;
int q[505];
Matrix ans, A, B, C;
int T, k, n;

void solve() {
    ans.n = A.n = B.n = C.n = n;
    LL g = qpow(13, (p - 1) / n);
    x[0] = 1;
    for (int i = 1;i < n;++i) x[i] = x[i - 1] * g % p;
    for (int i = 0;i < n;++i)
        for (int j = 0;j < n;++j) {
            A.a[i][j] = qpow(g, i * j);
            B.a[i][j] = qpow(x[i] + x[j] + k - 2, L);
        }
    C = A;
    for (int i = 0;i < n;++i)
        for (int j = 0;j < i;++j)
            swap(C.a[i][j], C.a[j][i]);
    ans = A * B;ans = ans * C;
    LL m = qpow(n * n, p - 2);
    for (int i = 1;i < n;++i) q[i] = n - i;
    for (int i = 0;i < n;++i) {
        for (int j = 0;j < n;++j) {
            printf("%lld", ans.a[q[i]][q[j]] * m % p);
            if (j < n - 1) printf(" ");
        }
        printf("
");
    }
    return;
}

int main() {
    scanf("%d", &T);
    while (T--) {
        scanf("%d%lld%d", &k, &L, &n);
        solve();
    }
    return 0;
}