Skip to content

Commit d8890e6

Browse files
zimphazimpha
committed
add arbitrary modulo FFT, common operation for polynomial, replace ll with int64
1 parent ce6fc4d commit d8890e6

20 files changed

+445
-158
lines changed

‎mathematics/Fibonacci.cc

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,15 @@
1-
void fib(LL n, LL &x, LL &y) {// store in x, n-th
1+
#include "basic.hpp"
2+
3+
void fib(int64 n, int64 &x, int64 &y) {// store in x, n-th
24
if (n == 1) {
35
x = y = 1;
46
return;
57
} else if (n & 1) {
68
fib(n - 1, y, x);
79
y += x;
810
} else {
9-
LL a, b;
10-
fib(n >> 1, a, b, M);
11+
int64 a, b;
12+
fib(n >> 1, a, b);
1113
y = a * a + b * b;
1214
x = a * b + a * (b - a);
1315
}

‎mathematics/Simplex.cc

Lines changed: 12 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1,18 +1,23 @@
1+
#include <vector>
2+
#include <cmath>
3+
4+
using std::size_t;
5+
16
// 单纯形解线性规划 By 与星独白/猛犸也钻地 @ 2014.09.06
27
// 给出m个这样的约束条件：sum(A[i]*X[i])<=B
38
// 求出X的解，在满足X[i]>=0的情况下，sum(C[i]*X[i])达到最大
49
static const double EPS = 1e-12;
510
class Simplex {
611
public:
712
static const int SIZE = 100 + 10; // 略大于约束条件的最大数量
8-
vector<double> A[SIZE],B; // 约束条件：sum(A[k][i]*X[i])<=B[k]
9-
vector<double> C,X; // 目标函数：max(sum(C[i]*X[i]))
10-
void init(const vector<double>& C){ // 初始化，传入目标函数C
13+
std::vector<double> A[SIZE],B; // 约束条件：sum(A[k][i]*X[i])<=B[k]
14+
std::vector<double> C,X; // 目标函数：max(sum(C[i]*X[i]))
15+
void init(const std::vector<double>& C){ // 初始化，传入目标函数C
1116
this->C=C;
1217
for(size_t i=0;i<=B.size();i++) A[i].clear();
1318
B.clear(); X.clear();
1419
}
15-
void constraint(const vector<double>& u, double v){ // 添加约束条件
20+
void constraint(const std::vector<double>& u, double v){ // 添加约束条件
1621
A[B.size()]=u; // 约束条件的数组长度必须和目标函数的长度相同
1722
B.push_back(v);
1823
}
@@ -21,7 +26,7 @@ class Simplex {
2126
A[m].resize(n+1);
2227
for(int i=0;i<m;i++) A[i].push_back(B[i]);
2328
for(int i=0;i<n;i++) A[m][i]=-C[i];
24-
vector<int> idx(n+m);
29+
std::vector<int> idx(n+m);
2530
for(int i=0;i<n+m;i++) idx[i]=i;
2631
for(int k=0;k<m;k++) while(cmp(A[k][n],0)<0){
2732
int r=k,c;
@@ -52,14 +57,14 @@ class Simplex {
5257
}
5358
private:
5459
double cmp(double a, double b){return (1+fabs(a))*EPS<fabs(a-b)?a-b:0;}
55-
void pivot(vector<int>& idx, int r, int c){
60+
void pivot(std::vector<int>& idx, int r, int c){
5661
int m=B.size()+1,n=A[0].size();
5762
double div=1/A[r][c];
5863
for(int i=0;i<m;i++) if(i!=r)
5964
for(int j=0;j<n;j++) if(j!=c)
6065
A[i][j]-=A[r][j]*A[i][c]*div;
6166
for(int i=0;i<m;i++) A[i][c]*=-div;
6267
for(int j=0;j<n;j++) A[r][j]*=+div;
63-
A[r][c]=div; swap(idx[c],idx[n-1+r]);
68+
A[r][c]=div; std::swap(idx[c],idx[n-1+r]);
6469
}
6570
};

‎mathematics/ToneLLi-Shanks.cc

Lines changed: 0 additions & 27 deletions
This file was deleted.

‎mathematics/basic.hpp

Lines changed: 19 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -1,35 +1,31 @@
1-
#include <bits/stdc++.h>
2-
typedef long long LL;
1+
#pragma once
2+
#include <cassert>
3+
4+
using int64 = long long;
35

46
// mod should be not greater than 4e18
5-
inline LL mul_mod(LL a, LL b, LL mod) {
7+
// return a * b % mod when mod is less than 2^31
8+
inline int64 mul_mod(int64 a, int64 b, int64 mod) {
69
assert(0 <= a && a < mod);
710
assert(0 <= b && b < mod);
811
if (mod < int(1e9)) return a * b % mod;
9-
LL k = (LL)((long double)a * b / mod);
10-
LL res = a * b - k * mod;
12+
int64 k = (int64)((long double)a * b / mod);
13+
int64 res = a * b - k * mod;
1114
res %= mod;
1215
if (res < 0) res += mod;
1316
return res;
1417
}
1518

16-
inline LL add_mod(LL x, LL y) {
17-
return (x % y + y) % y;
19+
inline int64 add_mod(int64 x, int64 y, int64 mod) {
20+
return (x + y) % mod;
1821
}
1922

20-
LL pow_mod(LL a, LL n, LL m) {
21-
LL res = 1;
22-
for (a %= m; n; n >>= 1) {
23-
if (n & 1) res = res * a % m;
24-
a = a * a % m;
25-
}
26-
return res;
23+
inline int64 sub_mod(int64 x, int64 y, int64 mod) {
24+
return (x - y + mod) % mod;
2725
}
2826

29-
// 适合m超过int的时候
30-
// 条件允许的话使用__int128会比较快
31-
LL pow_mod_LL(LL a, LL n, LL m) {
32-
LL res = 1;
27+
int64 pow_mod(int64 a, int64 n, int64 m) {
28+
int64 res = 1;
3329
for (a %= m; n; n >>= 1) {
3430
if (n & 1) res = mul_mod(res, a, m);
3531
a = mul_mod(a, a, m);
@@ -43,7 +39,7 @@ T gcd(T a, T b) {
4339
}
4440

4541
// ax + by = gcd(a, b), |x| + |y| is minimum
46-
void exgcd(LL a, LL b, LL &g, LL &x, LL &y) {
42+
void exgcd(int64 a, int64 b, int64 &g, int64 &x, int64 &y) {
4743
if (!b) x = 1, y = 0, g = a;
4844
else {
4945
exgcd(b, a % b, g, y, x);
@@ -52,16 +48,16 @@ void exgcd(LL a, LL b, LL &g, LL &x, LL &y) {
5248
}
5349

5450
// ax = 1 (mod m), gcd(a, m) = 1
55-
LL mod_inv(LL a, LL m) {
56-
LL d, x, y;
51+
int64 mod_inv(int64 a, int64 m) {
52+
int64 d, x, y;
5753
exgcd(a, m, d, x, y);
5854
return d == 1 ? (x % m + m) % m : -1;
5955
}
6056

6157
//ax + by = c,x >= 0, x is minimum
6258
//xx = x + t * b, yy = y - t * a
63-
bool linear_equation(LL a, LL b, LL c, LL &x, LL &y) {
64-
LL g;
59+
bool linear_equation(int64 a, int64 b, int64 c, int64 &x, int64 &y) {
60+
int64 g;
6561
exgcd(a,b,g,x,y);
6662
if (c % g) return false;
6763
b /= g, a /= g, c /= g;

‎mathematics/binomial-coefficient.cc

Lines changed: 13 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,22 +1,25 @@
11
#include "basic.hpp"
2+
#include <vector>
23

34
namespace BinomialCoefficient {
4-
std::vector<int> fac, inv;
5-
void lucas_init(int m) {
6-
fac.assign(m, 1);
7-
inv.assign(m, 1);
8-
for (int i = 2; i < m; ++i) {
9-
fac[i] = 1ll * i * fac[i - 1] % m;
10-
inv[i] = 1ll * (m - m / i) * inv[m % i] % m;
5+
std::vector<int> fac, inv, ifac;
6+
void lucas_init(int n, int mod) {
7+
fac.assign(n + 1, 1);
8+
inv.assign(n + 1, 1);
9+
ifac.assign(n + 1, 1);
10+
for (int i = 2; i <= n; ++i) {
11+
fac[i] = (int64)i * fac[i - 1] % mod;
12+
inv[i] = int64(mod - mod / i) * inv[mod % i] % mod;
13+
ifac[i] = (int64)ifac[i - 1] * inv[i] % mod;
1114
}
1215
}
1316
inline int lucas_binom(int n, int m, int p) {
14-
return 1ll * fac[n] * inv[fac[m]] * inv[fac[n - m]] % p;
17+
return (int64)fac[n] * inv[fac[m]] * inv[fac[n - m]] % p;
1518
}
16-
int lucas(int n, int m, int p) {
19+
int lucas(int64 n, int64 m, int p) {
1720
if (!n && !m) return 1;
1821
int a = n % p, b = m % p;
1922
if (a < b) return 0;
20-
return 1ll * lucas(n / p, m / p, p) * lucas_binom(a, b, p) % p;
23+
return (int64)lucas(n / p, m / p, p) * lucas_binom(a, b, p) % p;
2124
}
2225
}

‎mathematics/congruence-equation.cc

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,10 @@
11
#include "basic.hpp"
2+
#include <vector>
23

34
// ax = b (mod m)
4-
std::vector<LL> congruence_equation(LL a, LL b, LL m) {
5-
std::vector<LL> ret;
6-
LL g = gcd(a, m), x;
5+
std::vector<int64> congruence_equation(int64 a, int64 b, int64 m) {
6+
std::vector<int64> ret;
7+
int64 g = gcd(a, m), x;
78
if (b % g != 0) return ret;
89
a /= g, b /= g;
910
x = mod_inv(a, m / g);

‎mathematics/congruence-equations.cc

Lines changed: 19 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -2,11 +2,11 @@
22

33
// 中国剩余定理， 要求m[i]互质
44
// x = c[i] (mod m[i]) (0 <= i < n)
5-
LL crt(int n, LL *c, LL *m) {
6-
LL M = 1, ans = 0;
5+
int64 crt(int n, int64 *c, int64 *m) {
6+
int64 M = 1, ans = 0;
77
for (int i = 0; i < n; ++i) M *= m[i];
88
for (int i = 0; i < n; ++i) {
9-
LL x, y, g, tm = M / m[i];
9+
int64 x, y, g, tm = M / m[i];
1010
exgcd(tm, m[i], g, x, y);
1111
ans = (ans + tm * x * c[i] % M) % M;
1212
}
@@ -16,15 +16,20 @@ LL crt(int n, LL *c, LL *m) {
1616
// 同余方程合并， m[i]可以不互质
1717
// x = c[i] mod m[i], c[i] < m[i]
1818
// 返回-1表示解不存在
19-
LL solve(int n, LL *c, LL *m){
20-
LL ans = c[0], LCM = m[0];
21-
for (int i = 1; i < n; ++i) {
22-
LL g, x, y;
23-
exgcd(LCM, m[i], g, x, y);
24-
if ((c[i] - ans) % g) return -1;
25-
LL tmp = add_mod((c[i] - ans) / g * x, m[i] / g);
26-
ans = add_mod(ans + LCM * tmp, LCM / g * m[i]);
27-
LCM = LCM / g * m[i];
28-
}
29-
return add_mod(ans, LCM);
19+
int64 solve(int n, int64 *c, int64 *m){
20+
auto mod = [](int64 x, int64 y) {
21+
x %= y;
22+
if (x < 0) x += y;
23+
return x;
24+
};
25+
int64 ans = c[0], LCM = m[0];
26+
for (int i = 1; i < n; ++i) {
27+
int64 g, x, y;
28+
exgcd(LCM, m[i], g, x, y);
29+
if ((c[i] - ans) % g) return -1;
30+
int64 tmp = mod((c[i] - ans) / g * x, m[i] / g);
31+
ans = mod(ans + LCM * tmp, LCM / g * m[i]);
32+
LCM = LCM / g * m[i];
33+
}
34+
return mod(ans, LCM);
3035
}

‎mathematics/det-mod.cc

Lines changed: 13 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,29 +1,32 @@
11
#include "basic.hpp"
2+
#include <algorithm>
23

34
/**
45
* 高斯消元法行列式求模.
56
* n为行列式大小, 计算|mat| % m
67
*
78
* time complexity: O(n^3 log n)
89
**/
9-
const int MAXN = 200;
10-
int det_mod(int n, int m, int mat[][MAXN]) {
11-
for (int i = 0; i < n; ++i)
12-
for (int j = 0; j < n; ++j)
13-
mat[i][j] %= m;
14-
LL ret = 1;
10+
const int N = 200;
11+
int det_mod(int n, int mod, int mat[][N]) {
12+
for (int i = 0; i < n; ++i) {
13+
for (int j = 0; j < n; ++j) {
14+
mat[i][j] %= mod;
15+
}
16+
}
17+
int64 ret = 1;
1518
for (int i = 0; i < n; ++i) {
1619
for (int j = i + 1; j < n; ++j)
1720
for (; mat[j][i]; ret = -ret) {
18-
LL t = mat[i][i] / mat[j][i];
21+
int64 t = mat[i][i] / mat[j][i];
1922
for (int k = i; k < n; ++k) {
20-
mat[i][k] = (mat[i][k] - mat[j][k] * t) % m;
23+
mat[i][k] = (mat[i][k] - mat[j][k] * t) % mod;
2124
std::swap(mat[j][k], mat[i][k]);
2225
}
2326
}
2427
if (mat[i][i] == 0) return 0;
25-
ret = ret * mat[i][i] % m;
28+
ret = ret * mat[i][i] % mod;
2629
}
27-
if (ret < 0) ret += m;
30+
if (ret < 0) ret += mod;
2831
return (int)ret;
2932
}

0 commit comments

Comments
 (0)