万能欧几里得知识整理
题目
解析
以坐标原点为起点，将矩阵$A$看作向右走，$B$看作向上走，平面上的点$(x, y)$有贡献$A^xB^y$，答案就是$y=\left \lfloor \frac{Px+R}{Q} \right \rfloor(1 \leqslant x \leqslant L)$每一个$x$下第一个整点的贡献和。
来一发课件，下图中斜线段即为$y=\left \lfloor \frac{Px+R}{Q} \right \rfloor$
操作序列满⾜信息可合并。只需要两个操作序列的信息，即可计算出将这两个操作序列拼接后的信息。
操作可以被看成代码段：
初始化： matrix curA = $I$, curB = $B$ ^ floor(R/Q), sum = $0$
U： curB = curB * B;
R： curA = curA * A, sum += curA * curB;
因此一个操作序列段需要维护的信息有curA, curB, sum
实现一个函数$solve(P, Q, R, L, A, B)$， $A$、$B$为操作序列段，即由$U$和$R$组合而成的段。该函数代表了一个操作序列段，生成序列有$L$个$B$，编号为$1\sim L$。第$i-1$个$B$与第$i$个$B$之间有$\left \lfloor \frac{Pi+R}{Q} \right \rfloor-\left \lfloor \frac{P(i-1)+R}{Q} \right \rfloor$个$A$，也就是说这里是以原点为起点，将操作$B$看作向右走，操作$A$看作向上走，斜线为$y=\left \lfloor \frac{Pi+R}{Q} \right \rfloor$。
根据这个函数的定义，有如下两个性质：
R = R mod Q ：⽆影响
P = P mod Q ：B = A ^ floor(R/Q) * B
现在考虑$P<Q, R<Q$时的情况:
对于第$x$个$A$，若第$y$个$B$为第$x$个$A$后的$B$，则满足:$$x \leqslant \left \lfloor \frac{Py+R}{Q} \right \rfloor \\ Qx \leqslant Py+R \\ \left \lceil \frac{Qx-R}{P} \right \rceil \leqslant y \\ \left \lfloor \frac{Qx-R+P-1}{P} \right \rfloor \leqslant y $$
也就是说第$x$个$A$前有$\left \lfloor \frac{Qx-R-1}{P} \right \rfloor$个$B$，第$x-1$个$A$与第$x$个$A$之间有$\left \lfloor \frac{Qx-R-1}{P} \right \rfloor-\left \lfloor \frac{Q(x-1)-R-1}{P} \right \rfloor$个$B$，似乎可以递归了。但由于当$x$等于$1$时，$\left \lfloor \frac{Q(x-1)-R-1}{P} \right \rfloor < 0$，因此第一个$A$需要单独计算，假设有$m$个$A$，这样的话还需要计算的$x$的范围为$2\leqslant x \leqslant m$，但$solve$函数的$x$是从$1$开始枚举的，所以用$x-1$代替原来的$x$，于是就变成了，第$x-1$个$A$与第$x$个$A$之间有$\left \lfloor \frac{Qx+Q-R-1}{P} \right \rfloor-\left \lfloor \frac{Q(x-1)+Q-R-1}{P} \right \rfloor$个$B$，此时$x$的范围为$1\leqslant x \leqslant m - 1$，可以递归调用函数$solve(Q, P, Q - R - 1, m - 1, B, A)$。注意第$m$个$A$后可能还有一些$B$，需要在递归后统计贡献。
递归次数等于$gcd(P, Q)$
注意矩阵乘法的顺序。
代码：
#include<cstdio> #include<iostream> #include<algorithm> using namespace std; typedef long long ll; const int mod = 998244353; ll add(ll x, ll y) { return x + y < mod? x + y: x + y - mod; } int n; ll qpow(ll x, ll y) { ll ret = 1; while(y) { if(y&1) ret = ret * x % mod; x = x * x % mod; y >>= 1; } return ret; } struct Matrix{ ll a[25][25]; }dw, dw0; Matrix operator + (Matrix x, Matrix y) { for(int i = 1; i <= n; ++i) for(int j = 1; j <= n; ++j) x.a[i][j] = add(x.a[i][j], y.a[i][j]); return x; } Matrix operator * (Matrix x, Matrix y) { Matrix ret = dw0; for(int i = 1; i <= n; ++i) for(int j = 1; j <= n; ++j) for(int k = 1; k <= n; ++k) ret.a[i][j] = add(ret.a[i][j], x.a[i][k] * y.a[k][j] % mod); return ret; } Matrix mat_qpow(Matrix x, ll y) { Matrix ret = dw; while(y) { if(y&1) ret = ret * x; x = x * x; y >>= 1; } return ret; } struct node{ Matrix curA, curB, sum; node(){} node(Matrix x, Matrix y, Matrix z) { curA = x; curB = y; sum = z; } }; node operator + (node x, node y) { return node(x.curA * y.curA, x.curB * y.curB, x.sum + x.curA * y.sum * x.curB);//计算sum时注意顺序 } node qsum(node x, ll y) { node ret = node(dw, dw, dw0); while(y) { if(y&1) ret = ret + x; x = x + x; y >>= 1; } return ret; } ll calc(ll P, ll Q, ll R, ll L)//计算floor((P*L+R)/Q) { return (ll)((long double)P / Q * L + (long double)R / Q) ; } node solve(ll P, ll Q, ll R, ll L, node A, node B) { if(L == 0) return node(dw, dw, dw0); R %= Q; if(P >= Q) return solve(P % Q, Q, R, L, A, qsum(A, P / Q) + B); ll m = calc(P, Q, R, L); if(m <= 0) return qsum(B, L); ll cnt = L - calc(m, P, P - R - 1, Q) + 1; return qsum(B, (Q - R - 1) / P) + A + solve(Q, P, Q - R - 1, m - 1, B, A) + qsum(B, cnt); } int main() { ll P, Q, R, L; scanf("%lld%lld%lld%lld%d", &P, &Q, &R, &L, &n); Matrix A, B; for(int i = 1; i <= n; ++i) for(int j = 1; j <= n; ++j) scanf("%lld", &A.a[i][j]); for(int i = 1; i <= n; ++i) for(int j = 1; j <= n; ++j) scanf("%lld", &B.a[i][j]); for(int i = 1; i <= n; ++i) dw.a[i][i] = 1; node u = node(dw, B, dw0), r = node(A, dw, A); node ans = qsum(u, R / Q) + solve(P, Q, R, L, u, r); for(int i = 1; i <= n; ++i) { for(int j = 1; j <= n; ++j) printf("%lld ", ans.sum.a[i][j]); printf("\n"); } return 0; }
