万能欧几里得知识整理

题目

　　LOJ#6440

解析

　　以坐标原点为起点，将矩阵$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;
}