作者: KSkun

KS的OI模板库开始维护啦~

KS的OI模板库开始维护啦~

最近是省选复习阶段,就想边打板子边整理板子,于是就在GayHub上建了个仓库放模板,也分享 

Lucas定理及其扩展原理及应用

Lucas定理及其扩展原理及应用

Lucas定理(Lucas’s theorem) 内容 对于非负整数m和n和素 

[SCOI2012]喵星球上的点名 题解

[SCOI2012]喵星球上的点名 题解

题目地址:洛谷:【P2336】[SCOI2012]喵星球上的点名 – 洛谷、BZOJ:Problem 2754. — [SCOI2012]喵星球上的点名

题目描述

a180285 幸运地被选做了地球到喵星球的留学生。他发现喵星人在上课前的点名现象非常有趣。
假设课堂上有 N 个喵星人,每个喵星人的名字由姓和名构成。喵星球上的老师会选择M 个串来点名,每次读出一个串的时候,如果这个串是一个喵星人的姓或名的子串,那么这个喵星人就必须答到。
然而,由于喵星人的字码过于古怪,以至于不能用 ASCII 码来表示。为了方便描述,a180285 决定用数串来表示喵星人的名字。
现在你能帮助 a180285 统计每次点名的时候有多少喵星人答到,以及 M 次点名结束后每个喵星人答到多少次吗?

输入输出格式

输入格式:
现在定义喵星球上的字符串给定方法:
先给出一个正整数 L ,表示字符串的长度,接下来L个整数表示字符串的每个字符。
输入的第一行是两个整数 N 和 M。
接下来有 N 行, 每行包含第 i 个喵星人的姓和名两个串。 姓和名都是标准的喵星球上的字符串。
接下来有 M 行,每行包含一个喵星球上的字符串,表示老师点名的串。

输出格式:
对于每个老师点名的串输出有多少个喵星人应该答到。
然后在最后一行输出每个喵星人被点到多少次。

输入输出样例

输入样例#1:

2 3
6 8 25 0 24 14 8 6 18 0 10 20 24 0
7 14 17 8 7 0 17 0 5 8 25 0 24 0
4 8 25 0 24
4 7 0 17 0
4 17 0 8 25

输出样例#1:

2
1
0
1 2

说明

对于30%的数据,保证:
1 ≤ N ,M ≤ 1000 喵星人的名字总长不超过 4000,点名串的总长不超过 2000,
对于100%的数据,保证:
1 ≤ N ≤ 20000, 1 ≤ M ≤ 50000 喵星人的名字总长和点名串的总长分别不超过100000
保证喵星人的字符串中作为字符存在的数不超过10000

题解

AC自动机暴力

考虑一个AC自动机的暴力做法。首先因为字符集大小并不确定,而且可能很大,我们要考虑map存儿子。接着因为字符集大小不确定我们没办法加虚边了只能暴力跳fail。这导致了最后洛谷#11 TLE的惨案。
我们考虑把点名串扔进去,把姓名用一个非法字符连接起来来匹配。每匹配到一个点名串在点名串的计数器上+1,然后统计这个姓名匹配到了多少点名串。
记录该点是否被访问过的数组vis开的很大,最好不要memset,用一个栈记下修改过的节点最后改回去。
我加了一堆register,全部inline,能避免STL的改手写,加了fread,然而还是敌不过洛谷的#11 QAQ!
所以说这是个AC自动机暴力啦,我也不会fail树高论,也不会SA,就这样吧,以后学了SA再来A这题。
对不起,SA真的是能为所欲为的!

fail树优化

