CTime_Pup_314

不能再颓废了

0%

[HNOI2019]白兔之舞

发表于 更新于 分类于 ， 阅读次数： Valine：

很显然，设答案序列为 \(a\) 则有

\[ a_t=\sum_{i=0}^L\binom{L}{i}w^i_{x,y}[i\operatorname{mod}k=t] \]

由单位根反演

\[ [n|k]=\frac{1}{n}\sum_{i=0}^{n-1}\omega_n^{ki} \] 带入答案式子就有

\[ a_t=\sum_{i=0}^L\binom{L}{i}w_{x,y}^i\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{(i-t)j} \]

我们交换求和顺序

\[ a_t=\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-tj}\sum_{i=0}^L\binom{L}{i}w_{x,y}^i\omega_k^{ij} \]

后面的部分用二项式定理得

\[ a_t=\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-tj}(\omega_k^jw+e)^L_{x,y} \]

我们考虑多项式

\[ f(x)=\frac{1}{k}\sum_{j=0}^{k-1}x^j(\omega_k^jw+e)^L_{x,y} \]

那么有

\[ a_t=f(\omega_k^{-t}) \]

单位根的求值可以用 FFT 快速求出，但前提是在项数为 \(2\) 的整数次幂的情况下，如果我们要做任意长度的 DFT 那么就需要用 Bluestein's Algorithm 具体细节如下

考虑 DFT 的形式

\[ y_j=\sum_{i=0}^{n-1}a_i(\omega_n^j)^i \]

根据 \(ij=\binom{i+j}{2}-\binom{i}{2}-\binom{j}{2}\) 可以得到

\[ y_j=\omega_n^{-\binom{j}{2}}\sum_{i=0}^{n-1}a_i\omega_n^{-\binom{i}{2}}\omega_n^{\binom{i+j}{2}} \]

如果设

\[ f=\sum_{i=0}^{n-1}a_i\omega_n^{-\binom{i}{2}}x^i \]

\[ g=\sum_{i=0}^{2n-2}\omega_n^{\binom{2n-2-i}{2}}x^i \]

则

\[ y_j=\omega_n^{-\binom{j}{2}}[x^{2n-2-j}]{(f\times g)} \]

而后者就是一个普通的卷积

接着，我们考虑 IDFT

\[ c_j=\frac{1}{n}\sum_{i=0}^{n-1}a_i\omega_n^{-ij} \]

其实与上面同理，只不过是符号变了罢了

若我们考虑模 \(p\) 意义下的原根 \(g\)，用 \(g^{\frac{p-1}{k}}\) 来代替 \(\omega_k\) 的话，就用拆系数 FTT 即可

