矩阵与线性代数学习笔记

矩阵简介

矩阵(Matrix)可以看作一个二维数组，就是一个有n行m列的填满数的数字阵。一个矩阵可以写作：

A=⎛⎝⎜⎜⎜⎜a11a21⋮an1a12a22⋮an2…………a1ma2m⋮anm⎞⎠⎟⎟⎟⎟

关于矩阵论的最基础理论可以参考1，这本书从需要的基础知识讲起，通过线性方程和向量空间、线性无关组、包和维度、矩阵的秩等概念逐步深入矩阵理论，之后引出行列式的计算与构造。在第三章还给出了近世代数的一些基本概念。对与OI来说，这本书的内容（从知识上）应该已经足够。

求解线性方程

通过直观的高斯消元法可以在O(n3)的时间复杂度内求解含有n个方程的方程组。本质的说，这是一个解空间的过程，也就是求解标准基下的一个向量再矩阵所代表的基下的解。也就是说，我们可以通过方程组建立模型来解决一些问题。

bzoj1013球形空间产生器

题意：给你n+1个n维空间中的点，求他们所在的球的球心。

我们大概都还记得平面上不同线的三点确定一个圆。这是为什么呢？设三点坐标为(x0,y0),(x1,y1),(x2,y2)，圆心坐标(x,y)，半径为r，那么有：

⎧⎩⎨⎪⎪⎪⎪(x0−x)2+(y0−y)2=r2(x1−x)2+(y1−y)2=r2(x2−x)2+(y2−y)2=r2

也就是：

⎧⎩⎨⎪⎪⎪⎪x20−2xx0+x2+y20−2yy0+y20=r2x21−2xx1+x2+y21−2yy1+y21=r2x22−2xx2+x2+y22−2yy2+y22=r2

由于r2和x2很不优雅，我们通过(2)−(1),(3)−(2)得到两组新方程，不难发现其为线性方程组。C为常数，这里省略不写。

