题解 P5293 【[HNOI2019]白兔之舞】
officeyutong
2019-04-25 00:29:59
蒟蒻并不会单位根反演，但是从朴素的DP也能推出来一样的东西qwq
(如果图片无法显示可以去https://yutong.site/?p=1388)
题外话:洛谷什么时候能支持多行Latex啊qwq..
```cpp
#include <assert.h>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <vector>
// #include
using int_t = long long int;
using std::cin;
using std::cout;
using std::endl;
using real_t = double;
using cpx_t = struct complex;
const int_t LARGE = 5e5;
const real_t PI = acos(-1);
struct complex {
real_t real, imag;
complex(real_t a = 0, real_t b = 0) : real(a), imag(b) {}
complex operator+(const complex& rhs) const {
return complex{real + rhs.real, imag + rhs.imag};
}
complex operator-(const complex& rhs) const {
return complex{real - rhs.real, imag - rhs.imag};
}
complex operator*(const complex& rhs) const {
return complex{real * rhs.real - imag * rhs.imag,
real * rhs.imag + imag * rhs.real};
}
complex& operator*=(const complex& rhs) {
*this = (*this) * rhs;
return *this;
}
complex conj() .{ return complex{real, -imag}; }
};
void transform(cpx_t* A, int size2, int arg);
void poly_mul(const int* A, int n, const int* B, int m, int* A0, int p);
// void transformX(int* A, int len, int g, int mod);
void transformNTT(int* A, int size2, int arg, int mod, int g);
std::vector<int> poly_dc_mul(const std::vector<int>& A);
void poly_inv(const int* A, int n, int* result);
void poly_div(const int* A, int n, const int* B, int m, int* R);
std::vector<int> poly_eval(const std::vector<int>& poly,
const std::vector<int>& vals);
void transformX(int* A0, int len, int g);
int power(int_t base, int_t index, int_t mod) {
int result = 1;
// base = (base % mod + mod) % mod;
index = (index % (mod - 1) + mod - 1) % (mod - 1);
while (index) {
if (index & 1) result = (int_t)result * base % mod;
index >>= 1;
base = (int_t)base * base % mod;
}
assert(result >= 0);
return result;
}
int_t mod;
struct Matrix {
int_t n;
int_t data[4][4];
int_t* operator[](int_t r) { return data[r]; }
int_t at(int_t r, int_t c) const { return data[r][c]; }
Matrix(int_t n) {
this->n = n;
memset(data, 0, sizeof(data));
}
Matrix operator*(const Matrix& rhs) const {
Matrix result(n);
for (int_t i = 1; i <= n; i++) {
for (int_t j = 1; j <= n; j++) {
for (int_t k = 1; k <= n; k++) {
result[i][j] =
(at(i, k) * rhs.at(k, j) % mod + mod + result[i][j]) %
mod;
}
}
}
return result;
}
Matrix pow(int_t index) const {
Matrix base = *this, result(n);
for (int_t i = 1; i <= n; i++) result[i][i] = 1;
while (index) {
if (index & 1) result = result * base;
base = base * base;
index >>= 1;
}
return result;
}
};
int_t n, k, l, x, y;
int_t mat[4][4];
Matrix makeMatrix(int_t x0) {
Matrix result(n);
for (int_t i = 1; i <= n; i++)
for (int_t j = 1; j <= n; j++) {
result[i][j] = x0 * mat[i][j] % mod;
if (i == j) result[i][j] += 1;
result[i][j] %= mod;
}
return result;
}
std::ostream& operator<<(std::ostream& os, const Matrix& mat) {
for (int_t i = 1; i <= mat.n; i++) {
for (int_t j = 1; j <= mat.n; j++) {
os << mat.at(i, j) << " ";
}
cout << endl;
}
return os;
}
int flips[20][LARGE];
const auto flip = [=](int x, int size2) {
int result = 0;
for (int i = 1; i < size2; i++) {
result |= (x & 1);
x >>= 1;
result <<= 1;
}
return result | (x & 1);
};
int main() {
for (int i = 1; i < 19; i++) {
for (int j = 1; j < (1 << i); j++) flips[i][j] = flip(j, i);
}
cin >> n >> k >> l >> x >> y >> mod;
for (int_t i = 1; i <= n; i++)
for (int_t j = 1; j <= n; j++) cin >> mat[i][j];
Matrix M0(n);
static int A[LARGE + 1];
M0[1][x] = 1;
int_t g = 2;
std::vector<int_t> divs;
int_t px = mod - 1;
assert(px % k == 0);
for (int_t i = 2; i * i <= px; i++) {
if (px % i == 0) {
divs.push_back(i);
if (i * i != px) divs.push_back(px / i);
}
}
for (; g <= mod - 2; g++) {
bool ok = true;
for (int_t x : divs) {
if (power(g, x, mod) == 1 && x != mod - 1) {
ok = false;
break;
}
}
if (ok) break;
}
for (int_t i = 0; i < k; i++) {
int_t gx = power(g, (mod - 1) / k * i, mod);
auto T = makeMatrix(gx);
A[i] = (M0 * T.pow(l))[1][y];
}
transformX(A, k, g);
const int_t leninv = power(k, mod - 2, mod);
for (int_t i = 0; i < k; i++) printf("%lld\n", A[i] * leninv % mod);
return 0;
}
void transformX(int* A0, int len, int g) {
const auto C2 = [](int_t n) -> int_t { return n * (n - 1) / 2; };
int root = power(g, (mod - 1) / len, mod);
static int A[LARGE + 1], B[LARGE + 1], A1[LARGE + 1];
for (int i = 0; i <= len * 2; i++) {
A[i] = power(root, -C2(i), mod);
}
for (int i = 0; i < len; i++) {
B[i] = (int_t)A0[i] * power(root, C2(i), mod) % mod;
}
std::reverse(B, B + len + 1);
poly_mul(A, len * 2, B, len, A1, mod);
for (int i = 0; i < len; i++)
A0[i] = (int_t)A1[i + len] % mod * power(root, C2(i), mod) % mod;
}
void poly_mul(const int* A, int n, const int* B, int m, int* A0, int p) {
int size2 = 0;
while ((1 << size2) < n + m + 1) size2++;
int len = (1 << size2);
const int px = 3e4;
static cpx_t C[LARGE + 1], D[LARGE + 1], G[LARGE + 1], Px[LARGE + 1];
for (int i = 0; i < len; i++) {
if (i <= n) {
C[i] = cpx_t{A[i] / px};
D[i] = cpx_t{A[i] % px};
} else {
C[i] = D[i] = cpx_t();
}
if (i <= m) {
G[i] = cpx_t{B[i] / px, B[i] % px};
} else {
G[i] = cpx_t{0, 0};
}
}
transform(C, size2, 1);
transform(D, size2, 1);
transform(G, size2, 1);
for (int_t i = 0; i < len; i++) {
C[i] = C[i] * G[i];
D[i] = D[i] * G[i];
}
transform(C, size2, -1);
transform(D, size2, -1);
const auto make = [=](real_t x) -> int_t {
assert(x / len >= -1);
return (x / len) + 0.5;
};
for (int i = 0; i <= n + m; i++) {
A0[i] = ((int_t)make(C[i].real) % p * px % p * px % p +
(int_t)make(C[i].imag + D[i].real) % p * px % p +
make(D[i].imag) % p) %
p;
}
}
void transform(cpx_t* A, int size2, int arg) {
int len = (1 << size2);
for (int i = 0; i < len; i++) {
int x = flips[size2][i];
if (x > i) std::swap(A[i], A[x]);
}
for (int i = 2; i <= len; i *= 2) {
for (int j = 0; j < len; j += i) {
for (int k = 0; k < i / 2; k++) {
auto u = A[j + k],
t = cpx_t(cos(2 * PI / i * k), sin(2 * PI / i * k) * arg) *
A[j + k + i / 2];
A[j + k] = u + t;
A[j + k + i / 2] = u - t;
}
}
}
}
```