如果我们将AC自动机中所有fail指针(未路径压缩)反向建有向边,会得到一棵以0为根方向向叶子节点的有向树。利用fail树的有关性质,我们可以把字符串匹配问题转换成树上问题,优化复杂度。本题中也可以使用该种算法。
首先我们可以类似的在姓名中间加一个不会出现的字符,这里我用的是-1,插入Trie树中,然后计算fail指针建立fail树。需要注意的是fail指针无法路径压缩,因此只能暴力跳fail来找。为了方便,后面建树的时候对节点编号+1处理了一下,避免0的特判。事实证明,这一操作为后面的处理带来了不少的麻烦,许多地方忘记+1换算了。有了fail树以后,我们要利用fail树的相关性质解决题目问题。
对于问题1,是求一个模式串出现在哪些文本串中了,这个问题可以转化为一个模式串的末尾点在fail树的子树中有多少不同文本串的节点。这是由于,文本串中出现过的模式串位置末尾的fail指针会指向这一模式串的末尾,由于反向,就变成了模式串末尾指向文本串节点了,那就是一个子树问题。这个问题可以转化成一个DFS序序列上的数颜色种数问题,这让我们想起了[SDOI2009]HH的项链一题,采用相同的树状数组+扫描线的做法即可求出。因此,我们需要标注每个节点属于哪个文本串这样一个颜色信息,并且处理出一个每个颜色有哪些点的链表,方便扫描线维护信息。如果文本串很多,我们不可能每次都枚举一下每个颜色,因此不如再开个链表记录一下每个点上有几种颜色。链表空间紧凑,适合用来做这种不确定空间且不需要随机查找的统计问题。另外还要处理出子树大小,方便确定查询区间。上述信息处理完毕以后,就是上面那个省选题的做法了。这一部分的复杂度是O(n \log n)的。
对于问题2,求一个文本串中出现过几种模式串。如果我们默认每个节点初始权值为0,而对每个模式串结尾节点的权值加1,则模式串到树根路径上的点权和就是这个点为结尾的后缀为模式串的个数,如果将文本串上每个节点到树根的路径求并求权值和,就是我们要的答案了。树链的并我们可以用树链剖分+线段树区间染色求和的方式来求,这样做则是O(n \log^2 n)的。如果想只有一个log,可以用将节点按照DFS序排序后的顺序统计信息,加入第i个节点的时候,减去i与i-1节点的LCA到树根的路径权值和,这样可以避免这一段公共路径权值和被计入两次,起到了去重的效果。求LCA采用倍增法,则复杂度可以降为O(n \log n)了。
但是我们仔细想一下,这些内容都很长。我们需要Trie树的插入和统计信息,需要两个链表分别维护一个节点上的颜色种类和一个颜色涉及的节点数,需要一个树状数组,需要一个树DFS统计树信息,需要一个倍增LCA的处理,需要一个离线处理区间颜色数的算法,需要一个按DFS序排序去重的计数算法。以上就是写出这道题的这种解法的所有需求,接下来就是一步一步地去实现它们了。
在实现过程中,遇到了种种问题,最困惑我的一点就是写代码的时候DFS序号和节点序号搞混了,以及节点序号+1忘记了,不过如果能够自己实现一遍,相信是很锻炼代码能力和调试能力的。这个题今天花了一天来做,收获也很大。

代码

AC自动化暴力

// Code by KSkun, 2018/3
#include <cstdio>
#include <cstring>

#include <map>
#include <vector>

struct io {
    char buf[1 << 26], *s;

    io() {
        fread(s = buf, 1, 1 << 26, stdin);
    }

    inline int read() {
        register int res = 0;
        while(*s < '0' || *s > '9') s++;
        while(*s >= '0' && *s <= '9') res = res * 10 + *s++ - '0';
        return res;
    }
} ip;

#define readint ip.read

const int MAXN = 100005, MAXL = 100005;

int n, m;
std::vector<int> names[MAXN], val[MAXL];
int fail[MAXL], cnt;
std::map<int, int> ch[MAXL];
int que[MAXL], l, r;

inline void insert(int id) {
    register int len = readint(), p = 0;
    for(register int i = 0; i < len; i++) {
        int t = readint();
        if(!ch[p][t]) ch[p][t] = ++cnt;
        p = ch[p][t];
    }
    val[p].push_back(id);
}

inline void calfail() {
    l = r = 0;
    fail[0] = 0;
    for(std::map<int, int>::iterator it = ch[0].begin(); it != ch[0].end(); it++) {
        register int u = it->second;
        que[r++] = u;
        fail[u] = 0;
    }
    while(l != r) {
        register int u = que[l++];
        for(std::map<int, int>::iterator it = ch[u].begin(); it != ch[u].end(); it++) {
            int t = it->first, v = it->second, p = fail[u];
            while(p && !ch[p][t]) p = fail[p];
            fail[v] = ch[p][t];
            que[r++] = v;
        }
    }
}

