高斯消元+组合数

数学

一、高斯消元

时间复杂度：$O(n^3)$

运用场景：解n元方程组算法

实现方法：人脑模拟，构造上三角形，转化为最简阶梯型矩阵。

基本操作：

①把某一行乘一个非零的数；
②交换某两行；
③把某行的若干倍加到另一行。（初等行列变换）

实现方式：

①枚举每一列c;
②找到绝对值最大的一行；
③将该行移到最上面；
④将该行第一个数变为1；
⑤将下面的所有行的第c列减为0。

结果判断：

①完美阶梯型——唯一解
②左边零 $=$ 右边非零——无解
③左边零 $=$ 右边零——无穷多解

例题：解 $n$ 元方程组

代码实现：

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
const int N = 110;
const double eps = 1e-6;
int n;
double a[N][N];
int gauss() {
	int c, r;
	for (c = 0, r = 0; c < n; c++) {
		int t = r;
		for (int i = r; i < n; i++)
		    if (fabs(a[i][c]) > fabs(a[t][c]))
			    t = i;
		if (fabs(a[t][c]) < eps) continue;
		for (int i = c; i < n + 1; i++) swap(a[t][i], a[r][i]);
		for (int i = n; i >= c; i--) a[r][i] /= a[r][c];
		for (int i = r + 1; i < n; i++)
		    if (fabs(a[i][c]) > eps)
			    for (int j = n; j >= c; j--)
				    a[i][j] -= a[r][j] * a[i][c];
		r++;
	}
	if (r < n) {
		for (int i = r; i < n; i++)
		    if (fabs(a[i][n]) > eps)
		        return 2;
		return 1;
	}
	for (int i = n - 2; i >= 0; i--)
	    for (int j = i + 1; j < n; j++)
	        a[i][n] -= a[j][n] * a[i][j];
	return 0;
}
int main() {
	scanf("%d", &n);
	for (int i = 0; i < n; i++) 
	    for (int j = 0; j <= n; j++)
	        scanf("%lf", &a[i][j]);
	int t = gauss();
	if (t == 0) {
		for (int i = 0; i < n; i++) printf("%.2lf\n", a[i][n]);
	}
	else if (t == 1) puts("Infinite group solutions");
	else puts("No solution");
	return 0;
}

二、组合数（各种方式求组合数）

组合数

$C_a^b$ $=$ $\frac{a\ast (a-1)\ast ...\ast (a-b+1)}{1\ast 2\ast 3\ast ...\ast b}$ $=$ $\frac{a!}{b!(a-b)!}$

$C_b^a$ $=$ $C_{a-1}b$ $+$ $C_{a-1}^{b-1}$

1、询问次数大，询问范围小

求解方法：递推法

时间复杂度：O(n^2)

例题：询问次数 $1\ast 10^5$， 组合数大小：$1\le b\le a\le 2000$，模数固定。

代码实现：

#include <iostream>
#include <cstdio>
using namespace std;
const int N = 2010, mod = 1e9 + 7;
int n, c[N][N];
void init() {
	for (int i = 0; i < N; i++)
	    for (int j = 0; j <= i; j++)
	        if (!j) c[i][j] = 1;
	        else c[i][j] = (c[i-1][j] + c[i-1][j-1]) % mod;
}
int main() {
	init();
	scanf("%d", &n);
	for (int i = 0; i < n; i++) {
		int a, b;
		scanf("%d %d", &a, &b);
		printf("%d\n", c[a][b]);
	}
	return 0;
}

2、询问次数中等，询问范围中等

求解方法：预处理

时间复杂度：$O(nlogn)$

例题：询问次数 $1\ast 10^4$， $1\le b\le a\le 100000$，模数固定。

代码实现：

