NOIp 2018 前的数学板子总结
数论
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\)的公因子。
代码:
template<typename IntegerType>
inline IntegerType Euclidean(const IntegerType &a, const IntegerType &b) {
return b ? Euclidean(b, a % b) : a;
}
或
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 + (a - \lfloor \frac{a}{b} \rfloor b)y \\
&= ay + (x - \lfloor \frac{a}{b} \rfloor y)b
\end {aligned}
\]
根据前面的结论: \(a\)和\(b\)都在减小, 当\(b\)减小到\(0\)时, 就可以得出\(x = 1\), \(y = 0\)。然后递归回去就可以求出最终的\(x\)和\(y\)了。
代码:
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)\)
代码:
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)\)
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 = \lfloor \frac{p}{i} \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} = (p - \lfloor \frac{p}{i} \rfloor)(p \bmod i)^{-1} \bmod p
\]
其中\(i \leq p\)。
时间复杂度: \(\Theta(a)\)
代码:
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 = \lfloor \frac{N}{m_i} \rfloor y\), 方程组等价于解同余方程: \(\lfloor \frac{N}{m_i} \rfloor y \equiv 1 \pmod {m_i}\), 若要得到特解\(y_i\), 只要令:\(x_i = \lfloor \frac{N}{m_i} \rfloor y_i\), 则方程组的解为:\(x_0 = a_1x_1 + a_2x_2 + \cdots + a_rx_r \pmod N\), 在模\(N\)的意义下唯一。
时间复杂度: \(\Theta(n \lg a)\)
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
#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(\sqrt{n})\)
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
又称素数的线性筛法。
方法如下:
- 使\(p\)等于\(2\), 即最小的素数
- 将表中所有\(p\)的倍数标记为合数
- 使\(p\)等于表中大于\(p\)的最小素数, 若没有则结束
- 重复第\(2\)步
时间复杂度: \(\Theta(n \lg \lg n)\)
题目链接: Luogu P3383 【模板】线性筛素数
#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质数测试
#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的线性筛选法
下面几个性质我就不证了知道有用就行
- \(\phi(p) = p - 1\)
- 如果\(i \bmod p = 0\), 那么\(\phi(ip) = p\phi(i)\)
- 若\(i \bmod p \ne 0\), 那么\(\phi(ip) = (p - 1)\phi(i)\)
时间复杂度: 趋近于\(\Theta(n)\)
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(\sqrt n)\)
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的结论:
-
\(Lucas(n, m, p) = cm(n \bmod p, m \bmod p) \cdot Lucas(\lfloor \frac{n}{p} \rfloor, \lfloor \frac{m}{p} \rfloor, p)\);
\(Lucas(x, 0, p) = 1\);
其中, \(cm(a, b) = a!(b!(a - b)!)^{p - 2} \bmod p = \frac{a!(b!)^{p - 2}}{(a - b)!}\)。
-
对于非负整数\(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} \]以下同余关系成立:
\[{\binom {m}{n}}\equiv \prod _{i=0}^{k}{\binom {m_{i}}{n_{i}}}{\pmod {p}} \]
时间复杂度: \(\Theta(\log_pn)\)
题目链接: Luogu P3807 【模板】卢卡斯定理
#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));
}
}
矩阵
#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\)的生成树个数。首先明确几个概念:
- \(G\)的度数矩阵\(D_G\)是一个\(n×n\)的矩阵, 并且满足: 当\(i\ne j\)时, \(d_{ij} = 0\); 当\(i = j\)时, \(d_{ij}\)等于\(v_i\)的度数。
- \(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(n^3)\)
题目链接: SPOJ 104 HIGH - Highways
这里只给主函数代码, 其余部分在上面的矩阵类里。
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;
}
