LOJ #6440. 万能欧几里得

链接：

link

题解：

考虑处理 $\frac{Ax}{B}$ ，作射线 $y = \frac{A}{B}x$ ，从原点出发，碰到 $y = c$ 记一个 $1$ ，碰到 $x = c$ 记一个 $0$ ，然后要处理出每个 $0$ 之前的 $1$ 的个数。

当 $A < B$ 的时候，将 $01$ 互换变成 $10$ 。

当 $A\ge B$ 的时候，记 $k = \lfloor \frac{A}{B}\rfloor$ ，那么可以将 $k$ 个 $1$ 和一个 $0$ 替换成 $0$ ，递归变成 $(B, A\bmod B)$ 。

大概就建出了一个线段树之类的结构，就可以对着这个解构了。

对于 $\frac{Ax+C}{B}$ ，首先可以将 $(A, B)$ 变得互质，那么求出 $x = C\times A^{-1}\mod B$ ，那么可以认为是解构区间 $[x + 1, x + L]$ ，在那个类似线段树的结构上递归做就行了。

代码：

#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;
}