矩阵 模板（矩阵加减，矩阵乘法，矩阵求逆）

struct matrix
{
    int m[100][100];
    int w,h;
    matrix(int x,int y)
    {
        h=x,w=y;
        memset(m,0,sizeof(m));
    }
    matrix()
    {
        memset(m,0,sizeof(m));
    }
    void in()
    {
        for(int i=0;i<h;++i)
            for(int j=0;j<w;++j)
                scanf("%d",&m[i][j]);
        return ;
    }
    void print()
    {
        for(int i=0;i<h;++i)
        {
            printf("%d",(m[i][0]%mod+mod)%mod);
            for(int j=1;j<w;++j)
                printf(" %d",((m[i][j]%mod)+mod)%mod);
            printf("\n");
        }
        return ;
    }
    matrix operator +(const matrix &s)const
    {
        matrix ret(h,w);
        for(int i=0;i<h;++i)
            for(int j=0;j<w;++j)
                ret.m[i][j]=((long long)m[i][j]+s.m[i][j])%mod;
        return ret;
    }
    matrix operator -(const matrix &s)const
    {
        matrix ret(h,w);
        for(int i=0;i<h;++i)
            for(int j=0;j<w;++j)
                ret.m[i][j]=((long long)m[i][j]-s.m[i][j])%mod;
        return ret;
    }
    matrix operator *(const matrix &s)const
    {
        matrix ret(h,w);
        for(int i=0;i<h;++i)
            for(int j=0;j<s.w;++j)
                for(int k=0;k<w;++k)
                    ret.m[i][j]=(((long long)m[i][k]*s.m[k][j])%mod+(long long)ret.m[i][j])%mod;
        return ret;
    }
    matrix env()
    {
        matrix ret(h,w);
        for(int i=0;i<h;++i)
            ret.m[i][i]=1;
        for(int i=0;i<h;++i)
        {
            int j;
            for(j=i;j<h&&!m[j][i];++j);
            if(j>=h)
                continue;
            for(int k=0;k<w;++k)
                swap(m[i][k],m[j][k]),
                swap(ret.m[i][k],ret.m[j][k]);
            long long tmp=m[i][i];
            for(j=0;j<w;++j)
                m[i][j]=((long long)m[i][j]*pow(tmp,mod-2))%mod,
                ret.m[i][j]=((long long)ret.m[i][j]*pow(tmp,mod-2))%mod;
            for(j=i+1;j<h;++j)
            {
                long long tp=((long long)m[j][i]*pow(m[i][i],mod-2))%mod;
                for(int k=0;k<w;++k)
                {
                    m[j][k]=((long long)m[j][k]-((long long)m[i][k]*tp)%mod)%mod,
                    ret.m[j][k]=((long long)ret.m[j][k]-((long long)ret.m[i][k]*tp)%mod)%mod;
                }
            }
        }
        for(int i=h-1;i>=0;--i)
            for(int j=i-1;j>=0;--j)
            {
                long long tmp=((long long)m[j][i]*pow(m[i][i],mod-2))%mod;
                m[j][i]=((long long)m[j][i]-((long long)m[i][i]*tmp)%mod)%mod;
                for(int k=0;k<w;++k)
                    ret.m[j][k]=((long long)ret.m[j][k]-((long long)ret.m[i][k]*tmp)%mod)%mod;
            }
        return ret;
    }
    matrix operator ^(const long long &s)const
    {
        matrix ret(h,w);
        for(int i=0;i<h;++i)
            ret.m[i][i]=1;
        long long tmp=s;
        matrix zz;
        zz=(*this);
        while(tmp)
        {
            if(tmp&1)
                ret=ret*zz;
            zz=zz*zz;
            tmp>>=1;
        }
        return ret;
    }
};