数论
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 mod b)), 所以
[egin {aligned}
ax + by &= (a, b) \
&= (b, a mod b) \
&= bx + (a mod b)y \
&= bx + (a - lfloor frac{a}{b}
floor b)y \
&= ay + (x - lfloor frac{a}{b}
floor 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} floor), (k = p mod 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}
floor)(p mod i)^{-1} mod 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_2cdots m_r), 则同余方程组:
[egin{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)):
[egin{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} floor y), 方程组等价于解同余方程: (lfloor frac{N}{m_i} floor y equiv 1 pmod {m_i}), 若要得到特解(y_i), 只要令:(x_i = lfloor frac{N}{m_i} floor 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))都互素, 这时候就需要用到扩展中国剩余定理。
对于同余方程组:
[egin{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}
]
我们首先只考虑其中的两个方程, 可以得到
[egin {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) mid (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
", 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) e 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))的解为(egin{cases}x_1 = 1 \ x_2 = p - 1end{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 mod p = 0), 那么(phi(ip) = pphi(i))
- 若(i mod p e 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} mod p)的值。其中: (n)和(m)时非负整数, (p)是素数。
Lucas's theorem的结论:
-
(Lucas(n, m, p) = cm(n mod p, m mod p) cdot Lucas(lfloor frac{n}{p} floor, lfloor frac{m}{p} floor, p));
(Lucas(x, 0, p) = 1);
其中, (cm(a, b) = a!(b!(a - b)!)^{p - 2} mod 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} ]以下同余关系成立:
[{inom {m}{n}}equiv prod _{i=0}^{k}{inom {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
", 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 ? '
' : ' ');
}
}
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 e 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 == 1 ? 1.00000000 : fabs(graph.principal_minor(1).det()));
}
return 0;
}