链接:
题解:
考虑处理 ,作射线 ,从原点出发,碰到 记一个 ,碰到 记一个 ,然后要处理出每个 之前的 的个数。
当 的时候,将 互换变成 。
当 的时候,记 ,那么可以将 个 和一个 替换成 ,递归变成 。
大概就建出了一个线段树之类的结构,就可以对着这个解构了。
对于 ,首先可以将 变得互质,那么求出 ,那么可以认为是解构区间 ,在那个类似线段树的结构上递归做就行了。
代码:
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 30;
const int M = 10000;
const int mod = 998244353;
ll p, q, r, l;
int n, total;
ll add(ll x, ll y, ll mod) {
x += y;
if (x >= mod) {
x -= mod;
}
return x;
}
ll mul(ll x, ll y, ll mod) {
ll result = 0;
for (; y; y >>= 1, x = add(x, x, mod)) {
if (y & 1) {
result = add(result, x, mod);
}
}
return result;
}
void exgcd(ll a, ll b, ll &x, ll &y) {
if (!b) {
x = 1;
y = 0;
} else {
exgcd(b, a % b, y, x);
y -= a / b * x;
}
}
ll div(ll x, ll y, ll mod) {
ll a, b;
exgcd(y, mod, a, b);
a = (a % mod + mod) % mod;
return mul(x, a, mod);
}
int add(int x, int y) {
x += y;
if (x >= mod) {
x -= mod;
}
return x;
}
int sub(int x, int y) {
x -= y;
if (x < 0) {
x += mod;
}
return x;
}
int mul(int x, int y) {
return (ll)x * y % mod;
}
int power(int x, int y) {
int result = 1;
for (; y; y >>= 1, x = mul(x, x)) {
if (y & 1) {
result = mul(result, x);
}
}
return result;
}
struct matrix_t {
int a[N][N];
matrix_t() {
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= n; ++j) {
a[i][j] = 0;
}
}
}
matrix_t operator + (const matrix_t &b) const {
matrix_t c;
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= n; ++j) {
c.a[i][j] = add(a[i][j], b.a[i][j]);
}
}
return c;
}
matrix_t operator * (const matrix_t &b) const {
matrix_t c;
for (int k = 1; k <= n; ++k) {
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= n; ++j) {
c.a[i][j] = add(c.a[i][j], mul(a[i][k], b.a[k][j]));
}
}
}
return c;
}
void print() {
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= n; ++j) {
printf("%d%c", a[i][j], j == n ? '\n' : ' ');
}
}
}
} a, b;
matrix_t indentity() {
matrix_t result;
for (int i = 1; i <= n; ++i) {
result.a[i][i] = 1;
}
return result;
}
matrix_t power(matrix_t x, ll y) {
matrix_t result = indentity();
for (; y; y >>= 1, x = x * x) {
if (y & 1) {
result = result * x;
}
}
return result;
}
struct node_t {
matrix_t power_a, power_b, answer;
int l, r;
ll zero;
node_t() {
l = r = zero = 0;
answer = matrix_t();
power_a = power_b = indentity();
}
node_t(int type) {
if (type == 1) {
zero = 0;
answer = matrix_t();
power_a = indentity();
power_b = b;
} else {
zero = 1;
answer = power_a = a;
power_b = indentity();
}
l = r = 0;
}
node_t operator + (const node_t &b) const {
node_t result;
result.zero = zero + b.zero;
result.power_a = power_a * b.power_a;
result.power_b = power_b * b.power_b;
result.answer = answer + power_a * b.answer * power_b;
return result;
}
} pool[M];
void push_up(int x) {
int l = pool[x].l, r = pool[x].r;
pool[x] = pool[l] + pool[r];
pool[x].l = l;
pool[x].r = r;
}
int build_power(ll value, int id) {
if (value == 1) {
return id;
}
int node = ++total;
if (value & 1) {
pool[node].l = build_power(value - 1, id);
pool[node].r = id;
} else {
pool[node].l = pool[node].r = build_power(value >> 1, id);
}
push_up(node);
return node;
}
int build(ll p, ll q) {
int zero, one;
bool flag = false;
pool[zero = ++total] = node_t(0);
pool[one = ++total] = node_t(1);
while (p != q) {
if (p < q) {
swap(p, q);
swap(zero, one);
flag = !flag;
} else {
ll div = (p - 1) / q;
int node = ++total;
pool[node].l = build_power(div, one);
pool[node].r = zero;
push_up(node);
zero = node;
p -= q * div;
}
}
if (flag) {
swap(zero, one);
}
++total;
pool[total].l = one;
pool[total].r = zero;
push_up(total);
return total;
}
node_t power(node_t x, ll y) {
node_t result;
for (; y; y >>= 1, x = x + x) {
if (y & 1) {
result = result + x;
}
}
return result;
}
node_t solve_l(int x) {
if (!x) {
return node_t();
}
if (!pool[x].zero) {
return pool[x];
}
if (pool[pool[x].l].zero) {
return solve_l(pool[x].l);
} else {
return pool[pool[x].l] + solve_l(pool[x].r);
}
}
node_t solve_r(int x) {
if (!x) {
return node_t();
}
if (!pool[x].zero) {
return pool[x];
}
if (pool[pool[x].r].zero) {
return solve_r(pool[x].r);
} else {
return solve_r(pool[x].l) + pool[pool[x].r];
}
}
node_t solve(int x, ll l, ll r) {
if (l == 1 && r == pool[x].zero) {
return pool[x];
}
node_t answer;
if (r > pool[x].zero) {
if ((l - 1) / pool[x].zero == (r - 1) / pool[x].zero) {
if (l % pool[x].zero == 0 && l) {
answer = answer + solve_r(x);
}
answer = answer + solve(x, (l - 1) % pool[x].zero + 1, (r - 1) % pool[x].zero + 1);
return answer;
} else {
if (l % pool[x].zero == 0 && l) {
answer = answer + solve_r(x);
}
answer = answer + solve(x, (l - 1) % pool[x].zero + 1, pool[x].zero);
answer = answer + power(pool[x], (r - 1) / pool[x].zero - (l - 1) / pool[x].zero - 1);
answer = answer + solve(x, 1, (r - 1) % pool[x].zero + 1);
return answer;
}
}
if (l > pool[pool[x].l].zero) {
if (l == pool[pool[x].l].zero + 1) {
answer = answer + solve_r(pool[x].l);
}
answer = answer + solve(pool[x].r, l - pool[pool[x].l].zero, r - pool[pool[x].l].zero);
} else if (r <= pool[pool[x].l].zero) {
answer = answer + solve(pool[x].l, l, r);
if (r == pool[pool[x].l].zero) {
answer = answer + solve_l(pool[x].r);
}
} else {
answer = answer + solve(pool[x].l, l, pool[pool[x].l].zero);
answer = answer + solve(pool[x].r, 1, r - pool[pool[x].l].zero);
}
return answer;
}
int main() {
#ifdef wxh010910
freopen("input.txt", "r", stdin);
#endif
scanf("%lld %lld %lld %lld %d", &p, &q, &r, &l, &n);
pool[0] = node_t();
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= n; ++j) {
scanf("%d", &a.a[i][j]);
}
}
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= n; ++j) {
scanf("%d", &b.a[i][j]);
}
}
ll gcd = __gcd(p, q);
p /= gcd;
q /= gcd;
r /= gcd;
ll pos = div(r, p, q);
matrix_t answer = solve(build(p, q), pos + 1, pos + l).answer * power(b, r / q);
answer.print();
return 0;
}