{2(x0−x1)x+2(y0−y1)y=C12(x0−x2)x+2(y0−y2)y=C1

当且仅当∣∣2(x0−x1)2(y0−y1) 2(x0−x2)2(y0−y2)∣∣≠0时，方程有唯一解。用高斯消元可以O(n3)解出。发现其恰好为三点不共线的条件（从叉积考虑）。

n和n+1只要仿照这里就可以了。

#include <bits/stdc++.h>
using namespace std;

const int N = 20;
double A[N][N];

int n;

void solve()
{
    for (int i = 1; i <= n+1; i++) {
        if (A[i][i] == 0) {
            int k = i;
            while (A[k][i]) k++;
            swap(A[k], A[i]);
        }
        for (int j = i+1; j <= n; j++) {
            for (int k = i+1; k <= n+1; k++)
                A[j][k] -= A[i][k]*(A[j][i]/A[i][i]);
            A[j][i] = 0;
        }
    }
    for (int i = n; i > 1; i--) {
        for (int j = 1; j < i; j++)
            A[j][n+1] -= A[i][n+1]*(A[j][i]/A[i][i]), A[j][i] = 0;
    }
    for (int i = 1; i < n; i++)
        printf("%.3f ", A[i][n+1]/A[i][i]);
    printf("%.3f", A[n][n+1]/A[n][n]);
}

int main()
{
    scanf("%d", &n);
    double a[N], d[N];
    for (int i = 1; i <= n; i++) scanf("%lf", &a[i]);
    for (int i = 2; i <= n+1; i++) {
        double ans = 0;
        for (int j = 1; j <= n; j++)
            scanf("%lf", &d[j]), A[i-1][j] = 2*(d[j]-a[j]), ans += d[j]*d[j]-a[j]*a[j];
        A[i-1][n+1] = ans;
    }
    solve();
    return 0;
}

求解带环dp

当dp方程并不满足最优子结构，即构成的决策图上有环时，如果决策会收敛（比如概率的指数衰减可以收敛），那么可以通过列方程求解。

HNOI2013 游走

一个无向连通图，顶点从1编号到N，边从1编号到M。
小Z在该图上进行随机游走，初始时小Z在1号顶点，每一步小Z以相等的概率随机选 择当前顶点的某条边，沿着这条边走到下一个顶点，获得等于这条边的编号的分数。当小Z 到达N号顶点时游走结束，总分为所有获得的分数之和。
现在，请你对这M条边进行编号，使得小Z获得的总分的期望值最小。
N≤500

思路与解： 我们考虑用一个转移方程描述这个问题。设f(i)为从1到i的期望走过每一条边次数构成的向量(l(E1),l(E2),…,l(Em))，则方程为：

f(1)=(0,0,0,…,0)f(i)=∑j→i(…,(j→i)+1dj,…)

解出该方程f(n)之后贪心的编号即可。然而我们的矩阵方程是n×n2的，可以被卡到O(n4)。一个解决方法是将f写在矩阵一边，l写在另一边，如果用快速矩阵乘法的奥义可以优化到O(n3+T(n×n,n×n2))，应该可以通过此题！

然而我们显然不能随手写出O(nω)的矩阵乘法…因此考虑更换思路，设h(i)为期望经过节点i的次数，则可以列出以下方程：

h(1)=1h(n)=0h(i)=∑j→i1djh(j)

h(n)=0是为了防止对其他节点产生影响。求解完之后每条边i→j的期望经过次数就是1dih(i)+1djh(j)，复杂度成功的降到了O(n3)。

返回来看我们一开始的思路，其有大量优化余地的原因是右侧矩阵虽然是O(n×n2)的，但最多只有O(n2)个元素。因此我们可以不再维护和这个点对无关的边的信息，转而维护点的性质，从而解决了这个问题。

#include <bits/stdc++.h>
using namespace std;

const int N = 505;
int g[N][N], n, m;
double A[N][N];
const double eps = 1e-9;

void solve(double ans[])
{
    for (int i = 1; i < n; i++) {
        int pos = i;
        for (int k = i+1; k <= n; k++)
            if (abs(A[k][i]) > abs(A[pos][i]))
                pos = k;
        if (pos != i)
            swap(A[i], A[pos]);
        for (int j = i+1; j <= n; j++) {
            for (int k = i+1; k <= n+1; k++)
                A[j][k] -= A[i][k]*(A[j][i]/A[i][i]);
            A[j][i] = 0;
        }
    }
    for (int i = n; i >= 1; i--) {
        for (int j = 1; j < i; j++)
            A[j][n+1] -= A[i][n+1]*(A[j][i]/A[i][i]), A[j][i] = 0;
        ans[i] = A[i][n+1]/A[i][i];
    }
}

double ans[N];
double E[N*N];
int top = 0;

int main()
{
    scanf("%d%d", &n, &m);
    for (int i = 1; i <= m; i++) {
        int u, v; scanf("%d%d", &u, &v);
        g[u][v] = g[v][u] = 1;
        g[u][u]++, g[v][v]++;
    }
    for (int i = 1; i < n; i++) {
        A[i][i] = -1;
        for (int j = 1; j <= n; j++)
            if (g[i][j] && i != j)
                A[i][j] = 1.0/g[j][j];
    }
    A[n][n] = 1, A[n][n+1] = 0, A[1][n+1] = -1;
    solve(ans);
    ans[n] = 1;
    for (int i = 1; i < n; i++) {
        for (int j = 1; j < i; j++)
            if (g[i][j])
                E[++top] = ans[i]*1.0/g[i][i]+ans[j]*1.0/g[j][j];
    }
    for (int i = 1; i < n; i++)
        if (g[i][n])
            E[++top] = ans[i]*1.0/g[i][i];
    sort(E+1, E+m+1);
    double a = 0;
    for (int i = 1; i <= m; i++)
        a += (m-i+1)*E[i];
    printf("%.3f", a);
    return 0;
}

Usaco2010 Hot Dotp驱赶猪猡

给定一个无向图，炸弹从1号点出发，每个时刻有P/Q的概率爆炸，如果没有爆炸，就会等概率沿着随机一条出边走到下一个城市，求最终每个城市的爆炸概率。
n≤300

思路与解： 我们设f(i)为i时刻炸弹停留再每个城市概率形成的向量，容易构造一个n×n的转移矩阵A，使得

Af(i)=f(i+1)

我们考虑最终答案ans，有

ans=PQ∑i=0∞Aif(1)

由于我们知道答案随i指数收敛，因此必然有

∑i=0∞Ai=II−A

也就是无限递减几何级数。

然后只要求解矩阵方程：

(I−A)×ans=PQf(1)

嘴巴AC而已

行列式与Matrix_Tree定理

行列式是个奇怪的东西。它可以用许多方式定义，也有着丰富的应用，并且从每个应用都能引出对其的一种定义。矩阵树定理告诉我们：一个图的生成树个数等于其度数矩阵减去邻接矩阵，关于任一一个元素的余子式的绝对值。

bzoj4031 小Z的房间

Description

你突然有了一个大房子，房子里面有一些房间。事实上，你的房子可以看做是一个包含n*m个格子的格状矩形，每个格子是一个房间或者是一个柱子。在一开始的时候，相邻的格子之间都有墙隔着。
你想要打通一些相邻房间的墙，使得所有房间能够互相到达。在此过程中，你不能把房子给打穿，或者打通柱子（以及柱子旁边的墙）。同时，你不希望在房子中有小偷的时候会很难抓，所以你希望任意两个房间之间都只有一条通路。现在，你希望统计一共有多少种可行的方案。

Input

第一行两个数分别表示n和m。
接下来n行，每行m个字符，每个字符都会是.或者*，其中.代表房间，*代表柱子。

Output

一行一个整数，表示合法的方案数 Mod 10^9

思路与解： 裸题。随手A？？

并不..注意到模数并不是素数，不能简单的用逆元表示除法。
我们需要一个被称为“模意义下行列式”的神奇方法。具体来说就是利用辗转相除的思路消元。辗转相除的过程中我们可以把a,b→b,a mod b最终将b消成0，这里只需要模拟这个过程就可以了。可以参考2。

注意在消元的过程中只能严格的用两种初等变换，即

• 交换两行
• 将一行乘以常数，加到另一行

// heoi2015 小z的房间
#include <bits/stdc++.h>
using namespace std;

int p[10][10], num[15][15], d[一个101], top = 0;
long long g[101][101];
int n, m;
char str[10];
int dx[] = {1, 0, -1, 0}, dy[] = {0, 1, 0, -1};
long long mod = 1e9;

long long Gauss(long long M[101][101], int n)
{
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= n; j++)
            ((M[i][j] %= mod) += mod ) %= mod;
    int dsgn = 1;
    for (int i = 1; i <= n; i++) {
        if (M[i][i] == 0) {
            int k = i+1;
            while (k <= n && M[k][i] == 0) k++;
            if (k > n) return 0; // 无解
            swap(M[i], M[k]);
            dsgn *= -1;
        }
        for (int j = i+1; j <= n; j++) {
            if (M[j][i] > M[i][i])
                swap(M[i], M[j]), dsgn *= -1;
            long long A = M[i][i], B = M[j][i];
            while (B) {
                long long t = A/B; A %= B, swap(A, B);
                for (int k = 1; k <= n; k++)
                    ((M[i][k] -= t*M[j][k]) += mod) %= mod;
                swap(M[i], M[j]);
                dsgn *= -1;
            }
        }
    }
    long long ans = dsgn;
    for (int i = 1; i <= n; i++)
        (ans *= M[i][i]) %= mod;
    return (ans += mod) %= mod;
}

