TopCoder - 12118
单位根反演，指数生成函数。

题意

给你两个仅包含 a,b,c\texttt{a,b,c} 字符串 A,BA,B（∣A∣,∣B∣≤11|A|,|B|\le 11），以及一个长度为 33 的数组分别表示 a→b,b→c,c→a\texttt{a}\to\texttt{b},\texttt{b}\to\texttt{c},\texttt{c}\to\texttt{a} 的代价，和 Maxcost(≤109)Maxcost (\le 10^9)，问在不超过 MaxcostMaxcost 将 AA 变成 BB 的方案数。两种方案视为不同的，当存在一个位置，两个方案改变了不同的位置的字符。

题解

每个位置若 Ai≠BiA_i\ne B_i ，先要花费一定步数使其变到 BB 。然后每个位置再可以变三的倍数次。

题意转化为每个位置操作次数为 ci=3k+ai(ai∈{0,1,2},k∈N)c_i=3k+a_i(a_i\in\{0,1,2\},k\in \mathbb{N}) 。然后总次数 ∑ici≤N\sum_i c_i \le N 。且对于一个确定的操作次数序列 ⟨c1,c2,…,cn⟩\langle c_1,c_2,\dots,c_n\rangle 其贡献为 (∑ici)!∏i(ci!)\dfrac{\big(\sum_i c_i\big)!}{\prod_i(c_i!)} 。

设 Fj^(x)\hat{F_j}(x) 表示位置 jj 的指数生成函数，[xi]Fj^(x)[x^i]\hat{F_j}(x) 即第 jj 个位置操作次数为 ii 的方案数。

Fj^(x)=∑i≥0x3i+aj(3i+aj)!=∑i≥0[i mod 3=aj]xii!=∑i≥0[3∣(i−aj)]xii!=∑i≥013∑k=02ω3k(i−aj)xii!=13∑i≥0∑k=02ω3kiω3kajxii!=13∑k=021ω3kaj∑i≥0ω3kixii!\begin{aligned} \hat{F_j}(x)&=\sum_{i\ge 0} \frac{x^{3i+a_j}}{\big(3i+a_j\big)!}\\ &=\sum_{i\ge 0} [i\bmod 3=a_j] \frac{x^i}{i!}\\ &=\sum_{i\ge 0} [3\mid (i-a_j)] \frac{x^i}{i!}\\ &=\sum_{i\ge 0} \frac{1}{3}\sum_{k=0}^{2} \omega^{k(i-a_j)}_{3} \frac{x^i}{i!}\\ &=\frac{1}{3}\sum_{i\ge 0} \sum_{k=0}^{2} \frac{\omega_3^{ki}}{\omega_3^{ka_j}} \frac{x^i}{i!}\\ &=\frac{1}{3} \sum_{k=0}^{2}\frac{1}{\omega_3^{ka_j}}\sum_{i\ge 0} \omega_3^{ki} \frac{x^i}{i!}\\ \end{aligned}

令 G^(x)=∏j=1nFj^(x)\hat{G}(x)=\prod_{j=1}^n \hat{F_j}(x) ，gi=[xi]G^(x)g_i=[x^i]\hat{G}(x) 答案就等于 ∑i=0Ngi\sum_{i=0}^{N} g_i

G^(x)=∏j=1n(13∑k=021ω3kaj∑i≥0ω3kixii!)=13n∏j=1n(∑k=021ω3kaj∑i≥0ω3kixii!)=13n∏j=1n(∑k=02ω3−kajeω3kx)\begin{aligned} \hat{G}(x)&=\prod_{j=1}^n\left(\frac{1}{3} \sum_{k=0}^{2}\frac{1}{\omega_3^{ka_j}}\sum_{i\ge 0} \omega_3^{ki} \frac{x^i}{i!}\right)\\ &=\frac{1}{3^n}\prod_{j=1}^n\left( \sum_{k=0}^{2} \frac{1}{\omega_3^{ka_j}} \sum_{i\ge 0} \omega_{3}^{ki} \frac{x^i}{i!}\right)\\ &=\frac{1}{3^n}\prod_{j=1}^n\left(\sum_{k=0}^{2} \omega_3^{-ka_j} e^{\omega_{3}^{k}x} \right)\\ \end{aligned}