int ctn[MAXN], qur[MAXN];
bool vis[MAXL];
int sta[MAXL], stot;

inline int calcnt(int p) {
    register int res = 0;
    for(; p; p = fail[p]) {
        if(vis[p]) break;
        vis[p] = true;
        sta[++stot] = p;
        for(register int i = 0; i < val[p].size(); i++) {
            res++;
            qur[val[p][i]]++;
        }
    }
    return res;
}

inline void query(int id) {
    register int p = 0;
    stot = 0;
    for(register int i = 0; i < names[id].size(); i++) {
        int t = names[id][i];
        while(p && !ch[p][t]) p = fail[p];
        p = ch[p][t];
        ctn[id] += calcnt(p);
    }
    for(register int i = 1; i <= stot; i++) {
        vis[sta[i]] = false;
    }
}

int main() {
    n = readint(); m = readint();
    for(register int i = 1; i <= n; i++) {
        register int len = readint();
        while(len--) names[i].push_back(readint());
        names[i].push_back(-1);
        len = readint();
        while(len--) names[i].push_back(readint());
    }
    for(register int i = 1; i <= m; i++) {
        insert(i);
    }
    calfail();
    for(register int i = 1; i <= n; i++) {
        query(i);
    }
    for(register int i = 1; i <= m; i++) {
        printf("%d\n", qur[i]);
    }
    for(register int i = 1; i <= n; i++) {
        printf("%d ", ctn[i]);
    }
    return 0;
}

fail树优化(附注释)

// Code by KSkun, 2018/5
#include <cstdio>
#include <cctype>

#include <algorithm>
#include <map>
#include <queue>
#include <vector>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) 
        ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    register char c = fgc();
    while(!isdigit(c)) {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(isdigit(c)) {
        res = (res << 1) + (res << 3) + c - '0';
        c = fgc();
    }
    return res * neg;
}

const int MAXN = 200005;

// 双向链表 
struct LinkedList {
    int head[MAXN], tail[MAXN], val[MAXN], pre[MAXN], nxt[MAXN], tot;
    LinkedList() {
        tot = 0;
    }
    // 在p所属的链表头部插入v这个值 
    inline void insert(int p, int v) {
        tot++; 
        if(!tail[p]) tail[p] = tot;
        val[tot] = v; if(head[p]) pre[head[p]] = tot; nxt[tot] = head[p]; head[p] = tot;
    }
} col, col2;

// 树状数组 
struct BinaryIndexTree {
    int tr[MAXN];
    inline int lowbit(int x) {
        return x & -x;
    }
    inline void add(int x, int v) {
        for(int i = x; i < MAXN; i += lowbit(i)) {
            tr[i] += v;
        }
    }
    inline int query(int x) {
        int res = 0;
        for(int i = x; i; i -= lowbit(i)) {
            res += tr[i];
        }
        return res;
    }
} bit;

// 字符集太大,用map存儿子省空间 
std::map<int, int> ch[MAXN];
std::queue<int> que;
std::vector<int> gra[MAXN]; // fail树存在这里 
int fail[MAXN], ecnt[MAXN], end[MAXN], tot;

// Trie树的插入,分别处理文本串和模式串 
inline void insert(std::vector<int> &str, int id, bool isname) {
    int p = 0;
    for(int i = 0; i < str.size(); i++) {
        int t = str[i];
        if(!ch[p][t]) {
            ch[p][t] = ++tot;
        }
        if(isname) {
            // 统计每个节点上的颜色种类 
            col.insert(ch[p][t] + 1, id);
        }
        p = ch[p][t];
    }
    if(!isname) {
        // 统计一个节点是多少模式串的结尾,模式串可能有多个相同不可用bool 
        ecnt[p + 1]++; end[id] = p + 1;
    }
}

