数学(数论)相关模板


一些乱七八糟的数学相关小东西。

数论相关

逆元

$\gcd \& \; \mathrm{exgcd}$

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
int gcd(int a, int b) {
// 欧几里得算法
while(b) {
std::swap(a, b);
b %= a;
}
return a;
}

int Stein(int a, int b) {
// Stein 算法,一定程度上是更相减损术的优化
int temp = 0;
while (true) {
if (b == 0 || a == b)
return a << c;
else if (a == 0)
return b << c;
while (!(a&1) && !(b&1))
a >>= 1,
b >>= 1,
++c;
while (a & 1 == 0) a >>= 1;
while (b & 1 == 0) b >>= 1;
if (a < b) std::swap(a, b);
a -= b;
}
}

int exgcd(int a, int b, int &x, int &y) {
if(!b) {
x = 1, y = 0;
return a;
}
int ans = exgcd(b, a % b, y, x);
y -= a / b * x;
return ans;
}

线性递推

1
2
3
ans[1] = 1;
for(int i = 2; i <= n; ++i)
ans[i] = (p - p / i * ans[p % i] % p) % p;

费马小定理

$$a^{p - 2} \equiv a^{-1} (\mathrm{mod} \ p)$$
直接用快速幂做即可,代码略。

素数筛法(线性筛)求积性函数

一份代码求了三个东西(素数筛,欧拉函数,莫比乌斯函数),三个愿望一次满足!

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
phi[1] = mu[1] = 1; //phi是欧拉函数,mu是莫比乌斯函数
for (int i = 2; i <= n; ++i) {
if (!vis[i]) {
vis[i] = true;
prime[++cnt] = i;
phi[i] = i - 1;
mu[i] = -1;
}
for (int j = 1; j <= cnt && prime[j] * i <= n; ++i) {
vis[i * prime[j]] = true;
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0;
phi[i * prime[j]] = phi[i] * prime[j];
break;
} else {
mu[i * prime[j]] = -mu[i];
phi[i * prime[j]] = phi[i] * (prime[j] - 1);
}
}
}

模方程(组)

线性模方程组:$\mathrm{CRT \ \& \ ExCRT}$

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
long long quickMul(long long a, long long b, long long p) {
//快(gui)速乘,防止中间过程爆 long long
long long res = 0;
while (b) {
if (b & 1)
res = (res + a) % p;
a = (a + a) % p;
b >>= 1;
}
return res;
}

long long CRT() {
long long res = 0, lcm = 1;
// x = a[i] (mod p[i]),所有 p[i] 互质
for (int i = 1; i <= n; ++i)
lcm *= p[i];
for (int i = 1; i <= n; ++i) {
long long x, y, mi = lcm / p[i];
exGCD(p[i], mi, x, y);
res = (res + mi * y * a[i] % lcm + lcm) % lcm;
}
return res % lcm;
}

long long exCRT() {
long long ans = a[1], lcm = p[1];
// x = a[i] (mod p[i])
for (int i = 2; i <= n; ++i) {
long long temp = (a[i] - ans % p[i] + p[i]) % p[i];
long long x, y, g = exGCD(lcm, p[i], x, y), p2 = p[i] / g;
if (temp % g) return -1;
x = quickMul(x, temp / g, p[i]);
ans += x * lcm;
lcm *= p2;
ans = (ans % lcm + lcm) % lcm;
}
return ans;
}

模意义下对数:$\mathrm{BSGS}$

$$b^x \equiv n (\mathrm{mod} \ p)$$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
long long BSGS() {
long long m = int(sqrt(p + 0.5));
long long inv = quickPow(b, p - m - 1);
std::map<int, int> s; s[1] = 0;
for (int i = 1; i < m; ++i) {
long long temp = quickPow(b, i);
if (temp == n) return i;
if (!s.count(temp)) s[temp] = i;
}
for (int i = 1; i <= m; ++i) {
n = n * inv % p;
if (s.count(n))
return s[n] + i * m;
}
return -1; //-1 表示无解
}

线性代数

高斯消元法

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
for (int i = 0; i < 115; ++i)
a[i] = new double[115];
for (int i = 1; i <= n; ++i)
for (int j = 1; j <= n + 1; ++j)
std::cin >> a[i][j];
for (int i = 1; i <= n; ++i) {
for (int j = i + 1; j <= n; ++j)
if (fabs(a[j][i]) - fabs(a[i][i]) > eps)
std::swap(a[i], a[j]);
if (fabs(a[i][i]) < eps) {
std::cout << "No Solution" << std::endl;
return 0;
}
for (int j = i + 1; j <= n + 1; ++j)
a[i][j] /= a[i][i];
a[i][i] = 1.0;
for (int j = i + 1; j <= n; ++j) {
double p = a[j][i];
for (int k = i; k <= n + 1; ++k)
a[j][k] -= a[i][k] * p;
}
}
ans[n] = a[n][n + 1];
for (int i = n - 1; i >= 1; --i) {
for (int j = i + 1; j <= n; ++j)
a[i][n + 1] -= ans[j] * a[i][j];
ans[i] = a[i][n + 1];
}

