题目
解析
以坐标原点为起点,将矩阵$A$看作向右走,$B$看作向上走,平面上的点$(x, y)$有贡献$A^xB^y$,答案就是$y=left lfloor frac{Px+R}{Q} ight floor(1 leqslant x leqslant L)$每一个$x$下第一个整点的贡献和。
来一发课件,下图中斜线段即为$y=left lfloor frac{Px+R}{Q} ight floor$
操作序列满⾜信息可合并。只需要两个操作序列的信息,即可计算出将这两个操作序列拼接后的信息。
操作可以被看成代码段:
初始化: 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$,编号为$1sim L$。第$i-1$个$B$与第$i$个$B$之间有$left lfloor frac{Pi+R}{Q} ight floor-left lfloor frac{P(i-1)+R}{Q} ight floor$个$A$,也就是说这里是以原点为起点,将操作$B$看作向右走,操作$A$看作向上走,斜线为$y=left lfloor frac{Pi+R}{Q} ight floor$。
根据这个函数的定义,有如下两个性质:
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}
ight
floor \ Qx leqslant Py+R \ left lceil frac{Qx-R}{P}
ight
ceil leqslant y \ left lfloor frac{Qx-R+P-1}{P}
ight
floor leqslant y $$
也就是说第$x$个$A$前有$left lfloor frac{Qx-R-1}{P}
ight
floor$个$B$,第$x-1$个$A$与第$x$个$A$之间有$left lfloor frac{Qx-R-1}{P}
ight
floor-left lfloor frac{Q(x-1)-R-1}{P}
ight
floor$个$B$,似乎可以递归了。但由于当$x$等于$1$时,$left lfloor frac{Q(x-1)-R-1}{P}
ight
floor < 0$,因此第一个$A$需要单独计算,假设有$m$个$A$,这样的话还需要计算的$x$的范围为$2leqslant x leqslant m$,但$solve$函数的$x$是从$1$开始枚举的,所以用$x-1$代替原来的$x$,于是就变成了,第$x-1$个$A$与第$x$个$A$之间有$left lfloor frac{Qx+Q-R-1}{P}
ight
floor-left lfloor frac{Q(x-1)+Q-R-1}{P}
ight
floor$个$B$,此时$x$的范围为$1leqslant 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(" "); } return 0; }