// 计算fail数组 
inline void calfail() {
    for(std::map<int, int>::iterator it = ch[0].begin(); it != ch[0].end(); it++) {
        que.push(it->second);
        gra[1].push_back(it->second + 1);
    }
    while(!que.empty()) {
        int p = que.front(); que.pop();
        for(std::map<int, int>::iterator it = ch[p].begin(); it != ch[p].end(); it++) {
            int t = fail[p], q = it->second;
            // 暴力跳fail计算 
            while(!ch[t][it->first] && t) t = fail[t];
            fail[q] = ch[t][it->first];
            gra[ch[t][it->first] + 1].push_back(q + 1);
            que.push(q);
        }
    }
}

int anc[MAXN][20], siz[MAXN], dfn[MAXN], dep[MAXN], ptn[MAXN], sum[MAXN], clk;

// DFS处理fail树信息 
inline void dfs(int u) {
    siz[u] = 1; dfn[u] = ++clk; ptn[dfn[u]] = u;
    sum[u] = sum[anc[u][0]] + ecnt[u];
    for(int i = 0; i < gra[u].size(); i++) {
        int v = gra[u][i];
        if(v == anc[u][0]) continue;
        dep[v] = dep[u] + 1;
        anc[v][0] = u;
        dfs(v);
        siz[u] += siz[v];
    }
}

int n, m;

struct Query {
    int id, l, r;
    inline bool operator<(const Query &rhs) const {
        return l == rhs.l ? r < rhs.r : l < rhs.l;
    }
} qur[MAXN];

// 第一个问题

int ans1[MAXN];

inline void main1() {
    // 按DFS序处理出每个颜色的链表 
    for(int i = 1; i < MAXN; i++) {
        for(int j = col.head[ptn[i]]; j; j = col.nxt[j]) {
            int c = col.val[j];
            col2.insert(c, i);
        }
    }
    // 初始化询问以及对询问排序 
    for(int i = 1; i <= m; i++) {
        qur[i] = Query {i, dfn[end[i]], dfn[end[i]] + siz[end[i]] - 1};
    }
    // 初始化树状数组,在每个颜色的初位置+1 
    for(int i = 1; i <= n; i++) {
        bit.add(col2.val[col2.tail[i]], 1);
    }
    std::sort(qur + 1, qur + m + 1);
    int lst = 1;
    // 同[SDOI2009]HH的项链
    for(int i = 1; i <= m; i++) {
        if(lst != qur[i].l) {
            for(int j = lst; j < qur[i].l; j++) {
                for(int k = col.head[ptn[j]]; k; k = col.nxt[k]) {
                    int c = col.val[k];
                    bit.add(j, -1);
                    col2.tail[c] = col2.pre[col2.tail[c]];
                    if(col2.val[col2.tail[c]]) bit.add(col2.val[col2.tail[c]], 1);
                }
            }
            lst = qur[i].l;
        }
        ans1[qur[i].id] = bit.query(qur[i].r);
    }
    for(int i = 1; i <= m; i++) {
        printf("%d\n", ans1[i]);
    }
}

// 倍增LCA
inline void preparelca() {
    for(int j = 1; j < 20; j++) {
        for(int i = 1; i < MAXN; i++) {
            anc[i][j] = anc[anc[i][j - 1]][j - 1];
        }
    }
}

inline int querylca(int x, int y) {
    if(dep[x] > dep[y]) std::swap(x, y);
    int del = dep[y] - dep[x];
    for(int i = 19; i >= 0; i--) {
        if(del & (1 << i)) {
            y = anc[y][i];
        }
    }
    if(x == y) return x;
    for(int i = 19; i >= 0; i--) {
        if(anc[x][i] != anc[y][i]) {
            x = anc[x][i]; y = anc[y][i];
        }
    }
    return anc[x][0];
}

// 第二个问题 

std::vector<int> name[MAXN];
int ans2[MAXN];

inline bool cmp(int a, int b) {
    return dfn[a] < dfn[b];
}

inline void main2() {
    for(int i = 1; i <= n; i++) {
        int p = 0;
        // 将字符串换成该字符串的节点 
        for(int j = 0; j < name[i].size(); j++) {
            p = ch[p][name[i][j]];
            name[i][j] = p + 1;
        }
        // 按DFS序计算树链的并的权值和 
        std::sort(name[i].begin(), name[i].end(), cmp);
        for(int j = 0; j < name[i].size(); j++) {
            ans2[i] += sum[name[i][j]];
            if(j) ans2[i] -= sum[querylca(name[i][j], name[i][j - 1])];
        }
    }
    for(int i = 1; i <= n; i++) {
        printf("%d ", ans2[i]);
    }
}

