BJ模拟: 装饰地板（矩阵快速幂+拉格朗日插值+特征多项式）

题意：
给一个$6*R$的矩阵，有两种瓷砖（$1*2$,$2*1$），贡献分别为$s_1,s_2$，一种铺满矩阵的方案的贡献为所有瓷砖的贡献乘积之和，所有方案的贡献和（$\log_2 R \le 8e3$）。

题解：
非常套路的题，裸的矩乘+状态压缩应该都会。

不过因为$\log R$很大，不能直接状压， 这时候可以先$n^4$把矩阵的特征多项式插值出来，然后直接取模即可。

不过细节还是很多的，而且因为全用vector封装导致常数很大，所以跑起来没有优化过常数的矩阵快，不过当做特征多项式练手的题还是不错的。

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int mod=998244353;
inline int add(int x,int y) {return (x+y>=mod)?(x+y-mod):x+y;}
inline int mul(int x,int y) {return (LL)x*y%mod;}
inline int power(int a,int b) {
    int rs=1;
    for(;b;b>>=1,a=mul(a,a)) if(b&1) rs=mul(rs,a);
    return rs;
}
int m,s1,s2,len,lim,bin[10050],pw[10050],yc[10050],g[10050];
struct Bignum {
    int a[10050],len;
}l,r;
istream& operator >>(istream& a,Bignum &b) {
    string s; a>>s; b.len=s.length();
    for(int i=0;i<b.len;i++) b.a[b.len-i-1]+=s[i]-'0';
    for(int i=0;i<b.len;i++) b.a[i+1]+=b.a[i]/10, b.a[i]%=10;
    while(b.a[b.len]) {
        b.a[b.len+1]+=b.a[b.len]/10; 
        b.a[b.len]%=10; 
        ++b.len;
    } --b.len;
    return a;
}
struct matrix {
    int a[65][65];
    matrix(int t=0) {
        memset(a,0,sizeof(a));
        for(int i=0;i<=len;i++) a[i][i]=t;
    }
    friend inline matrix operator *(const matrix &a,const matrix &b) {
        matrix c;
        for(int i=0;i<=len;i++) 
            for(int k=0;k<=len;k++)
                for(int j=0;j<=len;j++) 
                    c.a[i][j]=add(c.a[i][j],mul(a.a[i][k],b.a[k][j]));
        return c;
    }
    friend inline matrix operator *(const matrix &a,const int &b) {
        matrix c=a;
        for(int i=0;i<=len;i++) 
            for(int j=0;j<=len;j++)
                c.a[i][j]=mul(c.a[i][j],b);
        return c;
    }
    inline int det() {
        int sgn=1;
        for(int i=0;i<=len;i++) {
            int p;
            for(p=i;!a[p][i] && p<=len;++p);
            if(p==len+1) return 0;
            if(p!=i) {
                for(int j=i;j<=len;j++) swap(a[i][j],a[p][j]);
                sgn=mod-sgn;
            }
            const int inv=power(a[i][i],mod-2);
            for(int j=i+1;j<=len;j++) {
                int t=mul(a[j][i],inv);
                for(int k=i;k<=len;k++) a[j][k]=add(a[j][k],mod-mul(t,a[i][k]));
            }
        }
        int rs=1;
        for(int i=0;i<=len;i++) rs=mul(rs,a[i][i]);
        return mul(rs,sgn);  
    }
}st;
struct poly {
    int deg;
    vector <int> a;
    poly(int d=0,int t=0) {
        deg=d; a.resize(d+1); a[0]=t;
    }
    inline void init() {
        a.assign(lim,0); deg=lim-1;
        for(int i=1;i<=lim;i++) {
            memset(g,0,sizeof(g)); g[0]=1;
            for(int j=1;j<=lim;j++) {
                if(i==j) continue;
                const int inv=power(add(i,mod-j),mod-2);
                for(int k=lim-1;~k;k--) {
                    g[k]=mul(g[k],mul(mod-j,inv));
                    if(k) g[k]=add(g[k],mul(g[k-1],inv));
                }
            }
            for(int j=0;j<lim;j++) a[j]=add(a[j],mul(g[j],yc[i]));
        }
    }
    inline poly poly_mul(int val) {
        poly c=*this;
        for(int i=0;i<=c.deg;i++) c.a[i]=mul(c.a[i],val);
        return c;
    }
    inline void poly_add(const poly &b) {
        deg=max(deg,b.deg); a.resize(deg+1); 
        for(int i=0;i<=b.deg;i++) a[i]=add(a[i],b.a[i]);
    }
    inline int getval(int x) {
        int rs=1, ans=0;
        for(int i=0;i<=deg;i++)
            ans=add(ans,mul(rs,a[i])), rs=mul(rs,x);
        return ans;
    }
}st2;
inline poly operator %(poly a,const poly &b) {
    int deg=min(a.deg,b.deg-1);
    const int inv=power(b.a[b.deg],mod-2);
    for(int i=a.deg;i>=b.deg;i--) {
        int t=mul(a.a[i],inv);
        for(int j=i,k=b.deg;~k;--j,--k) 
            a.a[j]=add(a.a[j],mod-mul(b.a[k],t));
    }
    a.a.resize(deg+1); a.deg=deg; return a;
}
inline poly mul(const poly &a,const poly &b) {
    int deg=a.deg+b.deg; poly c(deg,0);
    for(int i=0;i<=deg;i++) 
        for(int j=0;j<=a.deg && j<=i;j++) {
            if(i-j<=b.deg) c.a[i]=add(c.a[i],mul(a.a[j],b.a[i-j]));
        }
    return c;
}
inline poly power(poly a,int b) {
    poly rs(0,1);
    for(;b;b>>=1, a=mul(a,a)%st2) 
        if(b&1) rs=mul(rs,a)%st2;
    return rs;
}
inline void dfs(int sta,int p,int s,int ori) {
    if(p==m) {st.a[sta][ori]=s; return;}
    if(sta&bin[p]) dfs(sta,p+1,s,ori);
    else {
        if(p!=m-1 && !(sta&bin[p+1])) dfs(sta|bin[p]|bin[p+1],p+2,mul(s,s2),ori);
        dfs(sta,p+1,s,ori);
    }
}
inline void pre() {
    memset(st.a,0,sizeof(st.a));
    len=1<<m; st.a[len][len-1]=1; st.a[len][len]=1; pw[0]=1; lim=len+10;
    for(int i=0;i<=m;i++) bin[i]=(1<<i)%mod;
    for(int i=1;i<=m;i++) pw[i]=mul(pw[i-1],s1);
    for(int i=0;i<len;i++) {
        int s=(len-1)^i, t=pw[__builtin_popcount(s)];
        dfs(s,0,t,i);
    }
    for(int i=1;i<=lim+10;i++) {
        matrix c;
        for(int j=0;j<=len;j++)
            for(int k=0;k<=len;k++)
                c.a[j][k]=mod-st.a[j][k];
        for(int j=0;j<=len;j++) c.a[j][j]=add(c.a[j][j],i);
        yc[i]=c.det();
    }
    st2.init(); 
    while(!st2.a.back()) st2.a.pop_back(), --st2.deg;
}
inline int solve(Bignum R) {
    poly t(1,0), f(0,1); t.a[1]=1;
    for(int j=0;j<=R.len;j++) {
        f=mul(f,power(t,R.a[j]))%st2;
        t=power(t,10);
    }
    matrix A; A.a[len-1][1]=1; int ans=0;
    for(int i=0;i<=f.deg;i++) {
        ans=add(ans,mul(f.a[i],A.a[len][1]));
        A=st*A; 
    }
    return ans;
}
int main() {
    r.a[0]=1;
    cin>>l>>r>>m>>s1>>s2;
    pre();
    cout<<add(solve(r),mod-solve(l));
}