int main()
{
    scanf("%d%d", &n, &m);
    for (int i = 1; i <= n; i++) {
        scanf("%s", str+1);
        for (int j = 1; j <= m; j++)
            if (str[j] == '.')
                p[i][j] = 1, num[i][j] = ++top;
            else
                p[i][j] = 0;
    }
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= m; j++)
            for (int k = 0; k < 4; k++)
                if (num[i+dx[k]][j+dy[k]]) {
                    g[num[i][j]][num[i+dx[k]][j+dy[k]]] = -1;
                    g[num[i][j]][num[i][j]]++;
                }
    long long ans = Gauss(g, top-1);
    cout << ans << endl;
    return 0;
}

逆矩阵

我们知道非奇异n阶矩阵和矩阵乘法构成一群，也就是说每个元素都有逆元。考虑如何求解逆矩阵。使用Gauss-Jordan方法，即对矩阵[A|I]运用消元，将A变成单位矩阵，则右侧矩阵就是其逆矩阵。

bzoj4128 Matrix

对于非奇异矩阵A，B和一质数p，求满足

Ax≡B(mod p)

的最小的x。

思路与解：一眼的BSGS。就是设s=p√，我们对于i∈[0,s)暴力计算并丢进哈希表，对于x≥s，它必然可以唯一确定的写成ls+r的形式且r<s(也就是带余除法)。我们只要枚举l，计算Als并查表即可。