int len;
std::vector<int> s;

int main() {
    n = readint(); m = readint();
    for(int i = 1; i <= n; i++) {
        s.clear();
        len = readint();
        while(len--) {
            s.push_back(readint());
        }
        s.push_back(-1);
        len = readint();
        while(len--) {
            s.push_back(readint());
        }
        insert(s, i, true);
        name[i] = s;
    }
    for(int i = 1; i <= m; i++) {
        s.clear();
        len = readint();
        while(len--) {
            s.push_back(readint());
        }
        insert(s, i, false);
    }
    calfail();
    dfs(1);
    preparelca();
    main1();
    main2();
    return 0;
}
AC自动机原理及实现

AC自动机原理及实现

概述 AC自动机算法是一种常见的多串匹配算法。理解本算法需要先理解当模式串只有一个的时候的 

[SDOI2008]Sandy的卡片 题解

[SDOI2008]Sandy的卡片 题解

题目地址:洛谷:【P2463】[SDOI2008]Sandy的卡片 – 洛谷、 

[NOI2014]动物园 题解

[NOI2014]动物园 题解

题目地址:洛谷:【P2375】[NOI2014]动物园 – 洛谷、BZOJ:Problem 3670. — [Noi2014]动物园

题目描述

近日,园长发现动物园中好吃懒做的动物越来越多了。例如企鹅,只会卖萌向游客要吃的。为了整治动物园的不良风气,让动物们凭自己的真才实学向游客要吃的,园长决定开设算法班,让动物们学习算法。
某天,园长给动物们讲解KMP算法。
园长:“对于一个字符串S,它的长度为L。我们可以在O(L)的时间内,求出一个名为next的数组。有谁预习了next数组的含义吗?”
熊猫:“对于字符串S的前i个字符构成的子串,既是它的后缀又是它的前缀的字符串中(它本身除外),最长的长度记作next[i]。”
园长:“非常好!那你能举个例子吗?”
熊猫:“例S为abcababc,则next[5]=2。因为S的前5个字符为abcab,ab既是它的后缀又是它的前缀,并且找不到一个更长的字符串满足这个性质。同理,还可得出next[1] = next[2] = next[3] = 0,next[4] = next[6] = 1,next[7] = 2,next[8] = 3。”
园长表扬了认真预习的熊猫同学。随后,他详细讲解了如何在O(L)的时间内求出next数组。
下课前,园长提出了一个问题:“KMP算法只能求出next数组。我现在希望求出一个更强大num数组一一对于字符串S的前i个字符构成的子串,既是它的后缀同时又是它的前缀,并且该后缀与该前缀不重叠,将这种字符串的数量记作num[i]。例如S为aaaaa,则num[4] = 2。这是因为S的前4个字符为aaaa,其中a和aa都满足性质‘既是后缀又是前缀’,同时保证这个后缀与这个前缀不重叠。而aaa虽然满足性质‘既是后缀又是前缀’,但遗憾的是这个后缀与这个前缀重叠了,所以不能计算在内。同理,num[1] = 0,num[2] = num[3] = 1,num[5] = 2。”
最后,园长给出了奖励条件,第一个做对的同学奖励巧克力一盒。听了这句话,睡了一节课的企鹅立刻就醒过来了!但企鹅并不会做这道题,于是向参观动物园的你寻求帮助。你能否帮助企鹅写一个程序求出num数组呢?
特别地,为了避免大量的输出,你不需要输出num[i]分别是多少,你只需要输出所有num[i]+1的乘积,对1,000,000,007取模的结果即可。

输入输出格式

输入格式:
第1行仅包含一个正整数n ,表示测试数据的组数。随后n行,每行描述一组测试数据。每组测试数据仅含有一个字符串S,S的定义详见题目描述。数据保证S 中仅含小写字母。输入文件中不会包含多余的空行,行末不会存在多余的空格。

