关于矩阵快速幂的若干优化

首先，我们复习一下矩阵乘法。

我们记3个矩阵A（a行b列）,B（b行c列）,C（a行c列）。我们要计算A*B，并把答案存到矩阵C中。

C[i][j]+=A[i][k]*B[k][j]（1<=i<=a，1<=j<=c，1<=k<=b），即新矩阵的第i行第j个元素是原1矩阵的第i行*原2矩阵的第j列得来的。

一般来说，我们的计算方法是for(int i=1;i<=a;i++)for(int j=1;j<=c;j++)for(int k=1;k<=b;k++)C[i][j]+=A[i][k]*B[k][j];

其次，让我们复习一下快速幂。

举个例子吧，计算a^101。

我们知道，于是：

a^101=a^(1*2^6)*a^(1*2^5)*a^(0*2^4)*a^(1*2^3)*a^(1*2^2)*a^(0*2^1)*a^(1*2^0)。

我们把101转成2进制：1101101。每个2^x前的系数就是二进制第x位的数。

a^ (2^x)=a^(2^(x-1))^2。我们可以通过a^(2^(x-1))来求得a^(2^x)。

这样，对于二进制下的第x位，该位如果为1，就把ans*=a（更新答案，初始化为1）。然后每次a*=a（用a^(2^x)更新出a^(2^(x+1))，准备处理下一位）。

我们便可以在O(logp)(p为指数)的时间复杂度内出解。

最后，让我们来复习一下矩阵快速幂。

我们要求A^B^B^B^B^B^B^B......（A，B为矩阵），即A^(B^p)的值。

就像ans初值=1一样，记一个单位矩阵（主对角线为1）Ans，结合上面两种做法，我们就可以求出A^(B^p)的值。

(1)对于稀疏矩阵的优化

稀疏矩阵，即为矩阵中有很多元素为0。

优化方法：改变循环顺序。改为for(int i=1;i<=a;i++)for(int k=1;k<=b;k++)for(int j=1;j<=c;j++)C[i][j]+=A[i][k]*B[k][j];

这样有什么好处呢？

我们可以发现，只要A[i][k]==0，那么对答案矩阵（C）不会有任何贡献。

所以我们可以进行优化，在第二个循环到第三个循环直接加一个if，若A[i][k]!=0，才进入第三个循环。

for(int i=1;i<=a;i++)for(int k=1;k<=b;k++)if(A[i][k])for(int j=1;j<=c;j++)C[i][j]+=A[i][k]*B[k][j];

题目：POJ 3735 Training little cats。

(2)预处理优化矩阵快速幂

主要针对多组数据。求A*B^k。给出A，B，T个询问k

在通常情况下，A是一个n行1列的矩阵，B是一个n行n列的矩阵。这样，我们的矩阵快速幂（求A^(B^k)）的复杂度就是O((n^3logk+Tn^2logk))。

具体来说，我们先用O(n^3log(maxk))预处理出B^(2^p)，再A*B^k=A*B^(2^a1)+A*B^(2^a2)+...算答案。复杂度O(Tn^2logk)

(3)优化快速幂过程

主要针对多组数据。

正常的快速幂的当次复杂度为O(log2(n))。看到那个2了吗，我们的工作就是要把这个2变大。

考虑一般的快速幂，一般的快速幂是以2进制为基础的，我们考虑用3进制为基础会怎么样。

对于每一个3进制位，如果该位是0，ans*=x^0，如果该位是1，ans*=x^1，如果该位是2，ans*=x^2

与2进制快速幂同理，每次x=x^3，p=p/3

所以复杂度是O(klogk(n))，k为进制

但是虽然这个2变大了，复杂度却一点也没变小

但是这并不能阻挡我们优化的决心，如果每次询问的底数都相同，我们是能优化的

预处理mi[a][b]=(x^(k^a))^b即可，每次ans*=mi[a][b]，a是当前做到第几位，b是当前这位的数

mi[a][1]=mi[a-1][k-1]*mi[a-1][1]

mi[a][b]=mi[a][b-1]*mi[a][1]

这样复杂度变为(klogk(n)+logk(n))。

(4)常数优化

ikj循环，循环展开 for(int i = 1; i <= n; i++) for(int k = 1; k <= n; k++) for(int j = 1; j <= n; j++) c[i][j] += a[i][k] * b[k][j];

这样能保证b数组的内存访问是连续的

 

拥有上面全部优化的模板题：https://www.luogu.org/problemnew/show/P5107

#include <cstdio>
#include <cstring>
#define mod 998244353
#define T 256
#include <algorithm>

struct xxx{
    int a[52][52];
};
struct xx{
    int a[52];
};
struct QQ{
    int x, id;
}q[50100];
int n, d[55];
xxx mi[4][T + 1];
xx ans;
long long Ans[50100];

bool cmp(QQ a, QQ b) {return a.x < b.x;}
 
xxx operator * (xxx a, xxx b)
{
    xxx c; memset(c.a, 0, sizeof(c.a));
    for(int i = 1; i <= n; i++)
        for(int k = 1; k <= n; k++)
            if(a.a[i][k])
            for(int j = 1; j <= n; j++)
                c.a[i][j] = (c.a[i][j] + 1ll * a.a[i][k] * b.a[k][j]) % mod;
    return c;
}

xx operator * (xx a, xxx b) 
{
    xx c; memset(c.a, 0, sizeof(c.a));
    for(int j = 1; j <= n; j++)
        for(int i = 1; i <= n; i++)
            c.a[j] = (c.a[j] + 1ll * a.a[i] * b.a[i][j]) % mod;
    return c;
}

int qpow(int x, int p)
{
    int ans = 1;
    while(p)
    {
        if(p & 1) ans = 1ll * ans * x % mod;
        x = 1ll * x * x % mod; p >>= 1;
    }
    return ans;
}

xx operator ^ (xx a, int p)
{
    int j = 0;
    while(p)
    {
        ans = ans * mi[j][p & 255];
        j++; p >>= 8;
    }
    return ans;
}

int main()
{
    int m, Q; scanf("%d%d%d", &n, &m, &Q);
    for(int i = 1; i <= n; i++) scanf("%d", &ans.a[i]), mi[0][1].a[i][i] = 1, d[i] = 1;
    for(int i = 1; i <= m; i++)
    {
        int u, v; scanf("%d%d", &u, &v);
        mi[0][1].a[u][v]++; d[u]++;
    }
    for(int i = 1; i <= n; i++)
        for(int j = 1; j <= n; j++)
            mi[0][1].a[i][j] = 1ll * mi[0][1].a[i][j] * qpow(d[i], mod - 2) % mod;
    for(int i = 0; i <= 3; i++)
    {
        for(int j = 0; j < T; j++)
        {
            if(i == 0 && j == 1) continue;
            if(j == 0) for(int k = 1; k <= n; k++) mi[i][j].a[k][k] = 1;
            else if(j == 1) mi[i][j] = mi[i - 1][T - 1] * mi[i - 1][1];
            else mi[i][j] = mi[i][j - 1] * mi[i][1];
        }
    }
    for(int i = 1; i <= Q; i++)
    {
        scanf("%d", &q[i].x);
        q[i].id = i;
    }
    std::sort(q + 1, q + Q + 1, cmp);
    for(int i = 1; i <= Q; i++)
    {
        ans = ans ^ (q[i].x - q[i - 1].x);
        for(int j = 1; j <= n; j++) Ans[q[i].id] = Ans[q[i].id] ^ ans.a[j];
        Ans[q[i].id] %= mod;
    }
    for(int i = 1; i <= Q; i++) printf("%lldn", Ans[i]);
}