NOIp 2018 前的数学板子总结

发表于 更新于 分类于
NOIp范围内的数学定理及算法

数论

Euclidean algorithm

用于求两个数的最大公因数, 也称辗转相除法。

证明:

设\(z \mid x\), \(z \mid y\), 则\(z \mid (y - x)\)。

设\(z\)不是\(x\)的因子, 则\(z\)不是\(x\), \(y - x\)的公因子。

设\(z \mid x\), \(z\)不是\(y\)的因子, 则\(z\)不是\(x\), \(y - x\)的公因子。

代码: 

1
2
3
4
template<typename IntegerType> 
inline IntegerType Euclidean(const IntegerType &a, const IntegerType &b) {
  return b ? Euclidean(b, a % b) : a;
}
或 
1
2
3
4
5
template<typename IntegerType>
inline IntegerType Euclidean(register IntegerType a, register IntegerType b) {
  if (b) while (b ^= a ^= b ^= a %= b) {}
  return a;
}

Extended Euclidean algorithm

用于求在已知\((a, b)\)时, 求解一组\((x, y)\), 使得\(ax+by=(a, b)\)。

首先, 书上说根据数论中的相关定理, 解一定存在。

其次, 因为\((a, b) = (b, a \bmod b)\), 所以

\[ \begin {aligned} ax + by &= (a, b) \\\\ &= (b, a \bmod b) \\\\ &= bx + (a \bmod b)y \\\\ &= bx + \left(a - \left\lfloor \frac{a}{b} \right\rfloor b\right)y \\\\ &= ay + \left(x - \left\lfloor \frac{a}{b} \right\rfloor y\right)b \end {aligned} \]

根据前面的结论: \(a\)和\(b\)都在减小, 当\(b\)减小到\(0\)时, 就可以得出\(x = 1\), \(y = 0\)。然后递归回去就可以求出最终的\(x\)和\(y\)了。

代码:

1
2
3
4
5
6
7
8
9
template<typename IntegerType>
inline void ExtendedEuclidean(const IntegerType &a, const IntegerType &b, register IntegerType &gcd, register IntegerType &x, register IntegerType &y) {
  if (!b) {
    x = 1, y = 0, gcd = a;
  } else {
    ExtendedEuclidean(b, a % b, gcd, y, x), 
    y -= a / b * x;
  }
}

Modular multiplicative inverse

若有\(ax \equiv 1 \pmod b\)(其中\(a\), \(b\)互素), 则称\(x\)为\(a\)的逆元, 记为\(a^{-1}\)。

因此逆元有如下性质:

\[a \cdot a^{-1} \equiv 1 \pmod b\]

逆元的一大应用是模意义下的除法, 除法在模意义下并不是封闭的，但我们可以根据上述公式，将其转化为乘法。

\[\frac{x}{a} = x \cdot a^{-1} \pmod b\]

Quick Power

根据Fermat's little theorem(即\(a^p \equiv a \pmod p\), 其中\(p\)为素数且\(a\)为\(p\)的倍数。若\(a\), \(p\)互素, 则有\(a^{p - 1} \equiv 1 \pmod p\))的公式, 变形可以得到

\[a \cdot a^{p - 2} \equiv 1 \pmod p\]

根据乘法逆元的定义, \(a^{p - 2}\)即为\(a\)的乘法逆元。

使用快速幂计算, 时间复杂度为\(\Theta(\lg p)\)

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
template<typename IntegerType>
inline IntegerType QuickPower(register IntegerType base, register IntegerType times, const IntegerType &kMod) {
  register IntegerType ret(1);
  while (times) {
    if (times & 1) ret = ret * base % kMod;
    base = base * base % kMod,
    times >>= 1;
  }
  return ret % kMod;
}
template<typename IntegerType>
inline IntegerType ModularMultiplicativeInverse(const IntegerType a, const IntegerType &p) {
  return QuickPower(a, p - 2, p);
}

Extended Euclidean algorithm

Extended Euclidean algorithm用来求解方程\(ax + by = (a, b)\)的一组解, 其中, 当\(b\)为素数时, 有\((a, b) = 1\), 则有

\[ax \equiv 1 \pmod b\]