输出格式:
包含 n 行,每行描述一组测试数据的答案,答案的顺序应与输入数据的顺序保持一致。对于每组测试数据,仅需要输出一个整数,表示这组测试数据的答案对 1,000,000,007 取模的结果。输出文件中不应包含多余的空行。

输入输出样例

输入样例#1:

3
aaaaa
ab
abcababc

输出样例#1:

36
1
32 

说明

测试点编号 约定
1 N ≤ 5, L ≤ 50
2 N ≤ 5, L ≤ 200
3 N ≤ 5, L ≤ 200
4 N ≤ 5, L ≤ 10,000
5 N ≤ 5, L ≤ 10,000
6 N ≤ 5, L ≤ 100,000
7 N ≤ 5, L ≤ 200,000
8 N ≤ 5, L ≤ 500,000
9 N ≤ 5, L ≤ 1,000,000
10 N ≤ 5, L ≤ 1,000,000

题解

参考资料:[省选前题目整理][BZOJ 3670][NOI 2014]动物园(KMP) – CSDN博客【题解】NOI2014动物园 – Twilight_Sx – 博客园
我们回想一下KMP不加优化的时候的fail数组的意义。它指的是当前位置之前的子串中最长的与某一前缀相同的后缀长度。我们利用这个来找num数组。
num数组指的是不重叠的与某一前缀相同的后缀数量,我们退一步,先不求不重叠,用一个数组cnt表示与某一前缀相同的后缀数量。我们可以把fail数组计算出来以后,利用fail算出cnt的值。
cnt[i] = cnt[fail[i]] + 1
其中cnt[1] = 1。
举个例子:ababdefghabab,abab是整个串的最长相同前后缀,cnt[3]=2,由于abab内部相同的部分a和ab在后面的abab中出现过,因此cnt[13]至少有cnt[3]中这么多,而abab本身也构成了相同前后缀,因此cnt[13]要比cnt[3]多1。
现在我们拿到了这个串的cnt数组,要怎么求num呢?我们考虑cnt数组是通过fail递推而来的,如果说cnt[i]中计入了某个重叠了的前后缀,那么某个j=若干层嵌套fail[i]的cnt[j]就是不重复的答案。这时候回到fail的前后缀长度的意义上,只要这个前后缀的长度的2倍不超过当前串长,答案是不是就不会重复了,那么j需要满足的条件便是2j≤i。
因此这个算法的复杂度是O(n)的。
注意本题的KMP写法与KMP算法原理与实现 | KSkun’s Blog中的并不相同,因为该文章中的算法字符串下标从0开始标,利用-1判断是否到头,这会给计算cnt造成麻烦,因此采用了其他同学的写法。

代码

// Code by KSkun, 2018/3
#include <cstdio>
#include <cstring>

typedef long long LL;

const int MAXN = 1000005, MO = 1e9 + 7;

int n, fail[MAXN], num[MAXN];
LL ans;
char str[MAXN];

inline void calfail() {
    memset(fail, 0, sizeof(fail));
    memset(num, 0, sizeof(num));
    int i = 2, j = 0;
    num[1] = 1;
    for(; str[i]; i++) {
        while(j && str[j + 1] != str[i]) {
            j = fail[j];
        }
        if(str[j + 1] == str[i]) j++;
        fail[i] = j;
        num[i] = num[j] + 1;
    }
}

inline void match() {
    calfail();
    ans = 1;
    int i = 2, j = 0;
    for(; str[i]; i++) {
        while(j && str[j + 1] != str[i]) {
            j = fail[j];
        }
        if(str[j + 1] == str[i]) j++;
        while(j << 1 > i) j = fail[j];
        ans = ((num[j] + 1) * ans) % MO;
    }
}

int main() {
    scanf("%d", &n);
    while(n--) {
        scanf("%s", str + 1);
        match();
        printf("%lld\n", ans);
    }
    return 0;
}
KMP算法原理与实现

KMP算法原理与实现

概述 KMP是一种字符串匹配算法,其复杂度已经达到了该类算法的下界,即,其中T是文本串,P 

[SDOI2010]外星千足虫 题解

[SDOI2010]外星千足虫 题解

题目地址:洛谷:【P2447】[SDOI2010]外星千足虫 – 洛谷、BZO 