#include <iostream>
#include <cstdio>
using namespace std;
const long long N = 100010, mod = 1e9 + 7;
long long n, fact[N], infact[N];
long long qmi(long long a, long long k, long long p) {
	long long res = 1;
	while (k) {
		if (k & 1) res = res * a % p;
		a = a * a % p;
		k >>= 1;
	}
	return res;
}
int main() {
	fact[0] = infact[0] = 1;
	for (int i = 1; i < N; i++) {
		fact[i] = fact[i-1] * i % mod;
		infact[i] = infact[i-1] * qmi(i, mod - 2, mod) % mod;
	}
	scanf("%lld", &n);
	for (int i = 0; i < n; i++) {
		long long a, b;
		scanf("%lld %lld", &a, &b);
		printf("%lld\n", fact[a] * infact[b] % mod * infact[a-b] % mod);
	}
	return 0;
}

3、询问次数小，询问范围超大

求解方法：卢卡斯定理（Lucas）

时间复杂度：$O(PlogPlogN)$

Lucas定理：$C_a^b$ $\equiv$ $C_{amodp}^{bmodp}$ $\ast$ $C_{a/p}^{b/p}$ $(mod$ $p)$

例题：询问次数 $20$， 询问范围：$1\le b\le a\le b\le 10^{18}$，模数不固定。

代码实现：

#include <iostream>
#include <cstdio>
using namespace std;
long long qmi(long long a, long long k, long long p) {
	long long res = 1;
	while (k) {
		if (k & 1) res = res * a % p;
		a = a * a % p;
		k >>= 1;
	} 
	return res;
}
long long C(long long a, long long b, long long p) {
	if (b > a) return 0;
	long long res = 1;
	for (int i = 1, j = a; i <= b; i++, j--) {
		res = res * j % p;
		res = res * qmi(i, p - 2, p) % p;
	}
	return res;
}
long long lucas(long long a, long long b, long long p) {
	if (a < p && b < p) return C(a, b, p);
	return C(a % p, b % p, p) * lucas(a / p, b / p, p) % p;
}
int main() {
	int n;
	scanf("%d", &n);
	for (int i = 0; i < n; i++) {
		long long a, b, p;
		scanf("%lld %lld %lld", &a, &b, &p);
		printf("%lld\n", lucas(a, b, p));
	}
	return 0;
}

4、询问次数为 $1$， 询问范围 $1\le b\le a\le 100000$， 无模数。

求解方法：高精度 $+$ 因式分解

时间复杂度：$O(?)$

代码实现：

#include <iostream>
#include <cstdio>
#include <vector>
using namespace std;
const int N = 5010;
int primes[N], cnt;
int sum[N];
bool st[N];
void get_primes(int n) {
	for (int i = 2; i <= n; i++) {
		if (!st[i]) primes[cnt++] = i;
		for (int j = 0; primes[j] <= n / i; j++) {
			st[primes[j] * i] = true;
			if (i % primes[j] == 0) break;
		}
	}
}
int get(int n, int p) {
	int res = 0;
	while (n) {
		res += n / p;
		n /= p;
	}
	return res;
}
vector<int> mul(vector<int> a, int b) {
	vector<int> c;
	int t = 0;
	for (int i = 0; i < a.size(); i++) {
		t += a[i] * b;
		c.push_back(t % 10);
		t /= 10;
	}
	while (t) {
		c.push_back(t % 10);
		t /= 10;
	}
	return c;
}
int main() {
	int a, b;
	scanf("%d %d", &a, &b);
	get_primes(a);
	for (int i = 0; i < cnt; i++) {
		int p = primes[i];
		sum[i] = get(a, p) - get(a - b, p) - get(b, p);
	}
	vector<int> res;
	res.push_back(1);
	for (int i = 0; i < cnt; i++)
	    for (int j = 0; j < sum[i]; j++) 
	        res = mul(res, primes[i]);
	for (int i = res.size() - 1; i >= 0; i--) printf("%d", res[i]);
    puts("");
    return 0;
} 

三、小结

数学写题学思想，数学巩固推公式。