然而这题卡常数…我每次循环使用三次乘法果断gg…使用有理有据的底层优化才能AC..

#include <bits/stdc++.h>
using namespace std;

const int MAXN = 75;
int n, mod;
int inv[20005];
int power(int a, int n)
{
    if (n == 0) return 1;
    int p = power(a, n>>1);
    p = p*p%mod;
    if (n&1) p = p*a%mod;
    return p;
}

void read(int A[MAXN][MAXN])
{
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= n; j++)
            scanf("%d", &A[i][j]);
}

void put(int A[MAXN][MAXN])
{
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++)
            printf("%d ", A[i][j]);
        puts("");
    }
}

int M[MAXN][MAXN], A[MAXN][MAXN], B[MAXN][MAXN], ret[MAXN][MAXN];
void inverse(int A[MAXN][MAXN], int ans[MAXN][MAXN])
{
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= n; j++)
            M[i][j] = A[i][j], ans[i][j] = 0;
    for (int i = 1; i <= n; i++) ans[i][i] = 1;
    for (int i = 1; i <= n; i++) {
        int pos = i;
        for (register int k = i; k <= n; k++)
            if (M[k][i] != 0) {
                pos = k; break;
            }
        swap(M[i], M[pos]);
        int iv = inv[M[i][i]];
        for (int j = 1; j <= n; j++)
            M[i][j] = M[i][j]*iv%mod, ans[i][j] = ans[i][j]*iv%mod;
        for (register int j = 1; j <= n; j++) {
            if (j == i) continue;
            int tab = M[j][i];
            for (register int k = 1; k <= n; k++) {
                M[j][k] = ((M[j][k]-M[i][k]*tab)%mod+mod)%mod;
                ans[j][k] = ((ans[j][k]-ans[i][k]*tab)%mod+mod)%mod;
            }
        }
    }
}

void mul(int A[MAXN][MAXN], int B[MAXN][MAXN], int C[MAXN][MAXN])
{
    memset(ret, 0, sizeof ret);
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= n; j++) {
            int mid = n/4, t = mid*4;
            for (register int k = 1; k < t; k += 4) {
                ret[i][j] = (ret[i][j] + A[i][k]*B[k][j])%mod;
                ret[i][j] = (ret[i][j] + A[i][k+1]*B[k+1][j])%mod;
                ret[i][j] = (ret[i][j] + A[i][k+2]*B[k+2][j])%mod;
                ret[i][j] = (ret[i][j] + A[i][k+3]*B[k+3][j])%mod;
            }
            ret[i][j] = (ret[i][j] + A[i][t+1]*B[t+1][j])%mod;
            ret[i][j] = (ret[i][j] + A[i][t+2]*B[t+2][j])%mod;
            ret[i][j] = (ret[i][j] + A[i][t+3]*B[t+3][j])%mod;
        }
    memcpy(C, ret, sizeof ret);
}

unsigned int hash_val(int A[MAXN][MAXN])
{
    unsigned int val = 0;
    for (int i = 1; i <= n; i++)
        val = (val+A[i][i])*mod;
    return val;
}
map<int,int> mp;

int I[MAXN][MAXN], C[MAXN][MAXN], D[MAXN][MAXN], T[MAXN][MAXN];

