矩阵快速幂【模板】

const int MAXN = 110;
struct Matrax
{
    int m[MAXN][MAXN];
}a,per;

int N,M;

void Init()
{
    for(int i = 0; i < N; ++i)
        for(int j = 0; j < N; ++j)
        {
            scanf("%d",&a.m[i][j]);
            a.m[i][j] %= M;
            per.m[i][j] = (i == j);
        }
}

Matrax Multi(Matrax a,Matrax b) //a * b
{
    Matrax c;
    for(int i = 0; i < N; ++i)
        for(int j = 0; j < N; ++j)
        {
            c.m[i][j] = 0;
            for(int k = 0; k < N; ++k)
                c.m[i][j] += a.m[i][k]*b.m[k][j];
            c.m[i][j] %= M;
        }
    return c;
}

Matrax Power(int k) //a^k % M
{
    Matrax c,p,ans = per;
    p = a;
    while(k)
    {
        if(k&1)
            ans = Multi(ans,p);
        p = Multi(p,p);
        k >>= 1;
    }
    return ans;
}

Matrax Add(Matrax a,Matrax b)   // a + b
{
    Matrax c;
    for(int i = 0; i < N; ++i)
        for(int j = 0; j < N; ++j)
            c.m[i][j] = (a.m[i][j] + b.m[i][j]) % M;
    return c;
}
Matrax MatraxSum(int k) //a + a^2 + a^3 + … + a^k
{
    if(k == 1)
        return a;
    Matrax temp,b;
    temp = MatraxSum(k/2);
    if(k&1)
    {
        b = Power(k/2+1);
        temp = Add(temp,Multi(temp,b));
        temp = Add(temp,b);
    }
    else
    {
        b = Power(k/2);
        temp = Add(temp,Multi(temp,b));
    }
    return temp;
}