矩阵相关操作

矩阵乘法

此处实现一个矩阵结构体,并定义了输入输出和乘法操作。

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
struct matrix{
int n, m;
long long **val;
matrix(int a, int b): n(a), m(b){
val = new long long* [n];
for(int i = 0; i < n; ++i)
val[i] = new long long [m],
memset(val[i], 0, sizeof(long long)*m);
}
long long& operator()(int x, int y){
return val[x][y];
}
const long long& operator()(int x, int y) const{
return val[x][y];
}
friend istream& operator>>(istream &ist, matrix &ma){
for(int i = 0; i < ma.n; ++i)
for(int j = 0; j < ma.m; ++i)
ist>> ma(i, j);
return ist;
}
friend ostream& operator<<(ostream &ost, const matrix &ma){
for(int i = 0; i < ma.n; ++i, ost<<endl)
for(int j = 0; j < ma.m; ++i)
ost<<ma(i, j)<<' ';
return ost;
}
matrix operator*(const matrix& ma){
matrix res(n, ma.m);
for(int i = 0; i < n; ++i)
for(int j = 0; j < m; ++j)
for(int k = 0; k < ma.m; ++k)
res(i, k) += (*this)(i, j)*ma(j, k)%1000000007,
res(i, k) %= 1000000007;
return res;
}
void reset(int a, int b){
if(n){
for(int i = 0; i < n; ++i)
delete [] val[i];
delete [] val;
}
n = a, m = b;
val = new long long* [n];
for(int i = 0; i < n; ++i)
val[i] = new long long [m],
memset(val[i], 0, sizeof(long long)*m);
}
matrix& operator=(const matrix &ma){
//定义赋值运算符:避免拷贝得到的矩阵和源矩阵指向相同内存
reset(ma.n, ma.m);
for(int i = 0; i < n; ++i)
memcpy(val[i], ma.val[i], sizeof(long long)*m);
return *this;
}
~matrix(){//定义析构函数:避免内存泄漏
for(int i = 0; i < n; ++i)
delete [] val[i];
delete [] val;
}
};

多项式相关

多项式乘法:$\mathrm{FFT}$

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
/* 递归实现
* void FFT(std::complex<double> *f, int len, int mode) {
* if (len == 1) return;
* std::complex<double> *a1 = new std::complex<double> [len >> 1],
* *a2 = new std::complex<double> [len >> 1];
* for (int i = 0; i < len; i += 2)
* a1[i >> 1] = f[i], a2[i >> 1] = f[i + 1];
* FFT(a1, len >> 1, mode), FFT(a2, len >> 1, mode);
* std::complex<double> unit(cos(2 * pi / len), mode * sin(2 * pi / len)), temp(1, 0);
* for (int i = 0; i < len >> 1; ++i, temp *= unit)
* f[i] = a1[i] + temp * a2[i],
* f[i + (len >> 1)] = a1[i] - temp * a2[i];
* delete[] a1, delete[] a2;
* }
*/

// 迭代实现(推荐)
void FFT(std::complex<double> *f, int len, int mode) {
for (int i = 0; i < len; ++i)
if (i < ord[i]) std::swap(f[i], f[ord[i]]);
for (int i = 1; i < len; i <<= 1) {
std::complex<double> unit(cos(pi / i), mode * sin(pi / i));
for (int j = 0; j < len; j += i << 1) {
std::complex<double> temp(1, 0);
for (int k = 0; k < i; ++k, temp *= unit) {
std::complex<double> left(f[j + k]), right(f[j + i + k]);
f[j + k] = left + temp * right;
f[j + k + i] = left - temp * right;
}
}
}
}

int len = 1, temp = 0;
while (len <= n + m)
len <<= 1, ++temp;
for (int i = 0; i < len; ++i)
ord[i] = (ord[i >> 1] >> 1) | ((i & 1) << (temp - 1));
FFT(a, len, 1), FFT(b, len, 1);
for (int i = 0; i < len; ++i)
c[i] = a[i] * b[i];
FFT(c, len, -1);
for (int i = 0; i <= n + m; ++i)
std::cout << int(c[i].real() / len + 0.5) << ' ';

 评论