LOJ6440 万能欧几里得

Link
我们可以这么看：
初始：(a=I,b=I,ans=0)。
(g(1):b=b*B)
(g(0):a=a*A,ansleftarrow ans+a*b)

#include<cstdio>
#include<cstring>
using i64=long long;
const int N=23,P=998244353;
int n;
i64 read(){i64 x;scanf("%lld",&x);return x;}
void inc(int&a,int b){a+=b-P,a+=a>>31&P;}
int add(int a,int b){return inc(a,b),a;}
int mul(int a,int b){return 1ll*a*b%P;}
struct matrix{int a[N][N];matrix(){memset(a,0,sizeof a);}}I;
matrix operator*(matrix a,matrix b)
{
    matrix c;
    for(int i=0;i<n;++i) for(int j=0;j<n;++j) for(int k=0;k<n;++k) inc(c.a[i][j],mul(a.a[i][k],b.a[k][j]));
    return c;
}
matrix operator+(matrix a,matrix b)
{
    matrix c;
    for(int i=0;i<n;++i) for(int j=0;j<n;++j) c.a[i][j]=add(a.a[i][j],b.a[i][j]);
    return c;
}
struct node{matrix a,b,s;node(){a=b=I;}};
node operator+(node a,node b){node c;c.a=a.a*b.a,c.b=a.b*b.b,c.s=a.s+a.a*b.s*a.b;return c;}
node operator*(node a,i64 k){node r;for(;k;k>>=1,a=a+a) if(k&1) r=r+a;return r;}
node f(i64 a,i64 b,i64 c,i64 n,node p,node q)
{
    if(!n) return node();
    if(!a) return (q*(c/b)+p)*n;
    c%=b;
    if(a>=b) return f(a%b,b,c,n,q*(a/b)+p,q);
    i64 k=((__int128)a*n+c)/b;
    if(!k) return p*n;
    return p*((b-c-1)/a)+q+f(b,a,b-c-1,k-1,q,p)+p*(n-((__int128)k*b-c-1)/a);
}

int main()
{
    i64 a=read(),b=read(),c=read(),m=read();n=read();
    matrix A,B;
    for(int i=0;i<n;++i) I.a[i][i]=1;
    for(int i=0;i<n;++i) for(int j=0;j<n;++j) A.a[i][j]=read();
    for(int i=0;i<n;++i) for(int j=0;j<n;++j) B.a[i][j]=read();
    node p,q;p.a=p.s=A,q.b=B;
    matrix ans=(q*(c/b)+f(a,b,c,m,p,q)).s;
    for(int i=0;i<n;++i,puts("")) for(int j=0;j<n;++j) printf("%d ",ans.a[i][j]);
}