高斯-约旦消元法原理及实现

高斯-约旦消元法原理及实现

概述

高斯-约旦消元是一种方便计算线性方程组的方法。下面我们介绍它的原理及实现。

原理

思想

各位还记得我们是怎么解多元一次方程的吗?核心思想就是消元二字。我们通过对方程进行等价变换来消元,再把得到的解回代解出所有未知数的解。这就是我们接下来要做的事情。

高斯消元

我们尝试把一个线性方程组写成矩阵的形式,比如
\begin{cases} 2x +y -z =8 \\ -3x -y +2z =-11 \\ -2x +y +2z =-3 \end{cases}
可以写成
\left[ \begin{array}{rrr|r} 2&1&-1&8 \\ -3&-1&2&-11 \\ -2&1&2&-3 \end{array} \right]
对于这个矩阵,我们尝试把左边的部分转化为一个行梯阵势。首先对于当前正在处理的行i,我们需要找到一个r>i且绝对值最大的a_{r, i},交换r行和当前行,再继续做处理,这样可以提高数值稳定性。
我们对于正在处理的行i,使下面的每一行j的a_{j, i}变为0,即整个方程j减去a_{j, i} / a_{i, i}倍的方程i。这样,对于每一行,下面的行的第一个非0项就在第i项之后了。如果这个方程组有解,最后一个方程可以解出最后一个未知数的值,再回代求所有未知数即可。

高斯-约旦消元

上面讲述了高斯消元的过程,那么高斯-约旦消元与高斯消元有何不同呢?就是下面要介绍的内容了。
如果我们能通过类似的等价变换得到一个除主对角线以外的地方全0的矩阵,那就不用回代解未知数了,岂不美哉?。我们考虑如何用变换实现这一目标。
如果我们处理的地方不只是i行下面的每一行,而是除了i行以外的每一行。我们尝试把i列上除了a_{i, i}以外的元素全部变为0。因为在计算到该行的时候,i以前的系数全部变为0了,因此只需要对所有行的i以后的系数做变化。这样计算一番后,矩阵中留下的系数就只有a_{i, i}了,也就是说,方程组变成了a_{i, i}x_i = a_{i, n + 1},做一下除法就求出了未知数。
这种方法的计算比高斯消元要多一些,但是省去了回代的过程,代码也短些。但是不可避免的就是复杂度仍然是O(n^3)

代码

本代码可以通过【P3389】【模板】高斯消元法 – 洛谷一题。

// Code by KSkun, 2018/3
#include <cstdio>
#include <cmath>

#include <algorithm>

const int MAXN = 105;
const double EPS = 1e-10;

int n;
double mat[MAXN][MAXN];
bool fail;

inline void gauss_jordan() {
    for(int i = 1; i <= n; i++) {
        int r = i;
        for(int j = i + 1; j <= n; j++) {
            if(std::fabs(mat[r][i]) < std::fabs(mat[j][i])) {
                r = j;
            }
        }
        if(r != i) {
            for(int j = 1; j <= n + 1; j++) {
                std::swap(mat[i][j], mat[r][j]);
            }
        }
        if(fabs(mat[i][i]) < EPS) {
            fail = true;
            return;
        }
        for(int j = 1; j <= n; j++) {
            if(j != i) {
                double t = mat[j][i] / mat[i][i];
                for(int k = i + 1; k <= n + 1; k++) {
                    mat[j][k] -= mat[i][k] * t;
                }
            }
        }
    }
    for(int i = 1; i <= n; i++) {
        mat[i][n + 1] /= mat[i][i];
    }
} 

int main() {
    scanf("%d", &n);
    for(int i = 1; i <= n; i++) {
        for(int j = 1; j <= n + 1; j++) {
            scanf("%lf", &mat[i][j]);
        }
    }
    gauss_jordan();
    if(fail) {
        puts("No Solution");
    } else {
        for(int i = 1; i <= n; i++) {
            printf("%.2lf\n", mat[i][n + 1]);
        }
    }
    return 0;
}

参考资料

数学笔记:逆元

数学笔记:逆元

逆元 一整数a对同余n之模逆元是满足以下公式的整数b 也可以写成 求逆元的方法 扩展欧几里