int main()
{
    scanf("%d%d", &n, &mod);
    for (register int i = 1; i < mod; i++)
        inv[i] = power(i, mod-2);
    for (int i = 1; i <= n; i++)
        I[i][i] = 1;
    read(A);
    read(B);
    int stp = sqrt(mod)+1;
    memcpy(C, I, sizeof I);
    for (int i = 0; i < stp; i++) {
        mp[hash_val(C)] = i;
        mul(C, A, C);
    }
    inverse(C, D);
    memcpy(T, I, sizeof I);
    for (int i = 0; i <= stp; i++) {
        int ht = hash_val(B);
        if (mp[ht]) {
            cout << i*stp+mp[ht] << endl;
            return 0;
        }
        mul(B, D, B);
    }
    return 0;
}

bzoj3168钙铁锌硒维生素

留一个坑。大概思路就是求是否存在一组非零向量使其通过A线性组合出一组向量。放在矩阵里一起处理即可。

线性基

就是对一个自然数集合S找一个最小的的集合R，使得∀x∈S,∃c,Rc1⊕Rc2…Rc|c|=x。

由于异或可以看成模2意义下的加法，线性基其实是S对应向量空间的极大线性无关组。由于矩阵行秩等于列秩，线性基中元素只有O(lgn)个

如何构造一个线性基？只要在线的加入即可。注意逐位考虑。O(lgRmax)

bool add(long long num)
{
    for (int i = 1; i <= tp; i++) {
        if ((num^k[i]) < num && (num^k[i]) < k[i])
            num ^= k[i];
    }
    if (num != 0) {
        k[++tp] = num;
        for (int i = tp; i-1 && k[i-1] < num; i--)
            swap(k[i-1], k[i]);
        return 1;
    }
    return 0;
}

2460: [BeiJing2011]元素

线性基裸题，只要贪心加入即可。
证明？由于不会有不相交的线性无关组，考虑两个序列，容易说明贪心策略不劣于任何策略。

#include <bits/stdc++.h>
using namespace std;

long long k[70], n, tp = 0;
struct p {
    long long num, mag;
} ps[1005];
bool cmp(const p &a, const p &b)
{ return a.mag > b.mag; }
bool add(long long num)
{
    for (int i = 1; i <= tp; i++) {
        if ((num^k[i]) < num && (num^k[i]) < k[i])
            num ^= k[i];
    }
    if (num != 0) {
        k[++tp] = num;
        for (int i = tp; i-1 && k[i-1] < num; i--)
            swap(k[i-1], k[i]);
        return 1;
    }
    return 0;
}

int main()
{
    scanf("%d", &n);
    for (int i = 1; i <= n; i++)
        scanf("%lld%d", &ps[i].num, &ps[i].mag);
    sort(ps+1, ps+n+1, cmp);
    int ans = 0;
    for (int i = 1; i <= n; i++)
        if (add(ps[i].num))
            ans += ps[i].mag;
    cout << ans << endl;
    return 0;
}

[Scoi2016]幸运数字

问树上路径走过的元素形成的集合的最大异或值。

思路与解： 由于所有元素的异或值和线性基的异或值等价，且线性基合并是满足结合律的，我们只要用树上倍增维护线性基合并即可。注意判断a=b。

#include <bits/stdc++.h>
using namespace std;

struct linear_basis {
    long long k[64];
    int top;
    inline void clear() { top = 0; }
    linear_basis() { clear(); }
    void push(long long num)
    {
        for (int i = 1; i <= top; i++)
            if ((k[i]^num) < k[i] && (k[i]^num) < num)
                num ^= k[i];
        if (num) {
            k[++top] = num;
            for (int i = top; i-1 && k[i-1] < k[i]; i--)
                swap(k[i-1], k[i]);
        }
    }
    friend linear_basis operator + (linear_basis a, linear_basis b)
    {
        for (int i = 1; i <= a.top; i++)
            b.push(a.k[i]);
        return b;
    }
    long long max()
    {
        long long ans = k[1];
        for (int i = 2; i <= top; i++)
            if ((ans^k[i]) >= ans)
                ans ^= k[i];
        return ans;
    }
    void print()
    {
        for (int i = 1; i <= top; i++)
            cout << k[i] << " ";
        cout << endl;
    }
} ;