计算一下总复杂度为 \(\operatorname{O}(\Delta+kn^3\log L+k\log k)\) 最前面是求原根的复杂度

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
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
#include <cstdio>
#include <cstring>
#include <cctype>
#include <vector>
#include <cmath>
#include <algorithm>
using namespace std;
typedef long long i64;
typedef double f64;
inline int read(int f = 1, int x = 0, char ch = ' ')
{
    while(!isdigit(ch = getchar())) if(ch == '-') f = -1;
    while(isdigit(ch)) x = x*10+ch-'0', ch = getchar();
    return f*x;
}
const int N = 3, K = 1<<18|5; const f64 PI = acos(-1);
int n, k, kinv, L, x, y, P, gn, g[K], f[K];
struct Mat
{
    i64 A[N][N];
    Mat() { memset(A, 0, sizeof(A)); }
    i64* operator [] (const int i) { return A[i]; }
    friend Mat operator * (Mat A, Mat B)
    {
        Mat C;
        for(int i = 0; i < n; ++i)
            for(int j = 0; j < n; ++j)
                for(int k = 0; k < n; ++k)
                    C[i][j] = (C[i][j]+A[i][k]*B[k][j]%P)%P;
        return C;
    }
    friend Mat operator * (int k, Mat A)
    {
        for(int i = 0; i < n; ++i)
            for(int j = 0; j < n; ++j)
                A[i][j] = k*A[i][j]%P;
        return A;       
    }
    friend Mat operator + (Mat A, Mat B)
    {
        for(int i = 0; i < n; ++i)
            for(int j = 0; j < n; ++j)
                A[i][j] = (A[i][j]+B[i][j])%P;
        return A;
    }
    void set() { for(int i = 0; i < n; ++i) A[i][i] = 1; }
    friend Mat operator ^ (Mat A, int b)
    {
        Mat ret; ret.set();
        for( ; b; b >>= 1, A = A*A) if(b&1) ret = ret*A;
        return ret;
    }
}w, e;
i64 qpow(i64 a, int b) { i64 ret = 1; for( ; b; b >>= 1, a = a*a%P) if(b&1) ret = ret*a%P; return ret; }
void groot()
{
    vector<int> p; int n = P-1;
    for(int i = 2, t = sqrt(n); i <= n; ++i)
        if(n%i == 0) { p.push_back(i); while(n%i == 0) n /= i; }
    if(n != 1) p.push_back(n);
    for(gn = 1; ; ++gn)
    {
        int i = 0; for( ; i < p.size()&&qpow(gn, (P-1)/p[i]) != 1; ++i);
        if(i == p.size()) break;
    }
    gn = qpow(gn, (P-1)/k), g[0] = 1;
    for(int i = 1; i < k; ++i) g[i] = 1ll*g[i-1]*gn%P;
}
namespace Bluestein
{
struct C
{
    f64 a, b;
    C(f64 a = 0, f64 b = 0):a(a), b(b) {};
    friend C operator + (C a, C b) { return C(a.a+b.a, a.b+b.b); }
    friend C operator - (C a, C b) { return C(a.a-b.a, a.b-b.b); }
    friend C operator * (C a, C b) { return C(a.a*b.a-a.b*b.b, a.a*b.b+a.b*b.a); }
    C operator ~ () { return C(a, -b); }
}w[K], S[K], T[K];
int lim, rev[K], C2[K];
void prepare(int ti)
{
    for(lim = 1; lim <= ti; lim <<= 1);
    for(int i = 0, j = 0; i < lim; ++i)
    {
        w[i] = C(cos(2*PI*i/lim), sin(2*PI*i/lim)), rev[i] = j;
        for(int k = lim>>1; (j ^= k) < k; k >>= 1);
    }
}
struct Poly
{
    vector<int> A; vector<C> B;
    int& operator [] (const int i) { return A[i]; }
    void set(int ti) { A.resize(ti+1); }
    int ti() { return A.size()-1; }
    void clear() { return B.clear(); } 
    void init() 
    { 
        int n = ti(); B.resize(n+1); 
        for(int i = 0; i <= n; ++i) B[i] = C(A[i]>>15, A[i]&32767);
    }
    void FFT(int t)
    {
        if(!t)
        {
            B.resize(lim); for(int i = 0; i < lim; ++i) if(i < rev[i]) swap(B[i], B[rev[i]]);
            for(int mid = 1; mid < lim; mid <<= 1)
                for(int j = 0, len = mid<<1; j < lim; j += len)
                    for(int k = 0, p = 0, q = lim/len; k < mid; ++k, p += q)
                    {
                        C x = B[j+k], y = w[p]*B[j+k+mid];
                        B[j+k] = x+y, B[j+k+mid] = x-y;
                    }
        }
        else
        {
            reverse(++B.begin(), B.end()), FFT(0); C v(1.0/lim, 0);
            for(int i = 0; i < lim; ++i) B[i] = B[i]*v;
        }
    } 
    friend Poly operator * (Poly A, Poly B)
    {
        int n = A.ti(), m = B.ti(); prepare(n+m), A.init(), B.init();
        A.FFT(0), B.FFT(0); C p, q, a, b, c, d;
        for(int i = 0; i < lim; ++i)
        {
            p = A.B[i], q = ~A.B[i?lim-i:0], a = (p+q)*C(0.5, 0), b = (p-q)*C(0, -0.5);
            p = B.B[i], q = ~B.B[i?lim-i:0], c = (p+q)*C(0.5, 0), d = (p-q)*C(0, -0.5);
            S[i] = a*c+b*d*C(0, 1), T[i] = a*d+b*c;
        }
        for(int i = 0; i < lim; ++i) A.B[i] = S[i], B.B[i] = T[i];
        A.FFT(1), B.FFT(1), A.set(n+m);
        for(int i = 0; i <= n+m; ++i)
        {
            i64 a = A.B[i].a+0.5, b = B.B[i].a+0.5, c = A.B[i].b+0.5;
            A[i] = (a%P*(1<<30)+b%P*(1<<15)+c)%P, A[i] = (A[i]+P)%P;
        }
        A.clear(), B.clear(); return A;
    }
};
void solve(int *a, int n)
{
    Poly F, G; F.set(n-1), G.set((n-1)<<1);
    for(int i = 1; i < (n-1)<<1; ++i) C2[i+1] = (C2[i]+i)%n;
    for(int i = 0; i < n; ++i) F[i] = 1ll*a[i]*g[C2[i]?n-C2[i]:0]%P;
    for(int i = 0; i <= (n-1)<<1; ++i) G[i] = g[C2[(n-1)*2-i]];
    F = F*G; for(int i = 0; i < n; ++i) a[i] = 1ll*g[C2[i]?n-C2[i]:0]*F[(n-1)*2-i]%P;
}
}
int main()
{
    n = read(), k = read(), L = read(), x = read()-1, y = read()-1, P = read();
    kinv = qpow(k, P-2), e.set(), groot();
    for(int i = 0; i < n; ++i) for(int j = 0; j < n; ++j) w[i][j] = read();
    for(int i = 0; i < k; ++i) f[i] = kinv*((g[i]*w+e)^L)[x][y]%P;   
    Bluestein::solve(f, k);
    for(int i = 0; i < k; ++i) printf("%d\n", f[i?k-i:0]);
    return 0;
}