1.快速幂
O ( l o g n ) O(logn) O(logn)
ll ksm(ll a, ll b, ll p) {
ll res = 1;
a % = p;
while (b) {
if (b & 1) res = res * a % p;
b >>= 1;
a = (a * a) % p;
}
return res;
}
2.矩阵快速幂
P3390 【模板】矩阵快速幂
给定
n
×
n
n\times n
n×n 的矩阵
A
A
A,求
A
k
A^k
Ak。
O ( N 3 l o g n ) O(N^3logn) O(N3logn)
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll INF = 1e18;
const ll N = 1e2 + 10;
const ll mod = 1e9 + 7;
ll n, k;
struct matrix {
ll m[N][N];
}a;
matrix operator*(const matrix& a, const matrix& b) {
matrix c;
memset(c.m, 0, sizeof(c.m));
for (int i = 1; i < N; ++i) {
for (int j = 1; j < N; ++j) {
for (int k = 1; k < N; ++k) {
c.m[i][j] =(c.m[i][j]+ a.m[i][k] * b.m[k][j])%mod;
}
}
}
return c;
}
matrix matrixksm(matrix a, ll k) {
matrix ans;
memset(ans.m, 0, sizeof(ans.m));
for (int i = 1; i < N; ++i) ans.m[i][i] = 1;
while (k) {
if (k & 1) ans =ans* a;
a = a * a;
k >>= 1;
}
return ans;
}
void solve() {
cin >> n>>k;
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= n; ++j) {
cin >> a.m[i][j];
}
}
a=matrixksm(a, k);
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= n; ++j) {
cout << a.m[i][j] << " ";
}
cout << '\n';
}
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin>>_t;
while (_t--) {
solve();
}
return 0;
}
3.矩阵递推
P1939 【模板】矩阵加速(数列)
已知一个数列
a
a
a,它满足:
a
x
=
{
1
x
∈
{
1
,
2
,
3
}
a
x
−
1
+
a
x
−
3
x
≥
4
a_x= \begin{cases} 1 & x \in\{1,2,3\}\\ a_{x-1}+a_{x-3} & x \geq 4 \end{cases}
ax={1ax−1+ax−3x∈{1,2,3}x≥4
求
a
a
a 数列的第
n
n
n 项对
1
0
9
+
7
10^9+7
109+7 取余的值。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll INF = 1e18;
const ll N = 1e5 + 10;
const ll mod = 1e9 + 7;
struct matrix {
ll m[4][4];
};
matrix operator*(const matrix& a, const matrix& b) {
matrix c;
memset(c.m, 0, sizeof(c.m));
for (int i = 1; i <= 3; ++i) {
for (int j = 1; j <= 3; ++j) {
for (int k = 1; k <= 3; ++k) {
c.m[i][j] = (c.m[i][j] + a.m[i][k] * b.m[k][j]) % mod;
}
}
}
return c;
}
matrix ksm(matrix a, ll k) {
matrix ans;
memset(ans.m, 0, sizeof(ans.m));
for (int i = 1; i <= 3; ++i) ans.m[i][i] = 1;
while (k) {
if (k & 1) ans = ans * a;
a = a * a;
k >>= 1;
}
return ans;
}
void solve() {
int n;
cin >> n;
if (n == 1 || n == 2 || n == 3) {
cout << 1 << '\n';
return;
}
matrix s,a;
memset(s.m, 0, sizeof(s.m));
memset(a.m, 0, sizeof(a.m));
a.m[1][1] = a.m[1][2] = a.m[1][3] = 1;
s.m[1][1] = s.m[1][2] = s.m[2][3] = s.m[3][1] = 1;
a =a* ksm(s, n - 3);
cout << a.m[1][1] << '\n';
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
//_t = 1;
cin>>_t;
while (_t--) {
solve();
}
return 0;
}
4.矩阵路径问题
两点间只经过n条边的总路径数量
计算邻接矩阵的幂
G
=
M
n
G=M^n
G=Mn,其元素
G
(
i
,
j
)
G(i,j)
G(i,j)的值是从点i到点j经过n条边(或n-1个点)的总路径数量
两点间只经过n条边的最短路径长度
用邻接矩阵M表示一个图,
M
(
i
,
j
)
M(i,j)
M(i,j)是其第i行第j列元素,表示点i和点j的直连关系,若点i和点j直连,令
M
[
i
,
j
]
M[i,j]
M[i,j]等于边
i
j
ij
ij的权值,否则等于无穷大;当i=j时,令
M
[
i
,
i
]
=
∞
M[i,i]=\infty
M[i,i]=∞。定义矩阵的广义乘法
M
×
M
=
min
a
=
1
N
[
M
(
i
,
a
)
+
M
(
a
,
j
)
]
M\times M=\min_{a=1}^N[M(i,a)+M(a,j)]
M×M=mina=1N[M(i,a)+M(a,j)],也就是把普通的矩阵乘法从求和改成了取最小值,把内部项相乘改成了相加。
计算邻接矩阵的广义幂
G
=
M
n
G=M^n
G=Mn,
G
(
i
,
j
)
G(i,j)
G(i,j)的值是从i到j经过n条边(或n-1个点)的最短路径长度。
5.高斯消元
P3389 【模板】高斯消元法
给定一个线性方程组,对其求解
O ( n 3 ) O(n^3) O(n3)
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 1e2 + 10;
double a[N][N];
void solve() {
int n;
cin >> n;
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= n+1; ++j) {
cin >> a[i][j];
}
}
for (int i = 1; i <= n; ++i) {
int maxn = i;
for (int j = i+1; j <= n; ++j) {
if (fabs(a[j][i]) > fabs(a[maxn][i])) maxn = j;
}
for (int j = 1; j <= n + 1; ++j) swap(a[i][j], a[maxn][j]);
if (fabs(a[i][i]) < eps) {
cout << "No Solution";
return;
}
for (int j = n + 1; j >= i; --j) a[i][j]/=a[i][i];
for (int j = 1; j <= n; ++j) {
if (j != i) {
for (int k = n+1; k >=i; --k) a[j][k] -= a[j][i] * a[i][k];
}
}
}
for (int i = 1; i <= n; ++i) cout << setiosflags(ios::fixed)<<setprecision(2)<<a[i][n+1] << '\n';
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin>>_t;
while (_t--) {
solve();
}
return 0;
}
高斯消元求异或方程组
poj1681 Painter’s Problem
有一面方形的墙,是由n*n块小方砖砌成的。有些砖是白色的,有些砖是黄色的。鲍勃是个油漆工,他想把所有的砖都漆成黄色。但是鲍勃的刷子有点问题。一旦他用这把刷子刷砖块(i, j), (i, j), (i+1, j), (i, j-1), (i, j-1)和(i, j+1)处的砖块都会改变颜色。你的任务是找到鲍勃应该涂的最少砖块数量,以便使所有砖块都变成黄色。
AcWing884高斯消元解异或线性方程组
输入一个包含 n个方程 n个未知数的异或线性方程组。
方程组中的系数和常数为0或1,每个未知数的取值也为0或1,求解这个方程组。
异或线性方程组示例如下
M[1][1]x[1] ^ M[1][2]x[2]……M[1][n]x[n]
M[2][1]x[1] ^ M[2][2]x[2]……M[2][n]x[n]
……
M[n][1]x[1] ^ M[n][2]x[2]……M[n][n]x[n]
其中^表示异或(XOR),M[i][j] 表示第i个式子中x[j]的系数,B[i]是第i个方程右端的常数,取值均为0或1。
#include<bits/stdc++.h>
using namespace std;
const int N = 110;
int n,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 (a[i][c])t = i;
}
if (!a[t][c]) continue;
for (int i = c; i <= n; i++) swap(a[r][i], a[t][i]);
for (int i = r + 1; i < n; i++) {
if (a[i][c]) {
for (int j = n; j >= c; j--) a[i][j] ^= a[r][j];
}
}
r++;
}
if (r < n){
for (int i = r; i < n; i++) {
if (a[i][n]) return 2;
}
return 1;
}
for (int i = n - 1; i >= 0; i--) {
for (int j = i + 1; j < n; j++) a[i][n] ^= a[i][j] * a[j][n];
}
return 0;
}
int main(){
cin >> n;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n + 1; j++) cin >> a[i][j];
}
int t = gauss();
if (t == 0){
for (int i = 0; i < n; i++) cout << a[i][n] << endl;
}
else if (t == 1) puts("Multiple sets of solutions");
else puts("No solution");
return 0;
}
高斯消元求取模方程组
hdu5755 Gambler Bo
赌徒博精通矩阵游戏。
你有一个N×M矩阵,每个单元格都有一个{0,1,2}的值。
在这个游戏中,你可以在矩阵中选择一个单元格,这个单元格加2,所有相邻的单元格加1。
例如,您选择单元格(x,y),则(x,y)的值将被加2,(x - 1,y)(x,y - 1)(x,y - 1)(x,y+1)的值将被加1。
如果你选择单元格(1,2),单元格(1,2)将被加2,单元格(2,2)(1,1)(1,3)将被加1,单元格(0,2)不会改变,因为它不在矩阵中。
如果某些单元格的值超过2,则这些值将取3的模。
赌徒Bo给了你这样一个矩阵,你的任务是通过不超过2NM的操作使这个矩阵的所有值为0。
矩阵求逆
6.线性基
最大异或和
P3812 【模板】线性基
给定
n
n
n 个整数(数字可能重复),求在这些数中选取任意个,使得他们的异或和最大。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 1e5 + 10;
ll d[60];
void insert(ll x) {
for (int i = 50; i >= 0; --i) {
if (x & (1ll << i)) {
if (d[i]) x ^= d[i];
else {
d[i] = x;
return;
}
}
}
}
void solve() {
int n;
cin >> n;
for (int i = 1; i <= n; ++i) {
ll a;
cin >> a;
insert(a);
}
ll ans = 0;
for (int i = 50; i >= 0; --i) {
if ((ans ^ d[i]) > ans) ans ^= d[i];
}
cout << ans;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin>>_t;
while (_t--) {
solve();
}
return 0;
}
最小异或和
k大/k小异或和
7.0/1分数规划
∑
i
=
1
n
a
i
s
i
∑
i
=
1
n
b
i
s
i
≥
x
\frac{\sum\limits_{i=1}^na_is_i}{\sum\limits_{i=1}^nb_is_i}\ge x
i=1∑nbisii=1∑naisi≥x
f
=
∑
i
=
1
n
(
a
i
−
x
b
i
)
s
i
≥
0
f=\sum\limits_{i=1}^n(a_i-xb_i)s_i\ge0
f=i=1∑n(ai−xbi)si≥0
基础
poj2976 Dropping tests
在一门课程中,你要做n次测试。如果你在测试1中答对了2道题中的1道,你的累积平均值定义为
100
∑
i
=
1
n
a
i
∑
i
=
1
n
b
i
100\frac{\sum_{i=1}^na_i}{\sum_{i=1}^nb_i}
100∑i=1nbi∑i=1nai
给定你的考试成绩和一个正整数k,确定如果允许你放弃考试成绩中的任意k,你的累积平均成绩能达到多高。
假设你参加了3次考试,分数分别是5/5、0/1和2/6。在不放弃任何考试的情况下,你的累积平均成绩是
100
5
+
0
+
2
5
+
1
+
6
=
50
100\frac{5+0+2}{5+1+6}=50
1005+1+65+0+2=50。但是,如果你放弃了第三次考试,你的累积平均分就会变成
100
5
+
0
5
+
1
≈
83.33
≈
83
100\frac{5+0}{5+1}\approx83.33\approx83
1005+15+0≈83.33≈83。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 1e3 + 10;
struct node {
int a, b;
double y;
bool operator<(const node a)const {
return y > a.y;
}
}p[N];
int n, k;
bool check(double x) {
for (int i = 1; i <= n; ++i) p[i].y = 1.0 * p[i].a - x * p[i].b;
sort(p + 1, p + n + 1);
double f = 0;
for (int i = 1; i <= k; ++i) f += p[i].y;
return f < 0;
}
void solve() {
while (cin >> n >> k && n + k) {
double l = 0, r = 0;
k = n - k;
for (int i = 1; i <= n; ++i) cin >> p[i].a, r += p[i].a;
for (int i = 1; i <= n; ++i) cin >> p[i].b;
for (int i = 1; i <= 50; ++i) {
double mid = (l + r) / 2;
if (check(mid)) r = mid;
else l = mid;
}
cout << (int)(100 * (l + 0.005)) << '\n';
}
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin>>_t;
while (_t--) {
solve();
}
return 0;
}
最优比率背包
P4377 [USACO18OPEN] Talent Show G
Farmer John 要带着他的
n
n
n 头奶牛,方便起见编号为
1
…
n
1\ldots n
1…n,到农业展览会上去,参加每年的达牛秀!他的第
i
i
i 头奶牛重量为
w
i
w_i
wi,才艺水平为
t
i
t_i
ti,两者都是整数。
在到达时,Farmer John 就被今年达牛秀的新规则吓到了:
(一)参加比赛的一组奶牛必须总重量至少为
W
W
W(这是为了确保是强大的队伍在比赛,而不仅是强大的某头奶牛),并且。
(二)总才艺值与总重量的比值最大的一组获得胜利。
FJ 注意到他的所有奶牛的总重量不小于
W
W
W,所以他能够派出符合规则(一)的队伍。帮助他确定这样的队伍中能够达到的最佳的才艺与重量的比值。
最优比率生成树
畅通工程
有 n 个城市,城市之间道路,现在请你建造 n−1 条道路使城市连通。但不是任意两个城市之间都可以建造道路,必须两个城市之间相互友好才可以建造道路(这种相互友好的关系不具有传递性)。根据你的了解,有 m 对城市之间相互友好,可以建造道路,并且预期建造的道路长度为 len, 建造这条道路的花费为 cost 。你希望建造的所有道路的长度之和与花费之和的比值最大。
最优比率环
P2868 [USACO07DEC]Sightseeing Cows G
作为对奶牛们辛勤工作的回报,
F
a
r
m
e
r
J
o
h
n
Farmer\ John
Farmer John决定带她们去附近的大城市玩一天。旅行的前夜,奶牛们在兴奋地讨论如何最好地享受这难得的闲暇。
很幸运地,奶牛们找到了一张详细的城市地图,上面标注了城市中所有
L
(
2
⩽
L
⩽
1000
)
L(2\leqslant L\leqslant1000)
L(2⩽L⩽1000)座标志性建筑物(建筑物按
1
…
L
1\dots L
1…L顺次编号),以及连接这些建筑物的
P
(
2
⩽
P
⩽
5000
)
P(2\leqslant P\leqslant5000)
P(2⩽P⩽5000)条道路。按照计划,那天早上
F
a
r
m
e
r
J
o
h
n
Farmer\ John
Farmer John会开车将奶牛们送到某个她们指定的建筑物旁边,等奶牛们完成她们的整个旅行并回到出发点后,将她们接回农场。由于大城市中总是寸土寸金,所有的道路都很窄,政府不得不把它们都设定为通行方向固定的单行道。
尽管参观那些标志性建筑物的确很有意思,但如果你认为奶牛们同样享受穿行于大城市的车流中的话,你就大错特错了。与参观景点相反,奶牛们把走路定义为无趣且令她们厌烦的活动。对于编号为
i
i
i的标志性建筑物,奶牛们清楚地知道参观它能给自己带来的乐趣值
F
i
(
1
⩽
F
i
⩽
1000
)
F_i (1\leqslant F_i\leqslant1000)
Fi(1⩽Fi⩽1000)。相对于奶牛们在走路上花的时间,她们参观建筑物的耗时可以忽略不计。
奶牛们同样仔细地研究过城市中的道路。她们知道第i条道路两端的建筑物
L
1
i
L1_i
L1i和
L
2
i
L2_i
L2i(道路方向为
L
1
i
→
L
2
i
L1_i \rightarrow L2_i
L1i→L2i),以及她们从道路的一头走到另一头所需要的时间
T
i
(
1
⩽
T
i
⩽
1000
)
T_i(1\leqslant T_i\leqslant1000)
Ti(1⩽Ti⩽1000)。
为了最好地享受她们的休息日,奶牛们希望她们在一整天中平均每单位时间内获得的乐趣值最大。当然咯,奶牛们不会愿意把同一个建筑物参观两遍,也就是说,虽然她们可以两次经过同一个建筑物,但她们的乐趣值只会增加一次。顺便说一句,为了让奶牛们得到一些锻炼,
F
a
r
m
e
r
J
o
h
n
Farmer\ John
Farmer John要求奶牛们参观至少
2
2
2个建筑物。
请你写个程序,帮奶牛们计算一下她们能得到的最大平均乐趣值。
最大密度子图
AcWing2324 生活的艰辛
约翰是一家公司的 CEO。公司的股东决定让他的儿子斯科特成为公司的经理。约翰十分担心,儿子会因为在经理岗位上表现优异而威胁到他 CEO 的位置。因此,他决定精心挑选儿子要管理的团队人员,让儿子知道社会的险恶。
已知公司中一共有 n 名员工,员工之间共有 m 对两两矛盾关系。如果将一对有矛盾的员工安排在同一个团队,那么团队的管理难度就会增大。一个团队的管理难度系数等于团队中的矛盾关系对数除以团队总人数。团队的管理难度系数越大,团队就越难管理。
约翰希望给儿子安排的团队的管理难度系数尽可能大。请帮帮他。
8.GCD与LCM
gcd(a,b)=gcd(a,ka+b)
gcd(ka,kb)=k*gcd(a,b)
若gcd(a,b)=d,gcd(a/d,b/d)=1
算术基本定理:任何大于1的正整数n都可以唯一分解为有限给素数乘积
lcm(a,b)=a/gcd(a,b)*b
9.裴蜀定理
设 a,b 是不全为零的整数,则存在整数 x,y, 使得 ax+by=gcd(a,b).
P4549 【模板】裴蜀定理
给定一个包含
n
n
n 个元素的整数序列
A
A
A,记作
A
1
,
A
2
,
A
3
,
.
.
.
,
A
n
A_1,A_2,A_3,...,A_n
A1,A2,A3,...,An。
求另一个包含
n
n
n 个元素的待定整数序列
X
X
X,记
S
=
∑
i
=
1
n
A
i
×
X
i
S=\sum\limits_{i=1}^nA_i\times X_i
S=i=1∑nAi×Xi,使得
S
>
0
S>0
S>0 且
S
S
S 尽可能的小。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 1e5 + 10;
int a[30];
int gcd(int a, int b) {
return abs(b == 0 ? a : gcd(b, a % b));
}
void solve() {
int n;
cin >> n;
for (int i = 1; i <= n; ++i) cin >> a[i];
int ans=a[1];
for (int i = 2; i <= n; ++i) {
ans = gcd(ans, a[i]);
}
cout << ans;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin>>_t;
while (_t--) {
solve();
}
return 0;
}
10. 扩展欧几里得
二元线性丢番图方程
P5656 【模板】二元一次不定方程 (exgcd)
给定不定方程
a
x
+
b
y
=
c
ax+by=c
ax+by=c
若该方程无整数解,输出
−
1
-1
−1。
若该方程有整数解,且有正整数解,则输出其正整数解的数量,所有正整数解中
x
x
x 的最小值,所有正整数解中
y
y
y 的最小值,所有正整数解中
x
x
x 的最大值,以及所有正整数解中
y
y
y 的最大值。
若方程有整数解,但没有正整数解,你需要输出所有整数解中
x
x
x 的最小正整数值,
y
y
y 的最小正整数值。
正整数解即为
x
,
y
x, y
x,y 均为正整数的解,
0
\boldsymbol{0}
0 不是正整数。
整数解即为
x
,
y
x,y
x,y 均为整数的解。
x
x
x 的最小正整数值即所有
x
x
x 为正整数的整数解中
x
x
x 的最小值,
y
y
y 同理。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll INF = 1e18;
const ll N = 1e5 + 10;
inline ll read() {
ll res = 0;
char ch = getchar();
while (ch < '0' || ch>'9') ch = getchar();
while (ch >= '0' && ch <= '9') { res = (res << 3) + (res << 1) + (ch ^ 48); ch = getchar(); }
return res;
}
ll exgcd(ll a, ll b, ll& x, ll& y) {
if (b == 0) {
x = 1, y = 0;
return a;
}
ll t = exgcd(b, a % b, x, y);
ll x0 = x, y0 = y;
x = y0;
y = x0 - (a / b) * y0;
return t;
}
void solve() {
ll a, b, c,x,y,t;
a = read(), b = read(), c = read();
t = exgcd(a, b, x, y);
if (c % t) cout << -1 << '\n';
else {
x *= (c / t), y *= (c / t);
ll minx, miny, maxx, maxy;
minx = x > 0 && x % (b/t) != 0 ? x % (b/t) : x % (b/t) + b/t;
maxy = (c - minx * a) / b;
miny = y > 0 && y % (a/t) != 0 ? y % (a/t) : y % (a/t) + a/t;
maxx = (c - miny * b) / a;
if (maxy > 0) {
ll num = (maxx - minx) / (b / t) + 1;
cout << num << " " << minx << " " << miny << " " << maxx << " " << maxy << '\n';
}
else cout << minx << " " << miny << '\n';
}
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
//_t = 1;
_t = read();
while (_t--) {
solve();
}
return 0;
}
多元线性丢番图方程
分解为n-1个二元方程,从后往前依次求解
11. 同余方程(扩欧求逆)
P1082 [NOIP2012 提高组] 同余方程
求关于
x
x
x 的同余方程
a
x
≡
1
(
m
o
d
b
)
a x \equiv 1 \pmod {b}
ax≡1(modb) 的最小正整数解。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
void exgcd(ll a, ll b, ll& x, ll& y) {
if (b == 0) {
x = 1;
y = 0;
return;
}
exgcd(b, a%b, x, y);
ll x0 = x, y0 = y;
x = y0;
y = x0 - (a / b) * y0;
}
void solve() {
ll a, b, x, y;
cin >> a >> b;
exgcd(a, b, x, y);
cout << (x % b + b) % b;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin >> _t;
while (_t--) {
solve();
}
return 0;
}
12.费马小定理
若p为素数,gcd(a,p)=1,则 a p − 1 ≡ 1 ( m o d p ) a^{p-1}\equiv1(mod\ p) ap−1≡1(mod p)
13.中国剩余定理
P1495 【模板】中国剩余定理(CRT)/ 曹冲养猪
自从曹冲搞定了大象以后,曹操就开始捉摸让儿子干些事业,于是派他到中原养猪场养猪,可是曹冲满不高兴,于是在工作中马马虎虎,有一次曹操想知道母猪的数量,于是曹冲想狠狠耍曹操一把。举个例子,假如有
16
16
16 头母猪,如果建了
3
3
3 个猪圈,剩下
1
1
1 头猪就没有地方安家了。如果建造了
5
5
5 个猪圈,但是仍然有
1
1
1 头猪没有地方去,然后如果建造了
7
7
7 个猪圈,还有
2
2
2 头没有地方去。你作为曹总的私人秘书理所当然要将准确的猪数报给曹总,你该怎么办?
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 1e5 + 10;
ll m[11], a[11],M[11],invM[11];
void exgcd(ll a, ll b, ll& x, ll& y) {
if (b == 0) {
x = 1, y = 0;
return ;
}
exgcd(b, a % b, x, y);
ll x0 = x, y0 = y;
x = y0;
y = x0 - (a / b) * y0;
}
void solve() {
int n;
cin >> n;
ll mul=1;
for (int i = 1; i <= n; ++i) {
cin >> m[i] >> a[i];
mul *= m[i];
}
ll ans = 0;
for (int i = 1; i <= n; ++i) {
M[i] = mul / m[i];
ll x , y;
exgcd(M[i], m[i], x, y);
ans = (ans + a[i] * M[i] * ((x + m[i]) % m[i])) % mul;
}
cout << ans;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin>>_t;
while (_t--) {
solve();
}
return 0;
}
14.扩展中国剩余定理
P4777 【模板】扩展中国剩余定理(EXCRT)
给定
n
n
n 组非负整数
a
i
,
b
i
a_i, b_i
ai,bi ,求解关于
x
x
x 的方程组的最小非负整数解。
{
x
≡
b
1
(
m
o
d
a
1
)
x
≡
b
2
(
m
o
d
a
2
)
.
.
.
x
≡
b
n
(
m
o
d
a
n
)
\begin{cases} x \equiv b_1\ ({\rm mod}\ a_1) \\ x\equiv b_2\ ({\rm mod}\ a_2) \\ ... \\ x \equiv b_n\ ({\rm mod}\ a_n)\end{cases}
⎩
⎨
⎧x≡b1 (mod a1)x≡b2 (mod a2)...x≡bn (mod an)
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 1e5 + 10;
ll m[N], a[N];
ll mul(ll a, ll b, ll p) {
ll res = 0;
while (b) {
if (b & 1) res = (res + a) % p;
a = (a + a) % p;
b >>= 1;
}
return res;
}
ll exgcd(ll a, ll b, ll& x, ll& y) {
if (b == 0) { x = 1; y = 0; return a;}
ll t = exgcd(b, a % b, x, y);
ll x0 = x, y0 = y;
x = y0, y = x0 - (a / b) * y0;
return t;
}
void solve() {
int n;
cin >> n;
for (int i = 1; i <= n; ++i) cin >> m[i] >> a[i];
ll m1 = m[1], a1 = a[1],ans=0;
for (int i = 2; i <= n; ++i) {//合并每两个等式
ll m2 = m[i], a2 = a[i];
ll a = m1, b = m2, c = (a2 - a1 % m2 + m2) % m2;
ll x, y;//求解aX+bY=c
ll t = exgcd(a, b, x, y);
x = mul(x, c / t, b / t);
ans = a1 + x * m1;//代回
m1 = m1 / t * m2;//合并后新m1
ans = (ans % m1 + m1) % m1;//最小正整数解
a1 = ans;//合并后新a1
}
cout << ans;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin>>_t;
while (_t--) {
solve();
}
return 0;
}
15. 素数筛
埃氏筛
O ( n l o g l o g n ) O(nloglogn) O(nloglogn)
int ans[N],cnt = 0;
bool tag[N]
int shai() {
for(int i = 2; i <= N; i++){
if(!tag[i]){
ans[++cnt]=i;
for(int j = 2 * i; j < N; j += i) tag[j] = 1;
}
}
}
欧拉筛
O ( n ) O(n) O(n)
int ans[N],cnt = 0;
bool tag[N];
void shai() {
for (int i = 2; i <= N; i++) {
if (!tag[i]) ans[++cnt] = i;
for (int j = 1; j <= cnt && i * ans[j] <= N; j++) {
tag[i * ans[j]] = 1;
if (i % ans[j] == 0) break;
}
}
return;
}
区间素数筛
求[a,b]内素数
2~sqrt(b)的倍数可以最多筛掉到b以内的所有非素数
bool tag1[N],tag2[N];//tag1[0,sqrt(b)],tag2[a,b]
ll shai(ll a, ll b){
for (ll i = 2; i * i <= b; i++){
if (!tag1[i]){
for (ll j = 2 * i; j * j <= b; j += i) tag1[j] = 1;//低区间筛选
for (ll j = max(2ll, (a - 1) / i + 1) * i; j <= b; j += i) tag2[j - a] = 1;//跳到大于a的值来筛选a~b
}
}
ll ans = 0;
for (int i = 0; i <= b - a; i++) {
if (!tag2[i]) ans++;
}
if (a == 1) ans--;
return ans;
}
16.米勒拉宾素性测试
hdu2138 How many prime numbers
O
(
7
l
o
g
3
n
)
O(7log^3n)
O(7log3n)
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 1e5 + 10;
ll n, a[7] = { 2,325,9375,28178,450775,9780504,1795265022 };
ll mul(ll a, ll b, ll p) {
a %= p, b %= p;
ll res = 0;
while (b) {
if (b & 1) {
res += a;
if (res >= p) res -= p;
}
a <<= 1;
if (a >= p) a -= p;
b >>= 1;
}
return res;
}
ll ksm(ll a, ll b, ll p) {
ll res = 1;
a %= p;
while (b) {
if (b & 1) res = mul(res,a,p) % p;
a = mul(a,a,p) % p;
b >>= 1;
}
return res;
}
bool miller_rabin(ll m) {
if (m == 2) return 1;
if (m <= 1 || m % 2 == 0)return 0;
ll u = m-1, t = 0;
while ((u & 1) == 0) {
u >>= 1;
++t;
}
for (int i = 0; i < 7; ++i) {
ll x = ksm(a[i], u, m);
if (x <= 1 || x == m - 1) continue;
for (int i = 1; i <= t; ++i) {
x = ksm(x, 2, m);
if (x == m - 1 && i != t) {
x = 1;
break;
}
if (x == 1) return 0;
}
if (x != 1) return 0;
}
return 1;
}
void solve() {
while (cin >> n) {
ll m, ans = 0;
for (int i = 1; i <= n; ++i) {
cin >> m;
if (miller_rabin(m)) ans++;
}
cout << ans << '\n';
}
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin>>_t;
while (_t--) {
solve();
}
return 0;
}
17.试除法求质因子
O ( n ) O(\sqrt n) O(n)
int prime[N], cnt = 0;
void getprime(int n) {
for (int i = 2; i <= sqrt(n); ++i) {
if (n % i == 0) {
prime[++cnt] = i;
while (n % i == 0) n /= i;
}
}
if (n > 1) prime[++cnt] = n;
}
18.Pollard_rho 算法
P4718 【模板】Pollard’s rho 算法
Miller Rabin 算法是一种高效的质数判断方法。虽然是一种不确定的质数判断法,但是在选择多种底数的情况下,正确率是可以接受的。
Pollard rho 是一个非常玄学的方式,用于在
O
(
n
1
/
4
)
O(n^{1/4})
O(n1/4) 的期望时间复杂度内计算合数
n
n
n 的某个非平凡因子。事实上算法导论给出的是
O
(
p
)
O(\sqrt p)
O(p),
p
p
p 是
n
n
n 的某个最小因子,满足
p
p
p 与
n
/
p
n/p
n/p 互质。但是这些都是期望,未必符合实际。但事实上 Pollard rho 算法在实际环境中运行的相当不错。
这里我们要写一个程序,对于每个数字检验是否是质数,是质数就输出 Prime
;如果不是质数,输出它最大的质因子是哪个。
O ( n 1 4 ) O(n^{\frac{1}{4}}) O(n41)
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 1e5 + 10;
ll maxn, a[7] = { 2,325,9375,28178,450775,9780504,1795265022 };
inline ll read() {
ll res = 0;
char ch = getchar();
while (ch < '0' || ch>'9')ch = getchar();
while (ch >= '0' && ch <= '9') { res = (res << 3) + (res << 1) + (ch ^ 48); ch = getchar(); }
return res;
}
inline ll gcd(ll a, ll b) {
return b == 0 ? a : gcd(b, a % b);
}
inline ll mul(ll a,ll b,ll MOD){
ll c = a*b-(ll)((long double)a/MOD*b+1e-8)*MOD;
return ((c %= MOD) < 0) ? (c+MOD) : c;
}
inline ll ksm(ll a, ll b, ll p) {
ll res = 1;
a %= p;
while (b) {
if (b & 1) res = mul(res, a, p) % p;
a = mul(a, a, p) % p;
b >>= 1;
}
return res;
}
inline bool mr(ll m) {
if (m == 2) return 1;
if (m <= 1 || m % 2 == 0)return 0;
ll u = m - 1, t = 0;
while ((u & 1) == 0) {
u >>= 1;
++t;
}
for (register int i = 0; i < 7; ++i) {
ll x = ksm(a[i], u, m);
if (x <= 1 || x == m - 1) continue;
for (int i = 1; i <= t; ++i) {
x = ksm(x, 2, m);
if (x == m - 1 && i != t) {
x = 1;
break;
}
if (x == 1) return 0;
}
if (x != 1) return 0;
}
return 1;
}
inline ll pr(ll n) {
if (n == 4) return 2;
ll x = rand() % n, y = x;
ll c = rand() % (n - 3) + 3;
ll d = 1;
x = (mul(x, x, n)+ c)%n;
y = (mul(y, y, n) + c) % n, y = (mul(y, y, n) + c) % n;
for (register int lim = 1; x != y; lim = min(128, lim << 1)) {
ll cnt = 1;
for (register int i = 0; i < lim; ++i) {
ll tmp = mul(cnt, abs(x - y), n);
if (!tmp) break;
cnt = tmp;
x = (mul(x, x, n) + c) % n;
y = (mul(y, y, n) + c) % n, y = (mul(y, y, n) + c) % n;
}
d = gcd(cnt, n);
if (d != 1) return d;
}
return n;
}
inline void findfac(ll x) {
if (x <= maxn || x < 2) return;
if (mr(x)) {
maxn = max(maxn, x);
return;
}
ll p = x;
while (p >= x) p = pr(p);
while (x % p == 0) x /= p;
findfac(x);
findfac(p);
}
inline void solve() {
ll n;
n = read();
maxn = 0;
findfac(n);
if (maxn == n) cout << "Prime\n";
else cout << maxn << '\n';
}
int main() {
//ios::sync_with_stdio(false);
//cin.tie(0);
//cout.tie(0);
int _t;
//_t = 1;
cin >> _t;
while (_t--) {
solve();
}
return 0;
}
19.威尔逊定理
若 p 为素数,则 p 可以整除
(
p
−
1
)
!
+
1
(p-1)!+1
(p−1)!+1
(
p
−
1
)
!
m
o
d
p
=
p
−
1
(p-1)!\ mod\ p=p-1
(p−1)! mod p=p−1
(
p
−
1
)
!
=
p
q
−
1
(p-1)! = pq-1
(p−1)!=pq−1,q为正整数
(
p
−
1
)
!
≡
−
1
(
m
o
d
p
)
(p-1)!\equiv-1(mod\ p)
(p−1)!≡−1(mod p)
20.积性函数
对于任意互素正整数p,q均有
f
(
p
q
)
=
f
(
p
)
f
(
q
)
f(pq)=f(p)f(q)
f(pq)=f(p)f(q),积性函数
对于任意正整数p,q均有
f
(
p
q
)
=
f
(
p
)
f
(
q
)
f(pq)=f(p)f(q)
f(pq)=f(p)f(q),完全积性函数
积性函数的和函数也是积性函数,
F
(
n
)
=
∑
d
∣
n
f
(
d
)
F(n)=\sum\limits_{d|n}f(d)
F(n)=d∣n∑f(d)
单位函数
i
d
(
n
)
=
n
id(n)=n
id(n)=n
幂函数
I
k
(
n
)
=
n
k
I_k(n)=n^k
Ik(n)=nk
元函数
ε
(
n
)
=
{
1
n
=
1
0
n
>
1
\varepsilon(n)=\begin{cases}1& n=1 \\0& n>1 \\\end{cases}
ε(n)={10n=1n>1
因数和函数
σ
(
n
)
=
∑
d
∣
n
d
\sigma(n)=\sum\limits_{d|n}d
σ(n)=d∣n∑d
约数个数
d
(
n
)
=
∑
d
∣
n
1
d(n)=\sum\limits_{d|n}1
d(n)=d∣n∑1
21.欧拉函数
不超过n且与n互素的正整数的个数,
ϕ
(
n
)
=
∑
i
=
1
n
[
g
c
d
(
i
,
n
)
=
1
]
\phi(n)=\sum\limits_{i=1}^n[gcd(i,n)=1]
ϕ(n)=i=1∑n[gcd(i,n)=1]
设p,q是互素正整数,那么
ϕ
(
p
q
)
=
ϕ
(
p
)
ϕ
(
q
)
\phi(pq)=\phi(p)\phi(q)
ϕ(pq)=ϕ(p)ϕ(q)
设
n
=
p
1
a
1
×
p
2
a
2
×
.
.
.
×
p
k
a
k
n=p_1^{a_1}\times p_2^{a_2}\times...\times p_k^{a_k}
n=p1a1×p2a2×...×pkak,则
ϕ
(
n
)
=
ϕ
(
p
1
a
1
)
×
ϕ
(
p
2
a
2
)
×
.
.
.
×
ϕ
(
p
k
a
k
)
\phi(n)=\phi(p_1^{a_1})\times\phi(p_2^{a_2})\times...\times\phi(p_k^{a_k})
ϕ(n)=ϕ(p1a1)×ϕ(p2a2)×...×ϕ(pkak)
设n为正整数,那么
n
=
∑
d
∣
n
ϕ
(
d
)
n=\sum\limits_{d|n}\phi(d)
n=d∣n∑ϕ(d)
欧拉定理:设m是正整数,gcd(a,m)=1,则
a
ϕ
(
m
)
≡
1
(
m
o
d
m
)
a^{\phi(m)}\equiv1(mod\ m)
aϕ(m)≡1(mod m)
设
n
=
p
1
a
1
×
p
2
a
2
×
.
.
.
×
p
k
a
k
n=p_1^{a_1}\times p_2^{a_2}\times...\times p_k^{a_k}
n=p1a1×p2a2×...×pkak为正整数n的质幂因数分解,那么
ϕ
(
n
)
=
n
∏
i
=
1
k
(
1
−
1
p
i
)
\phi(n)=n\prod\limits_{i=1}^k(1-\frac{1}{p_i})
ϕ(n)=ni=1∏k(1−pi1)
(1)若n是素数,
ϕ
(
n
)
=
n
−
1
\phi(n)=n-1
ϕ(n)=n−1
(2)若
n
=
p
k
n=p^k
n=pk,
ϕ
(
n
)
=
p
k
−
p
k
−
1
=
p
k
−
1
(
p
−
1
)
=
p
k
−
1
ϕ
(
p
)
\phi(n)=p^k-p^{k-1}=p^{k-1}(p-1)=p^{k-1}\phi(p)
ϕ(n)=pk−pk−1=pk−1(p−1)=pk−1ϕ(p)
试除法求单个欧拉函数
O ( n ) O(\sqrt n) O(n)
int getphi(int n) {
int res = n;
for (int i = 2; i * i <= n; ++i) {
if (n % i == 0) {
res = res / i * (i - 1);
while (n % i == 0) n /= i;
}
}
if (n != 1) res = res / n * (n - 1);
return res;
}
欧拉筛求欧拉函数
P2158 [SDOI2008] 仪仗队
作为体育委员,C 君负责这次运动会仪仗队的训练。仪仗队是由学生组成的
N
×
N
N \times N
N×N 的方阵,为了保证队伍在行进中整齐划一,C 君会跟在仪仗队的左后方,根据其视线所及的学生人数来判断队伍是否整齐。
现在,C 君希望你告诉他队伍整齐时能看到的学生人数。
O ( n ) O(n) O(n)
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 4e4 + 10;
int prime[N], tag[N], phi[N], sum[N], cnt = 0;
void getphi() {
phi[1] = 1;
for (int i = 2; i < N; ++i) {
if (!tag[i]) {
prime[++cnt] = i;
phi[i] = i - 1;
}
for (int j = 1; j <= cnt && i * prime[j] < N; ++j) {
tag[i * prime[j]] = 1;
if (i % prime[j] == 0) {
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
phi[i * prime[j]] = phi[i] * phi[prime[j]];
}
}
}
void solve() {
getphi();
for (int i = 1; i < N; ++i) sum[i] = sum[i - 1] + phi[i];
int n;
cin >> n;
if (n == 1) cout << 0;
else cout << 2 * sum[n-1] + 1;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin>>_t;
while (_t--) {
solve();
}
return 0;
}
22.扩展欧拉定理
a b = { a b m o d ϕ ( m ) g c d ( a , m ) = 1 a b g c d ( a , m ) ≠ 1 , b < ϕ ( m ) a ( b m o d ϕ ( m ) ) + ϕ ( m ) g c d ( a , m ) ≠ 1 , b ≥ ϕ ( m ) ( m o d m ) a^b=\begin{cases}a^{b\ mod\ \phi(m)}&gcd(a,m)=1\\a^b&gcd(a,m)\neq1,b<\phi(m)\\a^{(b\ mod\ \phi(m))+\phi(m)}&gcd(a,m)\neq1,b\geq\phi(m)\end{cases}(mod\ m) ab=⎩ ⎨ ⎧ab mod ϕ(m)aba(b mod ϕ(m))+ϕ(m)gcd(a,m)=1gcd(a,m)=1,b<ϕ(m)gcd(a,m)=1,b≥ϕ(m)(mod m)
P5091 【模板】扩展欧拉定理
给你三个正整数,
a
,
m
,
b
a,m,b
a,m,b,你需要求:
a
b
m
o
d
m
a^b \bmod m
abmodm
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 1e5 + 10;
ll ksm(ll a, ll b, ll p) {
a %= p;
ll res = 1;
while (b) {
if (b & 1) res = res * a % p;
a = a * a % p;
b >>= 1;
}
return res;
}
ll getphi(ll x) {
ll res = x;
for (ll i = 2; i * i <= x; ++i) {
if (x % i == 0) {
res = res / i * (i - 1);
while (x % i == 0) x /= i;
}
}
if (x > 1) res = res / x * (x - 1);
return res;
}
void solve() {
ll a, b = 0, m;
string s;
bool f = 0;
cin >> a >> m >> s;
ll phi = getphi(m);
int len = s.length();
for (int i = 0; i < len; ++i) {
b = b * 10 + s[i] - '0';
if (b >= phi) {
f = 1;
b %= phi;
}
}
if (f) b += phi;
cout << ksm(a, b, m);
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin>>_t;
while (_t--) {
solve();
}
return 0;
}
23.整除分块
分块数量少于 2 n , O ( n ) 2\sqrt n,O(\sqrt n) 2n,O(n)
#include<bits/stdc++.h>
using namespace std;
#define ll long long
int main() {
ll n, l, r, ans = 0;
cin >> n;
for (l = 1; l <= n; l = r + 1) {
r = n / (n / l);
ans += (r - l + 1) * (n / l);
}
cout << ans;
return 0;
}
24.狄利克雷卷积
(
f
∗
g
)
(
n
)
=
∑
d
∣
n
f
(
d
)
g
(
n
d
)
(f*g)(n)=\sum\limits_{d|n}f(d)g(\frac{n}{d})
(f∗g)(n)=d∣n∑f(d)g(dn)
(
I
∗
1
)
(
n
)
=
σ
(
n
)
(I*1)(n)=\sigma(n)
(I∗1)(n)=σ(n)
交换律,结合律,分配律
两个积性函数的狄利克雷卷积仍是积性函数
25.莫比乌斯反演
莫比乌斯函数
μ
(
n
)
=
{
1
n
=
1
(
−
1
)
r
n
=
p
1
p
2
.
.
.
p
r
,
其中
p
i
为不同素数
0
其他
\mu(n)=\begin{cases}1&n=1\\(-1)^r&n=p_1p_2...p_r,其中p_i为不同素数\\0&其他\end{cases}
μ(n)=⎩
⎨
⎧1(−1)r0n=1n=p1p2...pr,其中pi为不同素数其他
莫比乌斯函数的和函数在整数n处的值
F
(
n
)
=
∑
d
∣
n
μ
(
d
)
F(n)=\sum\limits_{d|n}\mu(d)
F(n)=d∣n∑μ(d)满足
F
(
n
)
=
{
1
n
=
1
0
n
>
1
F(n)=\begin{cases}1&n=1\\0&n>1\end{cases}
F(n)={10n=1n>1
莫比乌斯反演:若 f 是算术函数,F为 f 的和函数,对任意正整数n,满足
F
(
n
)
=
∑
d
∣
n
f
(
d
)
F(n)=\sum\limits_{d|n}f(d)
F(n)=d∣n∑f(d),则有
f
(
n
)
=
∑
d
∣
n
μ
(
d
)
F
(
n
d
)
f(n)=\sum\limits_{d|n}\mu(d)F(\frac{n}{d})
f(n)=d∣n∑μ(d)F(dn)
F
(
n
)
=
∑
n
∣
d
f
(
d
)
F(n)=\sum\limits_{n|d}f(d)
F(n)=n∣d∑f(d),则有
f
(
n
)
=
∑
n
∣
d
μ
(
d
n
)
F
(
d
)
f(n)=\sum\limits_{n|d}\mu(\frac{d}{n})F(d)
f(n)=n∣d∑μ(nd)F(d)
设 f 是算术函数,它的和函数为F(n),如果F是积性函数,则 f 也是积性函数
P2522 [HAOI2011]Problem b
对于给出的
n
n
n 个询问,每次求有多少个数对
(
x
,
y
)
(x,y)
(x,y),满足
a
≤
x
≤
b
a \le x \le b
a≤x≤b,
c
≤
y
≤
d
c \le y \le d
c≤y≤d,且
gcd
(
x
,
y
)
=
k
\gcd(x,y) = k
gcd(x,y)=k,
gcd
(
x
,
y
)
\gcd(x,y)
gcd(x,y) 函数为
x
x
x 和
y
y
y 的最大公约数。
f
(
n
)
=
∑
x
=
1
a
∑
y
=
1
b
g
c
d
(
x
,
y
)
=
n
f(n)=\sum\limits_{x=1}^a\sum\limits_{y=1}^bgcd(x,y)=n
f(n)=x=1∑ay=1∑bgcd(x,y)=n
F
(
n
)
=
∑
x
=
1
a
∑
y
=
1
b
n
∣
g
c
d
(
x
,
y
)
F(n)=\sum\limits_{x=1}^a\sum\limits_{y=1}^bn|gcd(x,y)
F(n)=x=1∑ay=1∑bn∣gcd(x,y)
F
(
d
)
=
⌊
a
d
⌋
⌊
b
d
⌋
F(d)=\lfloor\frac{a}{d}\rfloor\lfloor\frac{b}{d}\rfloor
F(d)=⌊da⌋⌊db⌋
f
(
n
)
=
∑
n
∣
d
μ
(
d
n
)
⌊
a
d
⌋
⌊
b
d
⌋
=
∑
d
′
μ
(
d
′
)
⌊
a
′
d
′
⌋
⌊
b
′
d
′
⌋
f(n)=\sum\limits_{n|d}\mu(\frac{d}{n})\lfloor\frac{a}{d}\rfloor\lfloor\frac{b}{d}\rfloor=\sum\limits_{d^{\prime}}\mu(d^{\prime})\lfloor\frac{a^{\prime}}{d^{\prime}}\rfloor\lfloor\frac{b^{\prime}}{d^{\prime}}\rfloor
f(n)=n∣d∑μ(nd)⌊da⌋⌊db⌋=d′∑μ(d′)⌊d′a′⌋⌊d′b′⌋
d
′
=
d
n
d^{\prime}=\frac{d}{n}
d′=nd
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 1e5 + 10;
int prime[N], tag[N], mu[N], cnt, pre[N];
void init() {
mu[1] = 1;
for (int i = 2; i < N; ++i) {
if (!tag[i]) prime[++cnt] = i, mu[i] = -1;
for (int j = 1; j <= cnt && i * prime[j] < N; ++j) {
tag[i * prime[j]] = 1;
if (i % prime[j] == 0) break;
mu[i * prime[j]] = -mu[i];
}
}
for (int i = 1; i < N; ++i) pre[i] = pre[i - 1] + mu[i];
}
ll f(int x, int y, int k) {
ll res = 0;
x /= k, y /= k;
int n = min(x, y);
for (int l = 1, r; l <= n; l = r + 1) {
r = min(n, min(x / (x / l), y / (y / l)));
res += 1ll * (pre[r] - pre[l - 1]) * (x / l) * (y / l);
}
return res;
}
void solve() {
int a, b, c, d, k;
cin >> a >> b >> c >> d >> k;
cout << f(b, d, k) - f(a - 1, d, k) - f(b, c - 1, k) + f(a - 1, c - 1, k) << '\n';
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
init();
int _t;
//_t = 1;
cin >> _t;
while (_t--) {
solve();
}
return 0;
}
26.二项式反演
形式1
g
n
g_n
gn表示至多n个/n种的方案数量,
f
n
f_n
fn恰好n个/n种的方案数量
g
n
=
∑
i
=
0
n
(
n
i
)
f
i
⟺
f
n
=
∑
i
=
0
n
(
−
1
)
n
−
i
(
n
i
)
g
i
g_n=\sum\limits_{i=0}^n\binom{n}{i}f_i\Longleftrightarrow f_n=\sum\limits_{i=0}^n(-1)^{n-i}\binom{n}{i}g_i
gn=i=0∑n(in)fi⟺fn=i=0∑n(−1)n−i(in)gi
形式2
g
k
g_k
gk表示至少k个/k种的方案数量,
f
k
f_k
fk恰好k个/k种的方案数量
g
k
=
∑
i
=
k
n
(
i
k
)
f
i
⟺
f
k
=
∑
i
=
k
n
(
−
1
)
i
−
k
(
i
k
)
g
i
g_k=\sum\limits_{i=k}^n\binom{i}{k}f_i\Longleftrightarrow f_k=\sum\limits_{i=k}^n(-1)^{i-k}\binom{i}{k}g_i
gk=i=k∑n(ki)fi⟺fk=i=k∑n(−1)i−k(ki)gi
27.杜教筛
记
S
(
n
)
=
∑
i
=
1
n
f
(
i
)
S(n)=\sum\limits_{i=1}^nf(i)
S(n)=i=1∑nf(i),构造两个积性函数h,g满足卷积h=g*f
g
(
1
)
S
(
n
)
=
∑
i
=
1
n
h
(
i
)
−
∑
i
=
2
n
g
(
i
)
S
(
⌊
n
i
⌋
)
g(1)S(n)=\sum\limits_{i=1}^nh(i)-\sum\limits_{i=2}^ng(i)S(\lfloor\frac{n}{i}\rfloor)
g(1)S(n)=i=1∑nh(i)−i=2∑ng(i)S(⌊in⌋)
欧拉函数
h
=
i
d
,
g
=
I
,
S
(
n
)
=
n
(
n
+
1
)
2
−
∑
i
=
2
n
S
(
⌊
n
i
⌋
)
h=id,g=I,S(n)=\frac{n(n+1)}{2}-\sum\limits_{i=2}^nS(\lfloor\frac{n}{i}\rfloor)
h=id,g=I,S(n)=2n(n+1)−i=2∑nS(⌊in⌋)
莫比乌斯函数
h
=
ε
,
g
=
I
,
S
(
n
)
=
1
−
∑
i
=
2
n
S
(
⌊
n
i
⌋
)
h=\varepsilon,g=I,S(n)=1-\sum\limits_{i=2}^nS(\lfloor\frac{n}{i}\rfloor)
h=ε,g=I,S(n)=1−i=2∑nS(⌊in⌋)
P4213 【模板】杜教筛(Sum)
给定一个正整数,求
a
n
s
1
=
∑
i
=
1
n
φ
(
i
)
ans_1=\sum_{i=1}^n\varphi(i)
ans1=i=1∑nφ(i)
a
n
s
2
=
∑
i
=
1
n
μ
(
i
)
ans_2=\sum_{i=1}^n \mu(i)
ans2=i=1∑nμ(i)
O ( n 2 3 ) O(n^{\frac{2}{3}}) O(n32)
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 5e6 + 10;
unordered_map<int, ll> sumphi, summu;
ll prime[N], tag[N], cnt, phi[N], mu[N];
void init() {
phi[1] = mu[1] = 1;
for (int i = 2; i < N; ++i) {
if (!tag[i]) {
prime[++cnt] = i;
phi[i] = i - 1;
mu[i] = -1;
}
for (int j = 1; j <= cnt && i * prime[j] < N; ++j) {
tag[i * prime[j]] = 1;
if (i % prime[j] == 0) {
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
phi[i * prime[j]] = phi[i] * phi[prime[j]];
mu[i * prime[j]] = -mu[i];
}
}
for (int i = 1; i < N; ++i) {
phi[i] += phi[i - 1];
mu[i] += mu[i - 1];
}
}
ll getsumphi(ll x) {
if (x < N) return phi[x];
if (sumphi[x]) return sumphi[x];
ll ans = x * (x + 1) / 2;
for (ll l = 2, r; l <= x; l = r + 1) {
r = x / (x / l);
ans -= (r - l + 1) * getsumphi(x / l);
}
return sumphi[x] = ans;
}
ll getsummu(ll x) {
if (x < N) return mu[x];
if (summu[x]) return summu[x];
ll ans = 1;
for (ll l = 2, r; l <= x; l = r + 1) {
r = x / (x / l);
ans -= (r - l + 1) * getsummu(x / l);
}
return summu[x] = ans;
}
void solve() {
ll n;
cin >> n;
cout << getsumphi(n) << " " << getsummu(n) << '\n';
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
init();
int _t;
//_t = 1;
cin >> _t;
while (_t--) {
solve();
}
return 0;
}
28.BSGS算法
a
i
∗
k
−
j
≡
b
(
m
o
d
p
)
a^{i*k-j}\equiv b(mod\ p)
ai∗k−j≡b(mod p)
a
i
∗
k
≡
b
∗
a
j
(
m
o
d
p
)
a^{i*k}\equiv b*a^j(mod\ p)
ai∗k≡b∗aj(mod p)
P3846 [TJOI2007] 可爱的质数/【模板】BSGS
给定一个质数
p
p
p,以及一个整数
b
b
b,一个整数
n
n
n,现在要求你计算一个最小的非负整数
l
l
l,满足
b
l
≡
n
(
m
o
d
p
)
b^l \equiv n \pmod p
bl≡n(modp)。
O ( p ) O(\sqrt p) O(p)
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 1e5 + 10;
ll ksm(ll a, ll b, ll p) {
ll res = 1;
a %= p;
while (b) {
if (b & 1) res = res * a % p;
a = a * a % p;
b >>= 1;
}
return res;
}
void bsgs(ll a, ll b, ll p) {
unordered_map<int, int> hash;
ll k = sqrt(p) + 1;
for (ll i = 0, j = b; i < k; ++i) {
hash[j] = i;
j = j * a % p;
}
ll t = ksm(a, k, p);
for (ll i = 1, j = t; i <= k; ++i) {
if (hash.count(j)) {
cout << i * k - hash[j];
return;
}
j = j * t % p;
}
cout << "no solution";
}
void solve() {
ll a, b, p;
cin >> p >> a >> b;
bsgs(a, b, p);
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin>>_t;
while (_t--) {
solve();
}
return 0;
}
29.扩展BSGS算法
a x − c n t ∗ a c n t ∏ d ≡ b ∏ d ( m o d p ∏ d ) a^{x-cnt}*\frac{a^{cnt}}{\prod d}\equiv\frac{b}{\prod d}(mod\ \frac{p}{\prod d}) ax−cnt∗∏dacnt≡∏db(mod ∏dp)
P4195 【模板】扩展 BSGS/exBSGS
给定
a
,
p
,
b
a,p,b
a,p,b,求满足
a
x
≡
b
(
m
o
d
p
)
a^x≡b \pmod p
ax≡b(modp) 的最小自然数
x
x
x 。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 1e5 + 10;
ll gcd(ll a, ll b) {
return b == 0 ? a : gcd(b, a % b);
}
ll exgcd(ll a, ll b, ll& x, ll& y) {
if (b == 0) {
x = 1, y = 0;
return a;
}
ll t = exgcd(b, a % b, x, y);
ll x0 = x, y0 = y;
x = y0;
y = x0 - a / b * y0;
return t;
}
ll ny(int a, int p) {
ll x, y;
exgcd(a, p, x, y);
return (x % p + p) % p;
}
ll ksm(ll a, ll b, ll p) {
ll res = 1;
a %= p;
while (b) {
if (b & 1) res = res * a % p;
a = a * a % p;
b >>= 1;
}
return res;
}
ll bsgs(ll a, ll b, ll p) {
unordered_map<int, int> hash;
ll k = ceil(sqrt(p));
for (ll i = 0, j = b; i < k; ++i) {
hash[j] = i;
j = j * a % p;
}
ll t = 1;
for (ll i = 0; i < k; ++i) t = t * a % p;
for (ll i = 1, j = t; i <= k; ++i) {
if (hash.count(j)) {
return i * k - hash[j];
}
j = j * t % p;
}
return -1;
}
ll exbsgs(ll a, ll b, ll p) {
a %= p, b %= p;
if (b == 1 || p == 1) return 0;
ll cnt = 0, d = gcd(a, p),t=1;
while (d != 1) {
if (b % d) return -1;
cnt++;
b /= d, p /= d;
t = t * a / d % p;
if (t == b) return cnt;
d = gcd(a, p);
}
ll ans=bsgs(a, b*ny(t,p)%p, p);
if (ans == -1) return -1;
else return ans + cnt;
}
void solve() {
ll a, b, p;
while (cin >> a >> p >> b && (a || p || b)) {
ll ans=exbsgs(a, b, p);
if (ans == -1) cout << "No Solution\n";
else cout << ans << '\n';
}
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin>>_t;
while (_t--) {
solve();
}
return 0;
}
30.FFT
f
(
ω
n
k
)
=
G
(
ω
n
/
2
k
)
+
ω
n
k
×
H
(
ω
n
/
2
k
)
f(\omega_n^k) = G(\omega_{n/2}^k) + \omega_n^k \times H(\omega_{n/2}^k)
f(ωnk)=G(ωn/2k)+ωnk×H(ωn/2k)
f
(
ω
n
k
+
n
/
2
)
=
G
(
ω
n
/
2
k
)
−
ω
n
k
×
H
(
ω
n
/
2
k
)
f(\omega_n^{k+n/2}) = G(\omega_{n/2}^k) - \omega_n^k \times H(\omega_{n/2}^k)
f(ωnk+n/2)=G(ωn/2k)−ωnk×H(ωn/2k)
O ( n l o g 2 n ) O(nlog_2n) O(nlog2n)
多项式乘法
P3803 【模板】多项式乘法(FFT)
给定一个
n
n
n 次多项式
F
(
x
)
F(x)
F(x),和一个
m
m
m 次多项式
G
(
x
)
G(x)
G(x)。
请求出
F
(
x
)
F(x)
F(x) 和
G
(
x
)
G(x)
G(x) 的卷积。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const double pi = acos(-1);
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 3e6 + 10;
struct complex {
double x, y;
complex operator+(complex a) { return { x + a.x,y + a.y }; }
complex operator-(complex a) { return { x - a.x,y - a.y }; }
complex operator*(complex a) { return { x * a.x - y * a.y,x * a.y + y * a.x }; }
}a[N], b[N];
int rev[N], bit, tot;
void fft(complex* a, int inv) {
for (int i = 0; i < tot; ++i) {
if (i < rev[i]) swap(a[i], a[rev[i]]);
}
for (int mid = 1; mid < tot; mid <<= 1) {
complex w1 = { cos(pi / mid),inv * sin(pi / mid) };
for (int i = 0; i < tot; i += 2 * mid) {
complex wk = { 1,0 };
for (int j = 0; j < mid; ++j, wk = wk * w1) {
complex x = a[i + j], y = wk * a[i + j + mid];
a[i + j] = x + y, a[i + j + mid] = x - y;
}
}
}
}
void solve() {
int n, m;
cin >> n >> m;
for (int i = 0; i <= n; ++i) cin >> a[i].x;
for (int i = 0; i <= m; ++i) cin >> b[i].x;
while ((1 << bit) - 1 < n + m) bit++;
tot = 1 << bit;
for (int i = 0; i < tot; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
fft(a, 1), fft(b, 1);
for (int i = 0; i < tot; ++i) a[i] = a[i] * b[i];
fft(a, -1);
for (int i = 0; i <= n + m; ++i) cout << (int)(a[i].x / tot + 0.5) << " ";
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin>>_t;
while (_t--) {
solve();
}
return 0;
}
高精度乘法
P1919 【模板】A*B Problem 升级版(FFT 快速傅里叶变换)
给你两个正整数
a
,
b
a,b
a,b,求
a
×
b
a \times b
a×b。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const double pi = acos(-1);
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 3e6 + 10;
struct complex {
double x, y;
complex operator+(complex a) { return { x + a.x,y + a.y }; }
complex operator-(complex a) { return { x - a.x,y - a.y }; }
complex operator*(complex a) { return { x * a.x - y * a.y,x * a.y + y * a.x }; }
}a[N], b[N];
int rev[N], ans[N], bit, tot;
void fft(complex* a, int inv) {
for (int i = 0; i < tot; ++i) {
if (i < rev[i]) swap(a[i], a[rev[i]]);
}
for (int mid = 1; mid < tot; mid <<= 1) {
complex w1 = { cos(pi / mid),inv * sin(pi / mid) };
for (int i = 0; i < tot; i += 2 * mid) {
complex wk = { 1,0 };
for (int j = 0; j < mid; ++j, wk = wk * w1) {
complex x = a[i + j], y = wk * a[i + j + mid];
a[i + j] = x + y, a[i + j + mid] = x - y;
}
}
}
}
void solve() {
string s1, s2;
cin >> s1 >> s2;
int n = s1.length() - 1, m = s2.length() - 1;
for (int i = 0; i <= n; ++i) a[i].x = s1[n - i] - '0';
for (int i = 0; i <= m; ++i) b[i].x = s2[m - i] - '0';
while ((1 << bit) - 1 < n + m) bit++;
tot = 1 << bit;
for (int i = 0; i < tot; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
fft(a, 1), fft(b, 1);
for (int i = 0; i < tot; ++i) a[i] = a[i] * b[i];
fft(a, -1);
int cnt = 0, tmp = 0;
for (int i = 0; i <= n + m || tmp; ++i) {
tmp += a[i].x / tot + 0.5;
ans[cnt++] = tmp % 10;
tmp /= 10;
}
for (int i = cnt - 1; i >= 0; --i) cout << ans[i];
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin>>_t;
while (_t--) {
solve();
}
return 0;
}
31.原根
阶:由欧拉定理可知,对 a ∈ Z , m ∈ N ∗ a\in \mathbb{Z},m\in\mathbb{N}^{*} a∈Z,m∈N∗若 gcd ( a , m ) = 1 \gcd(a,m)=1 gcd(a,m)=1则 a φ ( m ) ≡ 1 ( m o d m ) a^{\varphi(m)}\equiv 1\pmod m aφ(m)≡1(modm)。因此满足同余式 a n ≡ 1 ( m o d m ) a^n \equiv 1 \pmod m an≡1(modm)的最小正整数 n 存在,这个 n 称作 a 模 m 的阶,记作 δ m ( a ) \delta_m(a) δm(a)。
性质 1:
a
,
a
2
,
⋯
,
a
δ
m
(
a
)
a,a^2,\cdots,a^{\delta_m(a)}
a,a2,⋯,aδm(a)模 m 两两不同余。
性质 2:若
a
n
≡
1
(
m
o
d
m
)
a^n \equiv 1 \pmod m
an≡1(modm),则
δ
m
(
a
)
∣
n
\delta_m(a)\mid n
δm(a)∣n。
性质 3:设
m
∈
N
∗
,
a
,
b
∈
Z
,
gcd
(
a
,
m
)
=
gcd
(
b
,
m
)
=
1
m\in\mathbb{N}^{*},a,b\in\mathbb{Z},\gcd(a,m)=\gcd(b,m)=1
m∈N∗,a,b∈Z,gcd(a,m)=gcd(b,m)=1,则
δ
m
(
a
b
)
=
δ
m
(
a
)
δ
m
(
b
)
\delta_m(ab)=\delta_m(a)\delta_m(b)
δm(ab)=δm(a)δm(b)的充分必要条件是
gcd
(
δ
m
(
a
)
,
δ
m
(
b
)
)
=
1
\gcd\big(\delta_m(a),\delta_m(b)\big)=1
gcd(δm(a),δm(b))=1
性质 4:设
k
∈
N
,
m
∈
N
∗
,
a
∈
Z
,
gcd
(
a
,
m
)
=
1
k \in \mathbb{N},m\in \mathbb{N}^{*},a\in\mathbb{Z},\gcd(a,m)=1
k∈N,m∈N∗,a∈Z,gcd(a,m)=1,则
δ
m
(
a
k
)
=
δ
m
(
a
)
gcd
(
δ
m
(
a
)
,
k
)
\delta_m(a^k)=\dfrac{\delta_m(a)}{\gcd\big(\delta_m(a),k\big)}
δm(ak)=gcd(δm(a),k)δm(a)
原根:设 m ∈ N ∗ , a ∈ Z m \in \mathbb{N}^{*},a\in \mathbb{Z} m∈N∗,a∈Z。若 g c d ( a , m ) = 1 gcd(a,m)=1 gcd(a,m)=1,且 δ m ( a ) = φ ( m ) \delta_m(a)=\varphi(m) δm(a)=φ(m),则称 a 为模 m 的原根。
1.原根判定定理:设
m
⩾
3
,
gcd
(
a
,
m
)
=
1
m \geqslant 3, \gcd(a,m)=1
m⩾3,gcd(a,m)=1,则 a 是模 m 的原根的充要条件是,对于
φ
(
m
)
\varphi(m)
φ(m)的每个素因数 p,都有
a
φ
(
m
)
p
≢
1
(
m
o
d
m
)
a^{\frac{\varphi(m)}{p}}\not\equiv 1\pmod m
apφ(m)≡1(modm)。
2.原根个数:若一个数 m 有原根,则它原根的个数为
φ
(
φ
(
m
)
)
\varphi(\varphi(m))
φ(φ(m))。
3.原根存在定理:一个数 m 存在原根当且仅当
m
=
2
,
4
,
p
α
,
2
p
α
m=2,4,p^{\alpha},2p^{\alpha}
m=2,4,pα,2pα,其中 p 为奇素数,
α
∈
N
∗
\alpha\in \mathbb{N}^{*}
α∈N∗。
4.若 m 有原根 g,则若
gcd
(
k
,
φ
(
m
)
)
=
1
\gcd\big(k,\varphi(m)\big)=1
gcd(k,φ(m))=1,则
g
k
g^k
gk 也是模 m 的原根
P6091 【模板】原根
给定整数
n
n
n,求它的所有原根。
为了减小你的输出量,给出输出参数
d
d
d,设
n
n
n 的所有原根有
c
c
c 个,从小到大分别为
g
1
,
…
,
g
c
g_1,\ldots,g_c
g1,…,gc,你只需要依次输出
g
d
,
g
2
d
,
…
,
g
⌊
c
d
⌋
×
d
g_d,g_{2d},\ldots,g_{\lfloor\frac{c}{d}\rfloor\times d}
gd,g2d,…,g⌊dc⌋×d。
如果你不了解原根的定义,可以自行查找资料或阅读下列定义:
正整数
g
g
g 是正整数
n
n
n 的原根,当且仅当
1
≤
g
≤
n
−
1
1\leq g\leq n-1
1≤g≤n−1,且
g
g
g 模
n
n
n 的阶为
φ
(
n
)
\varphi(n)
φ(n)。
O ( n l o g n ) O(nlogn) O(nlogn)
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 1e6 + 10;
int phi[N], rt[N], prime[N], tag[N], cnt = 0, num = 0, fac[N], ans[N];
void init() {
phi[1] = 1;
for (int i = 2; i < N; ++i) {
if (!tag[i]) {
prime[++cnt] = i;
phi[i] = i - 1;
}
for (int j = 1; j <= cnt && i * prime[j] < N; ++j) {
tag[i * prime[j]] = 1;
if (i % prime[j] == 0) {
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
phi[i * prime[j]] = phi[i] * phi[prime[j]];
}
}
rt[2] = rt[4] = 1;
for (int i = 2; i <= cnt; ++i) {
for (ll j = prime[i]; j < N; j *= prime[i]) {
rt[j] = 1;
if (2 * j < N) rt[2 * j] = 1;
}
}
}
ll ksm(ll a, ll b, ll p) {
ll res = 1;
while (b) {
if (b & 1) res = res * a % p;
a = a * a % p;
b >>= 1;
}
return res;
}
int gcd(int a, int b) {
return b == 0 ? a : gcd(b, a % b);
}
void getfactor(int x) {
for (int i = 2; i * i <= x; ++i) {
if (x % i == 0) {
fac[++cnt] = i;
while (x % i == 0) x /= i;
}
}
if (x > 1) fac[++cnt] = x;
}
bool check(int p, int x) {
if (ksm(p, phi[x], x) != 1) return 0;
for (int i = 1; i <= cnt; ++i) {
if (ksm(p, phi[x] / fac[i], x) == 1) return 0;
}
return 1;
}
int getmin(int x) {
for (int i = 1; i < x; ++i) {
if (check(i, x)) return i;
}
}
void getrt(int x, int mn) {
int cur = 1;
for (int i = 1; i <= phi[x]; ++i) {
cur = 1ll * cur * mn % x;
if (gcd(i, phi[x]) == 1) ans[++num] = cur;
}
}
void solve() {
int n, d;
cin >> n >> d;
if (!rt[n]) {
cout << 0 << "\n\n";
return;
}
num = cnt = 0;
getfactor(phi[n]);
int mn = getmin(n);
getrt(n, mn);
sort(ans + 1, ans + num + 1);
cout << num << '\n';
for (int i = 1; i <= num / d; ++i) cout << ans[i * d] << " ";
cout << '\n';
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
init();
int _t;
//_t = 1;
cin >> _t;
while (_t--) {
solve();
}
return 0;
}
32.NTT
多项式乘法
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 998244353;
const ll INF = 1e18;
const ll N = 3e6 + 10;
const int g = 3;
ll a[N], b[N], rev[N], bit, tot;
ll ksm(ll a, ll b) {
ll res = 1;
while (b) {
if (b & 1) res = res * a % mod;
a = a * a % mod;
b >>= 1;
}
return res;
}
void ntt(ll a[], int inv) {
for (int i = 0; i < tot; ++i) {
if (i < rev[i]) swap(a[i], a[rev[i]]);
}
for (int mid = 1; mid < tot; mid <<= 1) {
ll w1 = ksm(g, (mod - 1) / (mid * 2));
if (inv == -1) w1 = ksm(w1, mod - 2);
for (int i = 0; i < tot; i += 2 * mid) {
ll wk = 1;
for (int j = 0; j < mid; ++j, wk = wk * w1 % mod) {
int x = a[i + j], y = wk * a[i + j + mid] % mod;
a[i + j] = (x + y) % mod, a[i + j + mid] = (x - y + mod) % mod;
}
}
}
}
void solve() {
int n, m;
cin >> n >> m;
for (int i = 0; i <= n; ++i) cin >> a[i];
for (int i = 0; i <= m; ++i) cin >> b[i];
while ((1 << bit) - 1 < n + m) bit++;
tot = 1 << bit;
for (int i = 0; i < tot; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
ntt(a, 1), ntt(b, 1);
for (int i = 0; i < tot; ++i) a[i] = a[i] * b[i] % mod;
ntt(a, -1);
ll ny = ksm(tot, mod - 2);
for (int i = 0; i <= n + m; ++i) cout << a[i] * ny % mod << " ";
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin >> _t;
while (_t--) {
solve();
}
return 0;
}
多项式逆元
P4238 【模板】多项式乘法逆
给定一个多项式
F
(
x
)
F(x)
F(x) ,请求出一个多项式
G
(
x
)
G(x)
G(x), 满足
F
(
x
)
∗
G
(
x
)
≡
1
(
m
o
d
x
n
)
F(x) * G(x) \equiv 1 \pmod{x^n}
F(x)∗G(x)≡1(modxn)。系数对
998244353
998244353
998244353 取模。
B ≡ 2 B ′ − A B ′ 2 ( m o d x z ) B\equiv2B^{\prime}−AB^{\prime2}(mod\ x^z) B≡2B′−AB′2(mod xz)
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 998244353;
const ll INF = 1e18;
const ll N = 3e6 + 10;
const int g = 3;
ll a[N], b[N], c[N], rev[N], bit, tot;
ll ksm(ll a, ll b) {
ll res = 1;
while (b) {
if (b & 1) res = res * a % mod;
a = a * a % mod;
b >>= 1;
}
return res;
}
void ntt(ll a[], int inv) {
for (int i = 0; i < tot; ++i) {
if (i < rev[i]) swap(a[i], a[rev[i]]);
}
for (int mid = 1; mid < tot; mid <<= 1) {
ll w1 = ksm(g, (mod - 1) / (mid * 2));
if (inv == -1) w1 = ksm(w1, mod - 2);
for (int i = 0; i < tot; i += 2 * mid) {
ll wk = 1;
for (int j = 0; j < mid; ++j, wk = wk * w1 % mod) {
int x = a[i + j], y = wk * a[i + j + mid] % mod;
a[i + j] = (x + y) % mod, a[i + j + mid] = (x - y + mod) % mod;
}
}
}
}
void work(int n, ll a[], ll b[]) {
if (n == 1) {
b[0] = ksm(a[0], mod - 2);
return;
}
work((n + 1) >> 1, a, b);
bit = 0;
while ((1 << bit) < (n<<1)) bit++;
tot = 1 << bit;
for (int i = 1; i < tot; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
for (int i = 0; i < n; ++i) c[i] = a[i];
for (int i = n; i < tot; ++i) c[i] = 0;
ntt(c, 1), ntt(b, 1);
for (int i = 0; i < tot; ++i) b[i] = (2 - c[i] * b[i] % mod + mod) % mod * b[i] % mod;
ntt(b, -1);
int ny = ksm(tot, mod - 2);
for (int i = 0; i < tot; ++i) b[i] = b[i] * ny % mod;
for (int i = n; i < tot; ++i) b[i] = 0;
}
void solve() {
int n;
cin >> n;
for (int i = 0; i < n; ++i) cin >> a[i];
work(n, a, b);
for (int i = 0; i < n ; ++i) cout << b[i] << " ";
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin >> _t;
while (_t--) {
solve();
}
return 0;
}
33.多项式幂函数
O ( n l o g n ) O(nlogn) O(nlogn)
常数项为1
P5245 【模板】多项式快速幂
给定一个
n
−
1
n-1
n−1 次多项式
A
(
x
)
A(x)
A(x),求一个在
m
o
d
x
n
\bmod\ x^n
mod xn 意义下的多项式
B
(
x
)
B(x)
B(x),使得
B
(
x
)
≡
(
A
(
x
)
)
k
(
m
o
d
x
n
)
B(x) \equiv (A(x))^k \ (\bmod\ x^n)
B(x)≡(A(x))k (mod xn)。保证
a
0
=
1
a_0=1
a0=1。
多项式的系数在
m
o
d
998244353
\bmod\ 998244353
mod 998244353 的意义下进行运算。
#include<bits/stdc++.h>
#define ll long long
const int mod = 998244353;
const int N = 262145;
const int g3 = (mod + 1) / 3;
int n, k, a[N], A[N], G[N];
inline void upd(int& a) { a += a >> 31 & mod; }
inline int pow(int a, int b) {
int r = 1;
for (; b; b >>= 1, a = (ll)a * a % mod)if (b & 1)r = (ll)r * a % mod; return r;
}
namespace poly {
int lim, rev[N], inv[N];
inline int poly_start() {
inv[1] = 1;
for (int i = 2; i < N; ++i)inv[i] = (mod - mod / i) * (ll)inv[mod % i] % mod;
return 1;
}
int __START__ = poly_start();
inline void init(int n) {
int l = -1;
for (lim = 1; lim < n; lim <<= 1)++l;
for (int i = 1; i < lim; ++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << l);
}
void NTT(int* a, int f) {
for (int i = 1; i < lim; ++i)if (i < rev[i])std::swap(a[i], a[rev[i]]);
for (int i = 1; i < lim; i <<= 1) {
const int gi = pow(f ? 3 : g3, (mod - 1) / (i << 1));
for (int j = 0; j < lim; j += i << 1)
for (int k = 0, g = 1; k < i; ++k, g = (ll)g * gi % mod) {
const int x = a[j + k], y = a[j + k + i] * (ll)g % mod;
upd(a[j + k] += y - mod), upd(a[j + k + i] = x - y);
}
}
if (!f) {
const ll iv = inv[lim];
for (int i = 0; i < lim; ++i)a[i] = a[i] * iv % mod;
}
}
void INV(const int* a, int* B, int n) {
if (n == 1)*B = pow(*a, mod - 2); else {
INV(a, B, n + 1 >> 1);
init(n << 1);
static int A[N];
for (int i = 0; i < n; ++i)A[i] = a[i];
for (int i = n; i < lim; ++i)A[i] = 0;
NTT(A, 1), NTT(B, 1);
for (int i = 0; i < lim; ++i)B[i] = B[i] * ((2 - (ll)A[i] * B[i] % mod + mod) % mod) % mod;
NTT(B, 0);
for (int i = n; i < lim; ++i)B[i] = 0;
}
}
void dao(const int* a, int* G, int n) {
G[n - 1] = 0;
for (int i = 1; i < n; ++i)G[i - 1] = (ll)i * a[i] % mod;
}
void INT(const int* a, int* G, int n) {
*G = 0;
for (int i = n - 1; ~i; --i)G[i + 1] = a[i] * (ll)inv[i + 1] % mod;
}
void LN(const int* a, int* B, int n) {
static int F[N];
dao(a, F, n), INV(a, B, n);
init(n << 1);
for (int i = n; i < lim; ++i)B[i] = F[i] = 0;
NTT(F, 1), NTT(B, 1);
for (int i = 0; i < lim; ++i)F[i] = (ll)B[i] * F[i] % mod;
NTT(F, 0);
INT(F, B, n);
for (int i = n; i < lim; ++i)B[i] = 0;
}
void EXP(const int* a, int* F, int n) {
if (n == 1)*F = 1; else {
EXP(a, F, n + 1 >> 1);
static int F0[N], A[N];
for (int i = 0; i <= n << 1; ++i)F0[i] = 0, A[i] = a[i];
LN(F, F0, n);
init(n << 1);
for (int i = n; i < lim; ++i)A[i] = 0;
NTT(A, 1), NTT(F0, 1), NTT(F, 1);
for (int i = 0; i < lim; ++i)F[i] = F[i] * (A[i] + 1ll - F0[i] + mod) % mod;
NTT(F, 0);
for (int i = n; i < lim; ++i)F[i] = 0;
}
}
}
int main() {
scanf("%d", &n);
int c = getchar(); while (isspace(c))c = getchar();
for (; isdigit(c); c = getchar())k = (k * 10ll + (c ^ '0')) % mod;
for (int i = 0; i < n; ++i)scanf("%d", a + i);
poly::LN(a, A, n);
for (int i = 0; i < n; ++i)A[i] = A[i] * (ll)k % mod;
for (int i = 0; i < n; ++i)a[i] = 0;
poly::EXP(A, a, n);
for (int i = 0; i < n; ++i)printf("%d ", a[i]);
return 0;
}
任意
P5273 【模板】多项式幂函数(加强版)
给定一个
n
−
1
n-1
n−1 次多项式
A
(
x
)
A(x)
A(x),求一个在
m
o
d
x
n
\bmod\ x^n
mod xn 意义下的多项式
B
(
x
)
B(x)
B(x),使得
B
(
x
)
≡
(
A
(
x
)
)
k
(
m
o
d
x
n
)
B(x) \equiv (A(x))^k \ (\bmod\ x^n)
B(x)≡(A(x))k (mod xn)。
多项式的系数在
m
o
d
998244353
\bmod\ 998244353
mod 998244353 的意义下进行运算。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 4e5 + 10;
const int mod = 998244353;
const int G = 3;
const int invG = 332748118;
int n, m, m1, m2;
int inv[N], f[N], g[N], h[N], A[N], B[N], C[N], D[N], pos[N];
inline int power(int x, int k){
int res = 1;
while (k){
if (k & 1)res = 1ll * res * x % mod;
x = 1ll * x * x % mod; k >>= 1;
}
return res;
}
inline void NTT(int* a, int op, int len){
int lim = 1, cnt = 0;
while (lim < len)lim <<= 1, cnt++;
for (int i = 0; i < lim; i++)pos[i] = (pos[i >> 1] >> 1) | ((i & 1) << (cnt - 1));
for (int i = 0; i < lim; i++)if (i < pos[i])swap(a[i], a[pos[i]]);
for (int l = 1; l < lim; l <<= 1){
int wn = power(op == 1 ? G : invG, (mod - 1) / (l << 1));
for (int i = 0; i < lim; i += l << 1){
int w = 1;
for (int j = 0; j < l; j++, w = 1ll * w * wn % mod){
int x = a[i + j], y = 1ll * w * a[i + l + j] % mod;
a[i + j] = (x + y) % mod; a[i + l + j] = (x - y + mod) % mod;
}
}
}
if (op == 1)return;
int inv = power(lim, mod - 2);
for (int i = 0; i < lim; i++)a[i] = 1ll * a[i] * inv % mod;
}
inline void getinv(int* a, int* b, int len){
if (len == 1) { b[0] = power(a[0], mod - 2); return; }
getinv(a, b, len >> 1);
for (int i = 0; i < len; i++)D[i] = a[i];
NTT(b, 1, len << 1); NTT(D, 1, len << 1);
for (int i = 0; i < (len << 1); i++)b[i] = 1ll * b[i] * (2ll - 1ll * b[i] * D[i] % mod + mod) % mod;
NTT(b, -1, len << 1);
for (int i = 0; i < len; i++)D[i] = 0;
for (int i = len; i < (len << 1); i++)b[i] = D[i] = 0;
}
inline void Direv(int* a, int* b, int len){
for (int i = 1; i < len; i++)b[i - 1] = 1ll * i * a[i] % mod;
b[len - 1] = 0;
}
inline void Inter(int* a, int* b, int len){
for (int i = 1; i < len; i++)b[i] = 1ll * a[i - 1] * inv[i] % mod;
b[0] = 0;
}
void ln(int* a, int* b, int len){
Direv(a, B, len); getinv(a, C, len);
NTT(B, 1, len << 1); NTT(C, 1, len << 1);
for (int i = 0; i < (len << 1); i++)B[i] = 1ll * B[i] * C[i] % mod;
NTT(B, -1, len << 1); Inter(B, b, len << 1);
for (int i = 0; i < (len << 1); i++)B[i] = C[i] = 0;
}
void exp(int* a, int* b, int len){
if (len == 1) { b[0] = 1; return; }
exp(a, b, len >> 1); ln(b, A, len);
A[0] = (a[0] + 1 - A[0] + mod) % mod;
for (int i = 1; i < len; i++)A[i] = (a[i] - A[i] + mod) % mod;
NTT(b, 1, len << 1); NTT(A, 1, len << 1);
for (int i = 0; i < (len << 1); i++)b[i] = 1ll * b[i] * A[i] % mod;
NTT(b, -1, len << 1);
for (int i = len; i < (len << 1); i++)b[i] = A[i] = 0;
}
int main(){
scanf("%d", &n);
char c = getchar(); c = getchar();
while (c >= '0' && c <= '9'){
m = (1ll * m * 10 % mod + (c - '0')) % mod;
m1 = (1ll * m1 * 10 % (mod - 1) + (c - '0')) % (mod - 1);
if ((ll)(m2 * 10 + c - '0') <= mod)m2 = m2 * 10 + c - '0';
c = getchar();
}
for (int i = 0; i < n; i++)scanf("%d", &f[i]);
int num = n;
for (int i = 0; i < n; i++)if (f[i]) { num = i; break; }
if (1ll * num * m2 >= 1ll * n){
for (int i = 0; i < n; i++)printf("%d ", 0);
return 0;
}
int tmpn = n; n -= num;
for (int i = 0; i < n; i++)f[i] = f[i + num];
int f0 = f[0], invf0 = power(f0, mod - 2);
for (int i = 0; i < n; i++)f[i] = 1ll * f[i] * invf0 % mod;
int lim = 1;
while (lim < n)lim <<= 1;
inv[1] = 1; for (int i = 2; i < lim; i++)inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
ln(f, g, lim);
for (int i = 0; i < n; i++)g[i] = 1ll * g[i] * m % mod;
for (int i = n; i < lim; i++)g[i] = 0;
exp(g, h, lim);
f0 = power(f0, m1);
for (int i = 0; i < n; i++)h[i] = 1ll * h[i] * f0 % mod;
n = tmpn; num = min(1ll * num * m2, 1ll * n);
for (int i = n - 1; i >= num; i--)h[i] = h[i - num];
for (int i = num - 1; ~i; i--)h[i] = 0;
for (int i = 0; i < n; i++)printf("%d ", h[i]);
return 0;
}
34.拉格朗日插值
n
+
1
n+1
n+1个
x
x
x坐标不同的点可以确定唯一的最高为
n
n
n次的多项式
f
(
k
)
=
∑
i
=
0
n
y
i
∏
i
≠
j
k
−
x
[
j
]
x
[
i
]
−
x
[
j
]
f(k) = \sum_{i = 0}^{n} y_i \prod\limits_{i \not = j} \frac{k - x[j]}{x[i] - x[j]}
f(k)=∑i=0nyii=j∏x[i]−x[j]k−x[j]
O
(
n
2
)
O(n^2)
O(n2)
x
x
x取值连续:
f
(
k
)
=
∑
i
=
0
n
y
i
p
r
e
i
−
1
∗
s
u
f
i
+
1
i
!
∗
(
n
−
i
)
!
f(k) = \sum\limits_{i=0}^n y_i \frac{pre_{i-1} * suf_{i+1}}{i! * (n - i)!}
f(k)=i=0∑nyii!∗(n−i)!prei−1∗sufi+1
p
r
e
i
=
∏
j
=
0
i
k
−
j
pre_i = \prod\limits_{j = 0}^{i} k - j
prei=j=0∏ik−j
s
u
f
i
=
∏
j
=
i
n
k
−
j
suf_i = \prod\limits_{j = i}^n k - j
sufi=j=i∏nk−j
当
n
−
i
n - i
n−i为奇数时,分母取负号
O
(
n
)
O(n)
O(n)
重心拉格朗日插值法:
设
g
=
∏
i
=
0
n
k
−
x
[
i
]
,
t
i
=
y
i
∏
j
≠
i
x
i
−
x
j
g = \prod\limits_{i=0}^n k - x[i]\ ,\ t_i = \frac{y_i}{\prod\limits_{j \not =i} x_i - x_j}
g=i=0∏nk−x[i] , ti=j=i∏xi−xjyi
f
(
k
)
=
g
∑
i
=
0
n
t
i
(
k
−
x
[
i
]
)
f(k) = g\sum\limits_{i = 0}^{n} \frac{t_i}{(k - x[i])}
f(k)=gi=0∑n(k−x[i])ti
这样每次新加入一个点的时候只需要计算它的
t
i
t_i
ti
应用:计算 ∑ i = 1 n i k \sum\limits_{i=1}^n i^k i=1∑nik
P4781 【模板】拉格朗日插值
n
n
n 个点
(
x
i
,
y
i
)
(x_i,y_i)
(xi,yi) 可以唯一地确定一个多项式
y
=
f
(
x
)
y = f(x)
y=f(x)。
现在,给定这
n
n
n 个点,请你确定这个多项式,并求出
f
(
k
)
m
o
d
998244353
f(k) \bmod 998244353
f(k)mod998244353 的值。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<int,int>
const double eps = 1e-8;
const ll mod = 998244353;
const ll INF = 1e18;
const ll N = 1e5 + 10;
ll x[N], y[N];
inline ll ksm(ll a, ll b, ll p) {
ll res = 1;
while (b) {
if (b & 1) res = res * a % p;
a = a * a % p;
b >>= 1;
}
return res;
}
inline ll ny(ll x) {
return ksm(x, mod - 2, mod);
}
void solve() {
ll n, k;
cin >> n >> k;
for (register int i = 1; i <= n; ++i) cin >> x[i] >> y[i];
ll ans = 0;
for (register int i = 1; i <= n; ++i) {
ll res = 1;
for (register int j = 1; j <= n; ++j) {
if (i == j) continue;
res = res * (k - x[j]) % mod * ny(x[i] - x[j]) % mod;
}
ans = (ans + res * y[i] % mod) % mod;
}
cout << (ans + mod) % mod;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t;
_t = 1;
//cin>>_t;
while (_t--) {
solve();
}
return 0;
}
35.Min_25筛
P5325 【模板】Min_25筛
定义积性函数
f
(
x
)
f(x)
f(x),且
f
(
p
k
)
=
p
k
(
p
k
−
1
)
f(p^k)=p^k(p^k-1)
f(pk)=pk(pk−1)(
p
p
p是一个质数),求
∑
i
=
1
n
f
(
i
)
\sum_{i=1}^n f(i)
i=1∑nf(i)
对
1
0
9
+
7
10^9+7
109+7取模。
36.行列式求值
P7112 【模板】行列式求值
给定一个
n
n
n 阶行列式
A
A
A,求
∣
A
∣
|A|
∣A∣。结果对
p
p
p 取模。
O
(
n
3
)
O(n^3)
O(n3)
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define P pair<ll,ll>
const double eps = 1e-8;
const ll mod = 1e9 + 7;
const ll INF = 1e18;
const ll N = 1e3 + 10;
ll a[N][N];
inline void solve() {
ll n, p, f = 1;
cin >> n >> p;
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= n; ++j) cin >> a[i][j];
}
for (int i = 1; i <= n; ++i) {
for (int j = i + 1; j <= n; ++j) {
while (a[i][i]) {
ll t = a[j][i] / a[i][i];
for (int k = i; k <= n && t; ++k) a[j][k] = (a[j][k] - t * a[i][k] % p + p) % p;
swap(a[i], a[j]);
f = -f;
}
swap(a[i], a[j]);
f = -f;
}
}
ll ans = 1;
for (int i = 1; i <= n; ++i) ans = (ans * a[i][i]) % p;
cout << (ans * f % p + p) % p;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int _t = 1;
//_t = read();
while (_t--) {
solve();
}
return 0;
}