const int MAXN = 20005;
linear_basis f[MAXN][16];
int fa[MAXN][16], depth[MAXN];
long long a[MAXN];
int n, q;
struct node {
    int to, next;
} edge[MAXN*2];
int head[MAXN], top, u, v;
inline void push(int i, int j)
{ ++top, edge[top] = {j, head[i]}, head[i] = top; }
void dfs(int nd, int ft = 0)
{
    fa[nd][0] = ft, f[nd][0].push(a[ft]), depth[nd] = depth[ft]+1;
    for (int i = head[nd]; i; i = edge[i].next)
        if (edge[i].to != ft)
            dfs(edge[i].to, nd);
}

linear_basis lca(int a, int b)
{
    if (depth[a] < depth[b]) swap(a, b);
    linear_basis lb;
    for (int i = 0; i < 16; i++)
        if ((depth[a]-depth[b])&(1<<i))
            lb = lb + f[a][i], a = fa[a][i];
    if (a == b) return lb;
    for (int i = 15; i >= 0; i--)
        if (fa[a][i] != fa[b][i]) {
            lb = lb + f[a][i] + f[b][i];
            a = fa[a][i], b = fa[b][i];
        }
    if (a != b) lb = lb + f[a][0] + f[b][0];
    return lb;
}

int main()
{
    scanf("%d%d", &n, &q);
    for (int i = 1; i <= n; i++) scanf("%lld", &a[i]), f[i][0].push(a[i]);
    for (int i = 1; i < n; i++) {
        scanf("%d%d", &u, &v);
        push(u, v), push(v, u);
    }
    dfs(1);
    for (int j = 1; j < 16; j++)
        for (int i = 1; i <= n; i++)
            fa[i][j] = fa[fa[i][j-1]][j-1];
    for (int j = 1; j < 16; j++)
        for (int i = 1; i <= n; i++)
            f[i][j] = f[i][j-1]+f[fa[i][j-1]][j-1];
    for (int i = 1; i <= q; i++) {
        scanf("%d%d", &u, &v);
        if (u == v) printf("%lld\n", a[u]);
        else printf("%lld\n", lca(u, v).max());
    }
    return 0;
}

线性齐次递推

我们考虑一个递推式f(n)=∑mi=1aif(n−i)，我们希望在更快的速度内完成计算。考虑用矩阵，我们发现：

A=⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜a11100⋮0a12010⋮0a13001⋮0a14000⋮0………………a1m−1000⋮1a1m000⋮0⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟

h(n)=⎛⎝⎜⎜⎜⎜⎜⎜⎜f(n)f(n−1)f(n−2)⋮f(n−m)⎞⎠⎟⎟⎟⎟⎟⎟⎟

则有

Ah(n)=h(n+1)

由于矩阵乘法满足结合律，可以用快速幂计算。

[HNOI2008]GT考试

阿申准备报名参加GT考试，准考证号为N位数X1X2....Xn(0≤Xi≤9),他不希望准考证号上出现不吉利的数字。
他的不吉利数学A1A2...Am(0≤Ai≤9)有M位，不出现是指X1X2...Xn中没有恰好一段等于A1A2...Am. A1和X1可以为0。

思路与解： 一个数位dp的套路，dp[i][j]表示还剩i位，匹配到第j位的方案数，则dp[i][j]=∑j→idp[i−1][k]。转移可以O(n)用kmp搞出来。然后就是矩阵乘起来了。

#include <bits/stdc++.h>
using namespace std;

const int MAXN = 30;
int n, m, p;
char str[30];
int pi[30];
int trans[30][10];
int M[MAXN][MAXN];

void kmp_init()
{
    pi[0] = pi[1] = 0;
    int p = 0;
    for (int i = 2; i <= m; i++) {
        while (p && str[p+1] != str[i]) p = pi[p];
        if (str[p+1] == str[i]) p++;
        pi[i] = p;
    }
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < 10; j++) {
            int k = i;
            while (k && str[k+1] != j+'0') k = pi[k];
            if (str[k+1] == j+'0')
                trans[i][j] = k+1, M[k+1][i]++;
            else
                M[0][i]++;
        }
    }
    M[m][m] = 10;
}

int R[MAXN][MAXN];