前面那项是个定值，观察后面这项就是 ∏j=1n(eω30x+ω3−ajeω31x+ω3−2ajeω32x)\prod _{j=1}^n\left(e^{\omega_3^0 x}+\omega_3^{-a_j}e^{\omega_3^1 x}+\omega_3^{-2a_j}e^{\omega_3^2 x} \right) 可以暴力展开 O(3n)O(3^n)，会得到若干 □e□x\Box e^{\Box x} 的和。然后每项可以这样展开：

p⋅eqx=∑i=0∞p qixii!p \cdot e^{qx}=\sum_{i=0}^{\infty} \frac{p\,q^ix^i}{i!}

然后对每个 p⋅eqxp \cdot e^{qx} 对应的指数生成函数求系数前缀和：

∑i=0Npqi=p×qN+1−1q−1\sum_{i=0}^N pq^i = p\times \frac{q^{N+1}-1}{q-1}

CODE

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
#include <bits/stdc++.h>
#define MOD(x) ((x)<mod?(x):((x)%mod))
template<typename _Tp>
inline void red(_Tp &x) {
    x=0;bool fg=0;char ch=getchar();
    while(ch<'0'||ch>'9') {if(ch=='-') fg^=1;ch=getchar();}
    while(ch>='0'&&ch<='9') x=(x<<1)+(x<<3)+(_Tp)(ch^48),ch=getchar();
    if(fg) x=-x;
}
template <typename _Tp,typename... Args>void red(_Tp &t,Args &... args){red(t);red(args...);}
using namespace std;
typedef long long ll;
const ll mod = 1e9+7;
const ll sq3 = 82062379; // 3 的二次剩余
const ll inv3 = 333333336;
struct com{
    ll x,y;
    com(ll _x=0,ll _y=0):x(_x),y(_y){}
    com operator+(const com &t) {return com((x+t.x)%mod,(y+t.y)%mod);}
    com operator+=(const com &t) {*this=*this+t;return *this;}
    com operator*(const com &t) {return com(((x*t.x-y*t.y)%mod+mod)%mod,(x*t.y+y*t.x)%mod);}
    com operator=(const ll &t) {*this=com(t%mod,0);return *this;}
};
ostream &operator<<(ostream &out,com a) {return out << "(" << a.x << ", " << a.y <<") ";}
class ConversionMachine {
private:
    ll fpow(ll a,ll b) {
        a%=mod;
        ll r=1; for(;b;b>>=1,a=MOD(a*a)) if(b&1) r=MOD(r*a);
        return r;
    }
    com inv(com t) {
        ll q=fpow((t.x*t.x+t.y*t.y)%mod,mod-2);
        return com(MOD(t.x*q),(mod-(t.y*q)%mod)%mod);
    }
    com fpow(com a,ll b) {
        com r=1;
        if(b<0) {a=inv(a); b=-b;}
        for(;b;b>>=1,a=a*a) if(b&1) r=r*a; 
        return r;
    }
    ll inv(ll x) {return fpow(x%mod,mod-2);}
    const com w=com((mod-1)/2,(mod+sq3)/2),iw=com((mod-1)/2,(mod-sq3)/2);
    com calc(com q,ll k) {
        if(q.x==1&&q.y==0) return com(k+1,0);
        return (fpow(q,k+1)+com(mod-1,0))*inv(q+com(mod-1,0));
    }
public:
    com f[2][32][32];
    int countAll(string A,string B,vector<int> cost,ll maxcost) {
        int n=A.size(),N=0;
        vector<int> a(n);
        com ans=com(1,0);
        for(int i=0;i<n;++i) {
            while(A[i]!=B[i]) {
                maxcost-=cost[A[i]-'a'];
                A[i]='a'+(A[i]-'a'+1)%3;
                a[i]++;
            }
            N+=a[i];
            ans=ans*com(inv3,0);
        }
        if(maxcost<0) return 0;
        N+=maxcost/(cost[0]+cost[1]+cost[2])*3;
        memset(f,0,sizeof(f));
        f[0][0][0]=1;
        bool V=0;
        for(int i=0;i<n;++i) {
            memset(f[V^1],0,sizeof(f[V^1]));
            for(int j=0;j<=i;++j) {
                for(int k=0;j+k<=i;++k) {
                    f[V^1][j+1][k]+=f[V][j][k];
                    f[V^1][j][k+1]+=fpow(w,-a[i])*f[V][j][k];
                    f[V^1][j][k]+=fpow(w,-2*a[i])*f[V][j][k];
                }
            }
            V^=1;
        }
        com sum;
        for(int i=0;i<=n;++i) {
            for(int j=0;i+j<=n;++j) {
                com q=com(i,0)+com(j,0)*w+com(n-i-j)*iw; 
                sum=sum+f[V][i][j]*calc(q,N);
            }
        }
        ans=ans*sum;
        return ans.x;
    }
public:
    void run_test(int Case) { if ((Case == -1) || (Case == 0)) test_case_0(); if ((Case == -1) || (Case == 1)) test_case_1(); if ((Case == -1) || (Case == 2)) test_case_2(); if ((Case == -1) || (Case == 3)) test_case_3(); }
private:
    template <typename T> string print_array(const vector<T> &V) { ostringstream os; os << "{ "; for (typename vector<T>::const_iterator iter = V.begin(); iter != V.end(); ++iter) os << '\"' << *iter << "\","; os << " }"; return os.str(); }
    void verify_case(int Case, const int &Expected, const int &Received) { cerr << "Test Case #" << Case << "..."; if (Expected == Received) cerr << "PASSED" << endl; else { cerr << "FAILED" << endl; cerr << "\tExpected: \"" << Expected << '\"' << endl; cerr << "\tReceived: \"" << Received << '\"' << endl; } }
    void test_case_0() { string Arg0 = "a"; string Arg1 = "b"; int Arr2[] = {100,2,3}; vector <int> Arg2(Arr2, Arr2 + (sizeof(Arr2) / sizeof(Arr2[0]))); int Arg3 = 205; int Arg4 = 2; verify_case(0, Arg4, countAll(Arg0, Arg1, Arg2, Arg3)); }
    void test_case_1() { string Arg0 = "abcba"; string Arg1 = "abcba"; int Arr2[] = {67,23,43}; vector <int> Arg2(Arr2, Arr2 + (sizeof(Arr2) / sizeof(Arr2[0]))); int Arg3 = 0; int Arg4 = 1; verify_case(1, Arg4, countAll(Arg0, Arg1, Arg2, Arg3)); }
    void test_case_2() { string Arg0 = "cccccccc"; string Arg1 = "aaaaaaaa"; int Arr2[] = {10000000,1,1}; vector <int> Arg2(Arr2, Arr2 + (sizeof(Arr2) / sizeof(Arr2[0]))); int Arg3 = 100; int Arg4 = 40320; verify_case(2, Arg4, countAll(Arg0, Arg1, Arg2, Arg3)); }
    void test_case_3() { string Arg0 = "aa"; string Arg1 = "cc"; int Arr2[] = {1,10,100}; vector <int> Arg2(Arr2, Arr2 + (sizeof(Arr2) / sizeof(Arr2[0]))); int Arg3 = 1787; int Arg4 = 123611681; verify_case(3, Arg4, countAll(Arg0, Arg1, Arg2, Arg3)); }
}G;

int main() {
    G.run_test(-1);
      return 0;
}