时间复杂度: \(\Theta(\lg a)\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
template<typename IntegerType>
inline void ExtendedEuclidean(const IntegerType &a, const IntegerType &b, register IntegerType &x, register IntegerType &y) {
  if (!b) {
    x = 1, y = 0;
  } else {
    ExtendedEuclidean(b, a % b, y, x), 
    y -= a / b * x;
  }
}
template<typename IntegerType>
inline IntegerType ModularMultiplicativeInverse(const IntegerType a, const IntegerType &p) {
  register IntegerType x, y;
  ExtendedEuclidean(a, p, x, y);
  return (x % p + p) % p;
}

Recurse

设\(t = \left\lfloor \frac{p}{i} \right\rfloor\), \(k = p \bmod i\), 那么

\[ti + k \equiv 0 \pmod p \Rightarrow -ti \equiv k \pmod p\]

两边同时除以\(ki\), 得到

\[-tk^{-1} \equiv i^{-1} \pmod p\]

把\(t\)和\(k\)替换回来, 得到递推式

\[i^{-1} = \left(p - \left\lfloor \frac{p}{i} \right\rfloor\right)(p \bmod i)^{-1} \bmod p\]

其中\(i \leq p\)。

时间复杂度: \(\Theta(a)\)

代码: 

1
2
3
4
5
6
7
8
template<typename IntegerType>
inline IntegerType ModularMultiplicativeInverse(const IntegerType a, const IntegerType &p) {
  register IntegerType inverse[100000] = {0, 1};
  for (register IntegerType i(2); i <= a; ++i) {
    inverse[i] = ((-p / i * inverse[p % i]) % p + p) % p;
  }
  return inverse[a];
}

Chinese remainder theorem

设自然数\(m_1, m_2, \cdots , m_r\)两两互素, 并记\(N = m_1m_2\cdots m_r\), 则同余方程组:

\[ \begin{cases} x \equiv a_1 \pmod {m_1} \\\\ x \equiv a_2 \pmod {m_2} \\\\ \cdots \\\\ x \equiv a_r \pmod {m_r} \end{cases} \]

在模\(N\)同余的意义下有唯一解。

解法:

考虑方程组\((1 \leq i \leq r)\):

\[ \begin{cases} x \equiv 0 \pmod {m_1} \\\\ \cdots \\\\ x \equiv 0 \pmod {m_{i - 1}} \\\\ x \equiv 1 \pmod {m_i} \\\\ x \equiv 0 \pmod {m_{i + 1}} \\\\ \cdots \\\\ x \equiv 0 \pmod {m_r} \end{cases} \]

由于各\(m_i\)(\(1 \leq i \leq r\))两两互素, 这个方程组作变量替换, 令\(x = \left\lfloor \frac{N}{m_i} \right\rfloor y\), 方程组等价于解同余方程: \(\left\lfloor \frac{N}{m_i} \right\rfloor y \equiv 1 \pmod {m_i}\), 若要得到特解\(y_i\), 只要令:\(x_i = \left\lfloor \frac{N}{m_i} \right\rfloor y_i\), 则方程组的解为:\(x_0 = a_1x_1 + a_2x_2 + \cdots + a_rx_r \pmod N\), 在模\(N\)的意义下唯一。

时间复杂度: \(\Theta(n \lg a)\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
template<typename IntegerType>
inline void ExtendedEuclidean(const IntegerType &a, const IntegerType &b, register IntegerType &x, register IntegerType &y) {
  if (!b) {
    x = 1, y = 0;
  } else {
    ExtendedEuclidean(b, a % b, y, x), 
    y -= a / b * x;
  }
}
template<typename IntegerType>
inline IntegerType ChineseRemainderTheorem(const std::vector<IntegerType> &a, const std::vector<IntegerType> &m) {
  register IntegerType N(1), x, y, ret(0);
  for (auto i : m) {
    N *= i;
  }
  for (register int i(0), maximum(a.size()); i < maximum; ++i) {
    register IntegerType b(N / m[i]);
    ExtendedEuclidean(b, m[i], x, y),
    ret = (ret + b * x * a[i]) % N;
  }
  return (ret % N + N) % N;
}

Extended Chinese remainder theorem

并不是所有的同余方程组的\(m_i\)(\(1 \leq i \leq r\))都互素, 这时候就需要用到扩展中国剩余定理。

对于同余方程组:

\[ \begin{cases} x \equiv a_1 \pmod {m_1} \\\\ x \equiv a_2 \pmod {m_2} \\\\ \cdots \\\\ x \equiv a_r \pmod {m_r} \end{cases} \]

我们首先只考虑其中的两个方程, 可以得到

\[ \begin{cases} x = a_1 + k_1m_1 \\\\ x = a_2 + k_2m_2 \end{cases} \Rightarrow a_1 + k_1m_1 = a_2 + k_2m_2 \\\\ \Rightarrow k_2m_2 - k_1m_1 = a_1 - a_2 \]

其形式类似于\(ax + by = c\)

则当\((m_1, m_2) \nmid (a_1 - a_2)\)时, 方程无解。

若\((m_1, m_2) \mid (a_1 - a_2)\), 就用Extended Euclidean algorithm求出\(m_1x + m_2y = (m_1, m_2)\)的\(y\), 两边同时乘上\(\frac{a_1 - a_2}{(m_1, m_2)}\), 就得到了\(k_1\)。

然后反推\(x\), 得到\(x = a_1 - k_1m_1\)。

这个\(x\)适用于这两个方程, 设它为\(x_0\), 就得到了通解: \(x = x_0 + k[m_1, m_2]\), 且原两个方程与该方程是等价的。

我们把这个方程稍微转化一下, 就得到了新的同余方程: \(x \equiv x_0 \pmod {[m_1, m_2]}\), 以此类推, 得到的最后的方程的最小非负数解就是我们要找的答案。

时间复杂度: \(\Theta(n \lg b)\)

题目链接: POJ 2891 Strange Way to Express Integers

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#include <cstdio>
#include <vector>
int n;
template<typename IntegerType>
inline void ExtendedEuclidean(const IntegerType &a, const IntegerType &b, register IntegerType &gcd, register IntegerType &x, register IntegerType &y) {
  if (!b) {
    x = 1, y = 0, gcd = a;
  } else {
    ExtendedEuclidean(b, a % b, gcd, y, x), 
    y -= a / b * x;
  }
}
template<typename IntegerType>
inline IntegerType ExtendedChineseRemainderTheorem(const std::vector<IntegerType> &a, const std::vector<IntegerType> &m) {
  register IntegerType N(m[1]), ret(a[1]), x, y, gcd;
  for (register int i(0), maximum(a.size()); i < maximum; ++i) {
    ExtendedEuclidean(N, m[i], gcd, x, y);
    if ((ret - a[i]) % gcd) return -1;
    x = (ret - a[i]) / gcd * x % m[i],
    ret -= N * x;
    N = N / gcd * m[i];
    ret %= N;
  }
  return (ret % N + N) % N;
}
std::vector<long long> a, r;
int main(int argc, char const *argv[]) {
  while (~scanf("%d", &n)) {
    a.clear(), r.clear();
    for (register long long i(1), _a, _r; i <= n; ++i) {
      scanf("%lld %lld", &_a, &_r);
      a.push_back(_a), r.push_back(_r);
    }
    printf("%lld\n", ExtendedChineseRemainderTheorem(r, a));
  }
  return 0;
}

素数判定

素数判定的暴力方法

时间复杂度: \(\Theta\left(\sqrt{n}\right)\)

1
2
3
4
5
6
7
template<typename IntegerType>
inline bool IsPrime(const IntegerType &n) {
  for (register int i(2); i * i <= n; ++i) {
    if (!(n % i)) return false;
  }
  return true;
}

Sieve of Eratosthenes

又称素数的线性筛法。

方法如下:

  1. 使\(p\)等于\(2\), 即最小的素数
  2. 将表中所有\(p\)的倍数标记为合数
  3. 使\(p\)等于表中大于\(p\)的最小素数, 若没有则结束
  4. 重复第\(2\)步

时间复杂度: \(\Theta(n \lg \lg n)\)

题目链接: Luogu P3383 【模板】线性筛素数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
#include <cstdio>
#include <vector>
int n, m;
bool isnt_prime[10000010] = {1, 1};
std::vector<int> prime;
inline void SieveOfEratosthenes() {
  for (register int i(2); i <= n; ++i) {
    if (!isnt_prime[i]) {
      prime.push_back(i);
    }
    for (auto j : prime) {
      if (i * j > n) break;
      isnt_prime[i * j] = 1;
      if (!(i % j)) break;
    }
  }
}
int main(int argc, char const *argv[]) {
  scanf("%d %d", &n, &m);
  SieveOfEratosthenes();
  for (register int i(0), x; i < m; ++i) {
    scanf("%d", &x),
    puts(isnt_prime[x] ? "No" : "Yes");
  }
  return 0;
}

Miller-Rabin primality test

这是一种随机性素性测试算法, 结果可能出错, 但可能性极小。

费马小定理:

若\(p\)为素数, \(a\)为正整数, 且\((a, p) \ne 1\), 则\(a^p \equiv a \pmod p\)。

若有\((a, p) = 1\), 则\(a^{p - 1} \equiv 1 \pmod p\)。

但对于一些数\(p\), 它们满足\(a^p \equiv a \pmod p\), 但却是合数, 我们将它们定义为伪素数(Pseudoprime)。

费马小定理的逆定理并不成立, 但很多时候是可行的。

Miller-Rabin正是基于费马小定理:

取多个底\(a\)使得\((a, n) = 1\), 测试是否有\(a^{p - 1} \equiv 1 \pmod p\), 若成立, 则可以近似地将\(n\)看作素数。

还有一类被称为Carmichael数的数, 它们满足\(a^{p - 1} \equiv 1 \pmod p\), 但是是合数。

为了减少Carmichael数对素性测试的影响, 我们引进一个引理:

\(1\)和\(-1\)的平方模\(p\)总得到\(1\), 称它们为\(1\)的平凡平方根; 同样地, 有\(1\)的非平凡平方根。使得\(x\)为一个与\(1\)关于模\(p\)同余的数的平方根, 那么:

\[ x^{2}\equiv 1{\pmod {p}} \\\\ (x-1)(x+1)\equiv 0{\pmod {p}} \]

换句话说, \(p \mid (x - 1)(x + 1)\)。根据Euclid's lemma, \(p \mid (x - 1)\)或\(p \mid (x + 1)\), 也就是\(x \equiv 1 \pmod p\)或\(x \equiv -1 \pmod p\)。

由此推出, 若\(p\)为素数, 那么\(x^2 \equiv 1 \pmod p\)(\(0 < x < p\))的解为\(\begin{cases}x_1 = 1 \\\\ x_2 = p - 1\end{cases}\)

时间复杂度: \(\Theta(s \log_3n)\)(\(s\)为测试次数)

题目链接: hihoCoder #1287 : 数论一·Miller-Rabin质数测试

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#include <ctime>
#include <cstdio>
#include <cstdlib>
template<typename IntegerType>
inline IntegerType QuickMultiplication(register IntegerType base, register IntegerType times, const IntegerType &kMod) {
	register IntegerType ret(0);
	base %= kMod;
	while (times) {
		if (times & 1) ret = (ret + base) % kMod;
		base = (base << 1) % kMod,
		times >>= 1;
	}
	return ret % kMod;
}
template<typename IntegerType>
inline IntegerType QuickPower(register IntegerType base, register IntegerType times, const IntegerType &kMod) {
  register IntegerType ret(1);
	base %= kMod;
	while (times) {
		if (times & 1) ret = QuickMultiplication(ret, base, kMod);
		base = QuickMultiplication(base, base, kMod),
		times >>= 1;
	}
	return ret % kMod;
}
template<typename IntegerType>
inline bool MillerRabin(const IntegerType &n) {
	if (n <= 2) return n == 2;
	if (!(n & 1)) return false;
	register int times(0);
	register IntegerType a, x, y, u(n - 1);
	while (!(u & 1)) ++times, u >>= 1;
	for (register int i(0); i < 10; ++i) {
		a = rand() % (n - 1) + 1,
		x = QuickPower(a, u, n);
		for (register int j(0); j < times; ++j) {
			y = QuickMultiplication(x, x, n);
			if (y == 1 && x != 1 && x != n - 1) return false;
			x = y;
		}
		if (x != 1) return false;
	}
	return true;
}
int T;
long long num;
int main(int argc, char const *argv[]) {
	srand(time(nullptr));
	scanf("%d", &T);
	while (T--) {
		scanf("%lld", &num);
		puts(MillerRabin(num) ? "Yes" : "No");
	}
	return 0;
}

Euler's totient function的线性筛选法

下面几个性质我就不证了知道有用就行！

其实是因为作为一个背板的蒟蒻OIer根本不会证

  1. \(\phi(p) = p - 1\)
  2. 如果\(i \bmod p = 0\), 那么\(\phi(ip) = p\phi(i)\)
  3. 若\(i \bmod p \ne 0\), 那么\(\phi(ip) = (p - 1)\phi(i)\)

时间复杂度: 趋近于\(\Theta(n)\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
template<typename IntegerType>
inline void SieveOfPhi() {
  phi[1] = 1;
  for (register IntegerType i(2); i <= n; ++i) {
    if (!iscomposite[i]) {
      prime.push_back(i), 
      phi[i] = i - 1;
    }
    for (auto j : prime) {
      if (i * j > n) break;
      iscomposite[i * j] = true;
      if (!(i % j)) {
        phi[i * j] = phi[i] * j;
        break;
      } else {
        phi[i * j] = phi[i] * (j - 1);
      }
    }
  }
}

求一个数的\(\phi\)函数值

减去与它不互素的就行。

时间复杂度: \(\Theta\left(\sqrt n\right)\)

1
2
3
4
5
6
7
8
9
10
11
template<typename IntegerType>
inline IntegerType Phi(register IntegerType n) {
  register IntegerType ret(n);
  for (register IntegerType i(2); i * i <= n; ++i) {
    if (!(n % i)) {
      ret -= ret / i;
      while (!(n % i)) n /= i;
    }
  }
  return n > 1 ? ret - ret / i : ret;
}

组合数学

Lucas's theorem

Lucas's theorem是用来求\({n \choose m} \bmod p\)的值。其中: \(n\)和\(m\)时非负整数, \(p\)是素数。

Lucas's theorem的结论:

1. $Lucas(n, m, p) = cm(n \bmod p, m \bmod p) \cdot Lucas\left(\left\lfloor \frac{n}{p} \right\rfloor, \left\lfloor \frac{m}{p} \right\rfloor, p\right)$; $Lucas(x, 0, p) = 1$; 其中, $cm(a, b) = a!(b!(a - b)!)^{p - 2} \bmod p = \frac{a!(b!)^{p - 2}}{(a - b)!}$。 2. 对于非负整数$m$和$n$以及素数$p$, 将$n$与$m$以$p$进制方式表达, 即 $$m=m_{k}p^{k}+m_{k-1}p^{k-1}+\cdots +m_{1}p+m_{0}$$ $$n=n_{k}p^{k}+n_{k-1}p^{k-1}+\cdots +n_{1}p+n_{0}$$ 以下同余关系成立: $${m\choose n}\equiv \prod _{i=0}^{k}{{m_{i}}\choose {n_{i}}}{\pmod {p}}$$

时间复杂度: \(\Theta(\log_pn)\)

题目链接: Luogu P3807 【模板】卢卡斯定理

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#include<cstdio>
long long a[100010];
template<typename IntegerType>
inline IntegerType pow(register IntegerType base, register IntegerType times, const IntegerType &kMod) {
	register IntegerType ans(1);
	base %= kMod;
	while (times) {
		if (times & 1) ans = ans * base % kMod;
		base = base * base % kMod,
		times >>= 1;
	}
	return ans;
}
template<typename IntegerType>
inline IntegerType C(const IntegerType &n, const IntegerType &m, const IntegerType &kMod) {
    return m > n ? 0 : a[n] * pow(a[m], kMod - 2, kMod) % kMod * pow(a[n - m], kMod - 2, kMod) % kMod;
}
template<typename IntegerType>
inline IntegerType Lucas(const IntegerType &n, const IntegerType &m, const IntegerType &kMod) {
    return m ? C(n % kMod, m % kMod, kMod) * Lucas(n / kMod, m / kMod, kMod) % kMod : 1;
}
int T;
long long n, m, p;
int main(int argc, char const *argv[]) {
    scanf("%d", &T);
    while (T--) {
			register long long n, m;
			scanf("%lld %lld %lld", &n, &m, &p);
        a[0] = 1;
        for (register long long i(1); i <= p; ++i) {
					a[i] = a[i - 1] * i % p;
				}
				printf("%lld\n", Lucas(n + m, n, p));
    }
}

矩阵

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
#include <cmath>
#include <vector>
#include <algorithm>
#include <iostream>
#include <cassert>
#define EPS 1e-8
class matrix {
 private:
  std::vector<std::vector<double> > mat;
  unsigned long swap_times;
 public:
  matrix();//初始化为空的构造函数
  matrix(const matrix&);//初始化为另一个矩阵的构造函数
  matrix(const unsigned long&, const unsigned long&);//同初始化的构造函数
  std::vector<double>& operator [](const unsigned long &lines) {//下标运算符
    return mat[lines];
  }
  void assign(const unsigned long&, const unsigned long&);//初始化
  friend std::ostream& operator <<(std::ostream&, const matrix&);//输出流
  matrix operator +(const matrix &another) const {//加法
    assert(this->mat.size() == another.mat.size() && this->mat.size() ? this->mat[0].size() == another.mat[0].size() : true);
    if (!this->mat.size() || !this->mat[0].size()) return matrix();
    matrix ret(this->mat.size() - 1, this->mat[0].size() - 1);
    for (unsigned long i(1), row(this->mat.size()); i < row; ++i) {
      for (unsigned long j(1), column(this->mat[0].size()); j < column; ++j) {
        ret[i][j] = this->mat[i][j] + another.mat[i][j];
      }
    }
    return ret;
  }
  matrix& operator +=(const matrix &another) {
    *this = *this + another;
    return *this;
  }
  matrix operator -(const matrix &another) const {//减法
    assert(this->mat.size() == another.mat.size() && this->mat.size() ? this->mat[0].size() == another.mat[0].size() : true);
    if (!this->mat.size() || !this->mat[0].size()) return matrix();
    matrix ret(this->mat.size() - 1, this->mat[0].size() - 1);
    for (unsigned long i(1), row(this->mat.size()); i < row; ++i) {
      for (unsigned long j(1), column(this->mat[0].size()); j < column; ++j) {
        ret[i][j] = this->mat[i][j] - another.mat[i][j];
      }
    }
    return ret;
  }
  matrix& operator -=(const matrix &another) {
    *this = *this - another;
    return *this;
  }
  matrix operator *(const double &number) const {//数乘
    if (!this->mat.size() || !this->mat[0].size()) return matrix();
    matrix ret(this->mat.size() - 1, this->mat[0].size() - 1);
    for (unsigned long i(1), row(this->mat.size()); i < row; ++i) {
      for (unsigned long j(1), column(this->mat[0].size()); j < column; ++j) {
        ret[i][j] = this->mat[i][j] * number;
      }
    }
    return ret;
  }
  friend matrix operator *(const double&, const matrix&);//数乘
  matrix& operator *=(const double &number) {
    *this = *this * number;
    return *this;
  }
  matrix operator ~() const {//转置
    if (!this->mat.size() || !this->mat[0].size()) return matrix();
    matrix ret(this->mat[0].size() - 1, this->mat.size() - 1);
    for (unsigned long i(1), row(this->mat.size()); i < row; ++i) {
      for (unsigned long j(1), column(this->mat[0].size()); j < column; ++j) {
        ret[j][i] = this->mat[i][j];
      }
    }
    return ret;
  }
  matrix operator *(const matrix &another) const {//矩阵乘法
    if (!this->mat.size() || !this->mat[0].size() || !another.mat.size() || !another.mat[0].size()) return matrix();
    assert(this->mat[0].size() == another.mat.size());
    matrix ret(this->mat.size() - 1, another.mat[0].size() - 1);
    for (unsigned long i(1), row(this->mat.size()); i < row; ++i) {
      for (unsigned long j(1), column(another.mat[0].size()); j < column; ++j) {
        ret.mat[i][j] = 0.00000000;
        for (unsigned long k(1), tmp(this->mat.size()); k < tmp; ++k) {
          ret.mat[i][j] += this->mat[i][k] * another.mat[k][j];
        }
      } 
    }
    return ret;
  }
  matrix operator *=(const matrix &another) {
    *this = *this * another;
    return *this;
  }
  matrix operator ^(unsigned long times) const {//矩阵快速幂
    if (!this->mat.size() || !this->mat[0].size()) return matrix();
    assert(this->mat[0].size() == this->mat.size());
    matrix ret(this->mat.size() - 1, this->mat.size() - 1), base(*this);
    while (times) {
      if (times & 1) ret = ret * base;
      base  = base * base,
      times >>= 1;
    }
    return ret;
  }
  matrix swap_row(const unsigned long&, const unsigned long&);//交换两行
  matrix swap_column(const unsigned long&, const unsigned long&);//交换两列
  matrix eliminate();//高斯消元
  double det();//行列式
  matrix cofactor(const unsigned long&, const unsigned long&);//余子式
  matrix algebraic_cofactor(const unsigned long&, const unsigned long&);//代数余子式
  matrix principal_minor(const unsigned long&);//主子式
};
matrix::matrix() {
  mat = std::vector<std::vector<double> >();
}
matrix::matrix(const matrix &mat) {
  *this = mat;
}
matrix::matrix(const unsigned long &row, const unsigned long &column) {
  mat.clear(), mat.assign(row + 1, std::vector<double>(column + 1, 0.000000));
  for (unsigned long i(0), maximum(std::min(row, column)); i <= maximum; mat[i][i] = 1.000000, ++i);
}
void matrix::assign(const unsigned long &row, const unsigned long &column) {
  mat.clear(), mat.assign(row + 1, std::vector<double>(column + 1, 0.000000));
  for (unsigned long i(0), maximum(std::min(row, column)); i <= maximum; mat[i][i] = 1.000000, ++i);
}
std::ostream& operator <<(std::ostream &os, const matrix &mat) {
  for (unsigned long i(1), row(mat.mat.size()); i < row; ++i) {
    for (unsigned long j(1), column(mat.mat[0].size()); j < column; ++j) {
      os << mat.mat[i][j] << (j == column - 1 ? '\n' : ' ');
    }
  }
  return os;
}
matrix operator *(const double &number, const matrix &mat) {
  if (!mat.mat.size() || !mat.mat[0].size()) return matrix();
  matrix ret(mat.mat.size() - 1, mat.mat[0].size() - 1);
  for (unsigned long i(1), row(mat.mat.size()); i < row; ++i) {
    for (unsigned long j(1), column(mat.mat[0].size()); j < column; ++j) {
      ret[i][j] = mat.mat[i][j] * number;
    }
  }
  return ret;
}
matrix matrix::swap_row(const unsigned long &index1, const unsigned long &index2) {
  assert(index1 > 0 && index1 < this->mat.size()),
  assert(index2 > 0 && index2 < this->mat.size());
  matrix ret(*this);
  if (index1 == index2) return ret;
  std::swap(ret[index1], ret[index2]);
  return ret;
}
matrix matrix::swap_column(const unsigned long &index1, const unsigned long &index2) {
  matrix ret(*this);
  return ~((~ret).swap_row(index1, index2));
}
matrix matrix::eliminate() {
  swap_times = 0;
  matrix ret(*this);
  unsigned long h(1), k(1), m(ret.mat.size() - 1), n(ret.mat[0].size() - 1), i_max;
  double f;
  while (h < m && k <= n) {
    i_max = h;
    for (unsigned long i(h + 1); i <= m; ++i) {
      i_max = fabs(ret[i_max][k]) > fabs(ret[i][k]) ? i_max : i;
    }
    if (fabs(ret[i_max][k]) < EPS) ++k;
    else {
      ret.swap_row(h, i_max), swap_times += h != i_max;
      for (unsigned long i(h + 1); i <= m; ++i) {
        f = ret[i][k] / ret[h][k];
        ret[i][k] = 0.00000000;
        for (unsigned long j(k + 1); j <= n; ++j) {
          ret[i][j] -= ret[h][j] * f;
        }
      }
      ++h, ++k;
    }
  }
  return ret;
}
double matrix::det() {
  matrix tmp(this->eliminate());
  double ret(1.00000000);
  if (!tmp.mat.size()) return ret;
  assert(tmp.mat.size() == tmp.mat[0].size());
  for (int i(1), lines(tmp.mat.size()); i < lines; ++i) {
    ret *= tmp[i][i];
  }
  return ret * (swap_times & 1 ? -1.00000000 : 1.00000000);
}
matrix matrix::cofactor(const unsigned long &row, const unsigned long &column) {
  matrix ret;
  if (!this->mat.size() || !this->mat[0].size()) return ret;
  ret.mat.assign(this->mat.size() - 1, std::vector<double>(1, 0.00000000));
  unsigned long h(1), m(this->mat.size()), n(this->mat[0].size() - 1);
  for (unsigned long i(1); i < m; ++i) {
    if (i == row) continue;
    if (column == 1) {
      ret.mat[h].insert(ret.mat[h].end(), this->mat[i].begin() + 2, this->mat[i].end());
    } else if (column == n) {
      ret.mat[h].insert(ret.mat[h].end(), this->mat[i].begin() + 1, this->mat[i].end() - 1);
    } else {
      ret.mat[h].insert(ret.mat[h].end(), this->mat[i].begin() + 1, this->mat[i].begin() + column);
      ret.mat[h].insert(ret.mat[h].end(), this->mat[i].begin() + column + 1, this->mat[i].end());
    }
    ++h;
  }
  matrix true_ret(this->mat.size() - 2, this->mat[0].size() - 2);
  for (unsigned long i(1); i < m - 1; ++i) {
    for (unsigned long j(1); j < n; ++j) {
      true_ret[i][j] = ret[i][j];
    }
  }
  return true_ret;
}
matrix matrix::algebraic_cofactor(const unsigned long &row, const unsigned long &column) {
  return this->cofactor(row, column) * ((row + column) & 1 ? -1.00000000 : 1.00000000);
}
matrix matrix::principal_minor(const unsigned long &lines) {
  return this->cofactor(lines, lines);
}
#undef EPS
int main(int argc, char const *argv[]) {
  return 0;
}

Kirchhoff's theorem

这个定理可以用来求一个无向图\(G\)的生成树个数。首先明确几个概念: 1. \(G\)的度数矩阵\(D_G\)是一个\(n×n\)的矩阵, 并且满足: 当\(i\ne j\)时, \(d_{ij} = 0\); 当\(i = j\)时, \(d_{ij}\)等于\(v_i\)的度数。 2. \(G\)的邻接矩阵\(A_G\)也是一个\(n×n\)的矩阵, 并且满足: 如果\(v_i\)、\(v_j\)之间有边直接相连, 则\(a_{ij} = 1\), 否则为\(0\)。 定义\(G\)的Laplacian matrix\(C_G\)为\(C_G = D_G - A_G\), 则Kirchhoff's theorem可以描述为: > \(G\)的所有不同的生成树的个数等于其Laplacian matrix任何一个\(n - 1\)阶主子式的行列式的绝对值。

时间复杂度: \(\Theta\left(n^3\right)\)

题目链接: SPOJ 104 HIGH - Highways

这里只给主函数代码, 其余部分在上面的矩阵类里。 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
matrix graph;
int T, n, m, u, v;
int main(int argc, char const *argv[]) {
  scanf("%d", &T);
  while (T--) {
    scanf("%d %d", &n, &m);
    graph.assign(n, n);
    for (unsigned long i(1); i <= n; ++i) graph[i][i] = 0;
    while (m--) {
      scanf("%d %d", &u, &v);
      ++graph[u][u], ++graph[v][v];
      graph[u][v] = graph[v][u] = -1.00000000;
    }
    printf("%.0lf\n", n == 1 ? 1.00000000 : fabs(graph.principal_minor(1).det()));
  }
  return 0;
}