void mul(int A[MAXN][MAXN], int B[MAXN][MAXN]) // A *= B
{
    memset(R, 0, sizeof R);
    for (int i = 0; i <= m; i++)
        for (int j = 0; j <= m; j++)
            for (int k = 0; k <= m; k++)
                (R[i][j] += A[i][k]*B[k][j]) %= p;
    for (int i = 0; i <= m; i++)
        for (int j = 0; j <= m; j++)
            A[i][j] = R[i][j];
}

int B[33][MAXN][MAXN], top = 0;

void quick_power(int A[MAXN][MAXN], int n)
{
    if (n == 1) return;
    ++top;
    for (int i = 0; i <= m; i++)
        for (int j = 0; j <= m; j++)
            B[top][i][j] = A[i][j];
    quick_power(B[top], n>>1);
    mul(B[top], B[top]);
    if (n&1) mul(B[top], A);
    for (int i = 0; i <= m; i++)
        for (int j = 0; j <= m; j++)
            A[i][j] = B[top][i][j];
}

int main()
{
    scanf("%d%d%d%s", &n, &m, &p, str+1);
    kmp_init();
    quick_power(M, n);
    int ans = 0;
    for (int i = 0; i < m; i++)
        (ans += M[i][0]) %= p;
    cout << ans << endl;
    return 0;
}

[HNOI2011]数学作业

给定n和m，求：

1234567891011...n mod m

思路与解： 可以不同位数分开考虑。考虑第k位上的情况：

f(i)=10if(i−1)+i

由于要记录i的信息，我们要扩充状态表示：

h(n)=⎛⎝⎜⎜⎜f(i)f(i−1)i+11⎞⎠⎟⎟⎟

可以构造矩阵：

A=⎛⎝⎜⎜⎜10k100000010100011⎞⎠⎟⎟⎟

即可。注意取模的技巧..否则会RE+WA+MLE.

#include <bits/stdc++.h>
using namespace std;

long long n, m;
long long D[5][5];

void mul(long long A[5][5], long long B[5][5], long long C[5][5])
{
    memset(D, 0, sizeof D);
    for (int i = 1; i <= 4; i++)
        for (int j = 1; j <= 4; j++)
            for (int k = 1; k <= 4; k++)
                (D[i][j] += A[i][k]*B[k][j]%m) %= m;
    for (int i = 1; i <= 4; i++)
        for (int j = 1; j <= 4; j++) {
            C[i][j] = D[i][j];
            if (C[i][j] < 0) throw;
        }
}

long long B[5][5];

void matrix_power(long long A[5][5], long long n)
{
    memset(B, 0, sizeof B);
    for (int i = 1; i <= 5; i++) B[i][i] = 1;
    while (n) {
        if (n&1) mul(B, A, B);
        mul(A, A, A); n >>= 1;
    }
    for (int i = 1; i <= 5; i++)
        for (int j = 1; j <= 5; j++)
            A[i][j] = B[i][j];
}

long long power(long long a, long long n, long long mod = 3e18)
{
    long long t = 1;
    while (n) {
        if (n&1) t = t*a%mod;
        a = a*a%mod;
        n >>= 1;
    }
    return t;
}

long long M[5][5];

int main()
{
    scanf("%lld%lld", &n, &m);
    long long ret = 0;
    for (int k = 1; k <= 19; k++) {
        if (power(10, k-1) > n) break;
        long long beg = power(10, k-1), tar = min(power(10, k)-1, n);
        memset(M, 0, sizeof M);
        M[1][1] = power(10, k, m), M[1][3] = 1;
        M[2][1] = M[3][3] = M[3][4] = M[4][4] = 1;
        matrix_power(M, tar-beg+1);
        ret = ret*power(power(10, (tar-beg+1), m)%m, k, m)%m;
        ret = (ret+M[1][3]*(beg%m)%m+M[1][4])%m;
    }
    cout << ret << endl;
    return 0;
}

// bug记录，函数默认参数忘记下传gg...
// 尽量不用默认参数...

参考资料3

---

1. 《代数学引论I(BA I)》，前苏联-柯斯特利金著，张英伯译，高等教育出版社 ↩
2. http://blog.csdn.net/u013010295/article/details/47451451 ↩
3. 某集训，妹主席ppt ↩