-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathalgorithm.tex
More file actions
2997 lines (2460 loc) · 105 KB
/
algorithm.tex
File metadata and controls
2997 lines (2460 loc) · 105 KB
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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
\documentclass[UTF8]{article}
\usepackage[UTF8]{ctex}
\usepackage[hidelinks]{hyperref}
\usepackage{listings}
\usepackage{xcolor}
\usepackage{needspace} % 新增:用于控制分页的宏包
\usepackage{enumitem}
\usepackage{amsmath, amssymb}
\usepackage[a4paper, margin=1in]{geometry}
% ===== 代码块样式设置 =====
\lstset{
language=C++,
basicstyle=\ttfamily\small,
keywordstyle=\color{blue},
commentstyle=\color{gray},
numbers=left,
numberstyle=\tiny,
frame=tb,
breaklines=true,
tabsize=4
}
\title{ACM 算法模板}
\author{Otter.put}
\date{2025.05.26}
\begin{document}
\maketitle
% ===== 自动生成目录 =====
\tableofcontents
% ===== 强制分页示例 =====
\section{基础算法}
\subsection{高精度}
\subsubsection{高精度加法}
C = A + B,满足 A $\ge$ 0, B $\ge$ 0;A、B、C 均为逆序 vector(即高位在前,低位在后)
\begin{lstlisting}[caption=高精度加法]
vector<int> add(vector<int> &A, vector<int> &B) {
if (A.size() < B.size()) return add(B, A);
vector<int> C;
int t = 0;
for (int i = 0; i < A.size(); i ++ )
{
t += A[i];
if (i < B.size()) t += B[i];
C.push_back(t % 10);
t /= 10;
}
if (t) C.push_back(t);
return C;
}
\end{lstlisting}
\subsubsection{高精度减法}
C = A - B,满足 A $\ge$ B, A $\ge$ 0, B $\ge$ 0;A、B、C 均为逆序 vector(即高位在前,低位在后)
\begin{lstlisting}[caption=高精度减法]
vector<int> sub(vector<int> &A, vector<int> &B) {
vector<int> C;
for (int i = 0, t = 0; i < A.size(); i ++ ) {
t = A[i] - t;
if (i < B.size()) t -= B[i];
C.push_back((t + 10) % 10);
if (t < 0) t = 1;
else t = 0;
}
while (C.size() > 1 && C.back() == 0) C.pop_back();
return C;
}
\end{lstlisting}
\subsubsection{高精度乘低精度}
C = A * B,满足 A $\ge$ 0, B $\textgreater$ 0;A、C 均为逆序 vector(即高位在前,低位在后)
\begin{lstlisting}[caption=高精度乘低精度]
vector<int> mul(vector<int> &A, int b) {
vector<int> C;
int t = 0;
for (int i = 0; i < A.size() || t; i ++ ) {
if (i < A.size()) t += A[i] * b;
C.push_back(t % 10);
t /= 10;
}
return C;
}
\end{lstlisting}
\subsubsection{高精度乘高精度}
C = A * B,满足 A $\ge$ 0, B $\textgreater$ 0;A、B、C 均为逆序 vector(即高位在前,低位在后)
\begin{lstlisting}[caption=高精度乘高精度][escapeinside=``]
bool is_zero(vector<int>& A) {
return A.size() == 1 && A[0] == 0;
}
vector<int> mul(vector<int>& A, vector<int>& B) {
if(is_zero(A) || is_zero(B))
return {0};
int m = A.size(), n = B.size();
vector<int> C(m + n, 0);
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
C[i + j] += A[i] * B[j]; // 累加到对应位置
// 统一处理进位
int carry = 0;
for (int i = 0; i < m + n; i++) {
int total = C[i] + carry;
C[i] = total % 10;
carry = total / 10;
}
if (carry)
C.push_back(carry);
while (C.size() > 1 && C.back() == 0)
C.pop_back();
return C;
}
\end{lstlisting}
\subsubsection{高精度除低精度}
A / b = C ... r,A $\ge$ 0, b $\textgreater$ 0;A、C 均为逆序 vector(即高位在前,低位在后)
\begin{lstlisting}[caption=高精度除低精度]
vector<int> div(vector<int> &A, int b, int &r) {
vector<int> C;
r = 0;
for (int i = A.size() - 1; i >= 0; i -- ) {
r = r * 10 + A[i];
C.push_back(r / b);
r %= b;
}
reverse(C.begin(), C.end());
while (C.size() > 1 && C.back() == 0) C.pop_back();
return C;
}
\end{lstlisting}
\subsubsection{封装为类 + 压位高精度}
注意,我们封装的高精度类为了实现的便捷没有符号的概念。同时,使用压位高精度的思想来提升运算速度。
\begin{lstlisting}[caption=压位高精度类]
const int BASE = 10000;
const int BASE_DIGITS = 4;
struct Big {
vector<int> digits;
// 当前数 > other,返回 1;小于返回 -1;等于返回 0
int cmp(const Big& other) const {
if (digits.size() != other.digits.size())
return digits.size() > other.digits.size() ? 1 : -1;
for (int i = digits.size() - 1; i >= 0; i--)
if (digits[i] != other.digits[i])
return digits[i] > other.digits[i] ? 1 : -1;
return 0;
}
Big() : digits({0}) {}
Big(long long num) {
if (num == 0) {
digits = {0};
return;
}
while (num > 0) {
digits.push_back(num % BASE);
num /= BASE;
}
}
Big(const string& str) {
for (int i = str.size() - 1; i >= 0; i -= BASE_DIGITS) {
int digit = 0;
int end = max(0, i - BASE_DIGITS + 1);
for (int j = end; j <= i; j++) {
digit = digit * 10 + (str[j] - '0');
}
digits.push_back(digit);
}
while (digits.size() > 1 && digits.back() == 0)
digits.pop_back();
}
Big(const Big& other) : digits(other.digits) {}
Big& operator=(const Big& other) {
digits = other.digits;
return *this;
}
Big operator+(const Big& other) const {
Big result;
result.digits.clear();
int carry = 0;
int maxSize = max(digits.size(), other.digits.size());
for (int i = 0; i < maxSize || carry; i++) {
int sum = carry;
if (i < digits.size()) sum += digits[i];
if (i < other.digits.size()) sum += other.digits[i];
carry = sum >= BASE;
if (carry) sum -= BASE;
result.digits.push_back(sum);
}
return result;
}
// 假设当前数 >= other
Big operator-(const Big& other) const {
Big result;
result.digits.clear();
int borrow = 0;
for (int i = 0; i < digits.size(); i++) {
int diff = digits[i] - borrow;
if (i < other.digits.size()) diff -= other.digits[i];
borrow = 0;
if (diff < 0) {
diff += BASE;
borrow = 1;
}
result.digits.push_back(diff);
}
// 去除前导零
while (result.digits.size() > 1 && result.digits.back() == 0)
result.digits.pop_back();
return result;
}
Big operator*(const int& other) const {
if (other == 0) return Big(0);
Big result;
result.digits.clear();
long long carry = 0;
for (int i = 0; i < digits.size() || carry; i++) {
if (i < digits.size()) carry += (long long)digits[i] * other;
result.digits.push_back(carry % BASE);
carry /= BASE;
}
while (result.digits.size() > 1 && result.digits.back() == 0)
result.digits.pop_back();
return result;
}
Big operator*(const Big& other) const {
if (other == 0) return Big(0);
Big result;
result.digits.resize(digits.size() + other.digits.size(), 0);
for (int i = 0; i < digits.size(); i++) {
long long carry = 0;
for (int j = 0; j < other.digits.size() || carry; j++) {
long long product = result.digits[i + j] + carry;
if (j < other.digits.size())
product += (long long)digits[i] * other.digits[j];
result.digits[i + j] = product % BASE;
carry = product / BASE;
}
}
// 去除前导零
while (result.digits.size() > 1 && result.digits.back() == 0)
result.digits.pop_back();
return result;
}
Big operator/(const int& divisor) const {
Big result;
result.digits.resize(digits.size());
long long remainder = 0;
for (int i = digits.size() - 1; i >= 0; i--) {
long long current = remainder * BASE + digits[i];
result.digits[i] = current / divisor;
remainder = current % divisor;
}
// 去除前导零
while (result.digits.size() > 1 && result.digits.back() == 0)
result.digits.pop_back();
return result;
}
int operator%(const int& divisor) const {
long long remain = 0;
for (int i = digits.size() - 1; i >= 0; i--) {
remain = (remain * BASE + digits[i]) % divisor;
}
return remain;
}
Big operator^(long long exponent) const {
Big base(*this);
Big result(1);
while (exponent > 0) {
if (exponent & 1)
result = result * base;
base = base * base;
exponent >>= 1;
}
return result;
}
bool operator<(const Big& other) const { return cmp(other) == -1; }
bool operator<=(const Big& other) const { return cmp(other) <= 0; }
bool operator>(const Big& other) const { return cmp(other) == 1; }
bool operator>=(const Big& other) const { return cmp(other) >= 0; }
bool operator==(const Big& other) const { return cmp(other) == 0; }
bool operator!=(const Big& other) const { return cmp(other) != 0; }
bool operator<(long long num) const { return *this < Big(num); }
bool operator>(long long num) const { return *this > Big(num); }
bool operator<=(long long num) const { return *this <= Big(num); }
bool operator>=(long long num) const { return *this >= Big(num); }
bool operator==(long long num) const { return *this == Big(num); }
bool operator!=(long long num) const { return *this != Big(num); }
void print() const {
cout << digits.back(); // 最高位不需要前导零
for (int i = digits.size() - 2; i >= 0; i--)
cout << setw(BASE_DIGITS) << setfill('0') << digits[i];
}
string toString() const {
string result;
result += to_string(digits.back());
for (int i = digits.size() - 2; i >= 0; i--) {
string segment = to_string(digits[i]);
result += string(BASE_DIGITS - segment.size(), '0') + segment;
}
return result;
}
friend ostream& operator<<(ostream& os, const Big& num) {
os << num.digits.back();
for (int i = num.digits.size() - 2; i >= 0; i--)
os << setw(BASE_DIGITS) << setfill('0') << num.digits[i];
return os;
}
};
\end{lstlisting}
% ===== 正确使用 newpage =====
\newpage % 此处强制分页
\section{搜索}
% 正常内容(删除 moca 123)
\newpage
\section{数学}
\subsection{数论}
\subsubsection{判定质数}
试除法判定某个数是不是质数,时间复杂度 $\mathrm{O(\sqrt{n})}$。
\begin{lstlisting}[caption=判定素数]
bool is_prime(int n) {
if(n < 2)
return false;
for(int i = 2; i <= sqrt(n); i++)
if(n % i == 0)
return false;
return true;
}
\end{lstlisting}
\subsubsection{埃氏筛}
时间复杂度 $\mathrm{O(nloglogn)}$。
\begin{lstlisting}[caption=埃氏筛]
const int N = 1e6 + 5;
int v[N], prime[N], m;
void primes(int n) {
memset(v, 0, sizeof(v));
for(int i = 2; i <= n; i++) {
if(v[i])
continue;
prime[++m] = i;
for(int j = i; j <= n / i; j++)
v[i * j] = 1;
}
}
\end{lstlisting}
\subsubsection{线性筛}
线性素数筛,最后 $\mathrm{prime[1 \sim m]}$ 中保存着 $\mathrm{1 \sim m}$ 范围内的所有素数,$\mathrm{v[i]}$ 保存着 $\mathrm{i}$ 的最小质因子。时间复杂度 $\mathrm{O(n)}$。
\begin{lstlisting}[caption=线性筛]
const int N = 1e6 + 5;
int v[N], prime[N], m;
void primes(int n) {
memset(v, 0, sizeof(v));
m = 0;
for(int i = 2; i <= n; i++) {
if(v[i] == 0)
v[i] = i, prime[++m] = i;
for(int j = 1; j <= m; j++) {
if(prime[j] > v[i] || prime[j] > n / i)
break;
v[i * prime[j]] = prime[j];
}
}
}
\end{lstlisting}
\subsubsection{分解质因数}
时间复杂度 $\mathrm{O(\sqrt{n})}$。
\begin{lstlisting}[caption=分解质因数]
const int N = 1e5 + 5;
int p[N], c[N], m;
void divide(int n) {
m = 0;
for(int i = 2; i <= sqrt(n); i++) {
if(n % i == 0) {
p[++m] = i, c[m] = 0;
while(n % i == 0)
n /= i, c[m]++;
}
}
if(n > 1)
p[++m] = n, c[m] = 1;
}
\end{lstlisting}
\subsubsection{求 n 的正约数集合}
试除法,时间复杂度 $\mathrm{O(\sqrt{n})}$,约数个数上界 $\mathrm{2\sqrt{n}}$。
\begin{lstlisting}[caption=试除法求 n 的正约数集合]
int factor[1600], m = 0;
void divide(int n) {
for(int i = 1; i * i <= n; i++)
if(n % i == 0) {
factor[++m] = i;
if(i != n / i)
factor[++m] = n / i;
}
}
\end{lstlisting}
\subsubsection{求 1 $\sim$ n 的正约数集合}
倍数法,时间复杂度 $\mathrm{NlogN}$。
\begin{lstlisting}[caption=倍数法求 1 $\sim$ n 的正约数集合]
vector<int> factor[500010];
void times(int n) {
for(int i = 1; i <= n; i++)
for(int j = 1; j <= n / i; j++)
factor[i * j].push_back(i);
}
\end{lstlisting}
\subsubsection{整除分块}
解决诸如求 $\sum_{i=1}^nf(i)\lfloor \frac{n}{i} \rfloor$ 一类的问题;首先对于 $\lfloor \frac{n}{i} \rfloor$ 有两条性质:
\textbf{性质 1:} 分块的块数 $\le 2\lfloor \sqrt{n} \rfloor$
\textbf{性质 2:}$i$ 所在块的右端点为 $\lfloor \frac{n}{\lfloor \frac{n}{i} \rfloor} \rfloor$
所以可以先预处理出 $f(i)$ 的前缀和 $s(i) = \sum^i_{j=1} f(j)$,再枚举每一个块 $[l, r]$,累加每块的贡献,时间复杂度 $O(\sqrt{n})$。
\begin{lstlisting}[caption=整除分块]
int div_block(int n) {
int r = 0, ret = 0;
for(int l = 1; l <= n; l = r + 1) {
if(n / l == 0)
break;
r = n / (n / l);
ret += (s[r] - s[l - 1]) * (n / l);
}
return ret;
}
\end{lstlisting}
也可以向二维扩展,求 $\sum_{i=1}^{min(n, m)}f(i) \lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor$,时间复杂度 $O(\sqrt{n} + \sqrt{m})$。
\begin{lstlisting}[caption=整数分块二维扩展]
int div_block(int n, int m) {
int r = 0, ret = 0;
int bound = min(n, m);
for(int l = 1; l <= bound; l = r + 1) {
if(n / l == 0)
break;
// 求两个块的公共区间
r = min(n / (n / l), m / (m / l));
ret += (s[r] - s[l - 1]) * (n / l) * (m / l);
}
return ret;
}
\end{lstlisting}
\subsubsection{欧几里得算法求 gcd}
时间复杂度 $\mathrm{log(a + b)}$。当需要做高精度运算时,建议使用更相减损术代替。
\begin{lstlisting}[caption=欧几里得算法求 gcd]
int gcd(int a, int b) {
return b ? gcd(b, a % b) : a;
}
\end{lstlisting}
\subsubsection{Stein 算法求大整数 gcd}
大整数取模的时间复杂度较高,而加减法时间复杂度较低。针对大整数,我们可以用加减代替乘除求出最大公约数。利用 Stein 算法,时间复杂度 $\mathrm{O(log(n))}$。为了优化常数,可以使用压位高精度。
若 $\mathrm{a = b}$,则 $\mathrm{gcd(a,b) = a}$,否则:
\begin{itemize}
\item 若 $\mathrm{a,b}$ 均为偶数,则 $\mathrm{gcd(a,b) = 2 gcd(\frac{a}{2}, \frac{b}{2})}$
\item 若 $\mathrm{a}$ 为偶数,$\mathrm{b}$ 为奇数,则 $\mathrm{gcd(a,b) = gcd(\frac{a}{2}, b)}$
\item 若 $\mathrm{a,b}$ 同为奇数,则 $\mathrm{gcd(a,b) = gcd(a -b, b)}$
\end{itemize}
\begin{lstlisting}[caption=Stein 算法]
// 基于压位高精度类实现
Big Stein(Big A, Big B) {
int shift = 0;
while(A != 0 && B != 0) {
bool f1 = (A.digits[0] % 2) == 0, f2 = (B.digits[0] % 2) == 0;
if(f1 && f2)
A = A / 2, B = B / 2, shift++;
else if(f1)
A = A / 2;
else if(f2)
B = B / 2;
else {
if(A >= B)
A = A - B;
else
B = B - A;
}
}
Big result = Big(A == 0 ? B : A);
return result * (Big(2) ^ shift);
}
\end{lstlisting}
\subsubsection{裴蜀定理}
对任意整数 $\mathrm{a,b}$,存在一对整数 $\mathrm{x,y}$,$\mathrm{s.t. \ ax + by = gcd(a,b)}$
\noindent \textbf{逆定理:}
设 $\mathrm{a,b}$ 是不全为 $\mathrm{0}$ 的整数,若 $\mathrm{d > 0}$ 是 $\mathrm{a,b}$ 的公因数,且存在整数 $\mathrm{x,y}$,$\mathrm{s.t. \ ax + by = d}$,则 $\mathrm{d = gcd(a,b)}$。
设 $\mathrm{a,b}$ 是不全为 $\mathrm{0}$ 的整数,且存在整数 $\mathrm{x,y}$,$\mathrm{s.t. \ ax + by = 1}$,则 $\mathrm{a,b}$ 互质。
\subsubsection{扩展欧几里得算法}
用于求解 $\mathrm{ax + by = gcd(a,b)}$ 的一组特解。
\begin{lstlisting}[caption=扩展欧几里得算法]
// 返回 a b 的最大公约数 d,并找到一组特解 x0 y0
int exgcd(int a, int b, int &x, int &y) {
if(b == 0) {
x = 1;
y = 0;
return a;
}
int d = exgcd(b, a % b, x, y);
int z = x;
x = y;
y = z - y * (a / b);
return d;
}
\end{lstlisting}
更一般地,对于方程 $\mathrm{ax + by = c}$,它有解当且仅当 $\mathrm{d \mid c}$。我们可以先求出 $\mathrm{ax + by = d}$ 的一组特解 $\mathrm{x_0, y_0}$,然后就得到了 $\mathrm{ax + by = c}$ 的一组特解 $\mathrm{(c/d)x_0, (c/d)y_0}$。则 $\mathrm{ax + by = c}$ 的通解可表示为:
$$
\mathrm{x = \frac{c}{d}x_0 + k\frac{b}{d},\ y = \frac{c}{d}y_0 - k\frac{a}{d}}
$$
其中 $\mathrm{k}$ 取遍整数集合。
\subsubsection{线性同余方程}
求一个整数 $\mathrm{x}$ 满足 $\mathrm{ax \equiv b\ (\bmod m) }$,或者给出无解。
方程可改写为 $\mathrm{ax + my = b}$,设 $\mathrm{d = gcd(a,m)}$ 故线性同余方程有解当且仅当 $\mathrm{d \mid b}$,利用丢番图方程相关知识,知道 $\mathrm{x = x_0 \times b / d}$ 是原线性同余方程的一个特解;方程恰有 $\mathrm{d}$ 个模 $\mathrm{m}$ 不同余的解,为 $\mathrm{x_i = x + i \times (m/d), \ 0 \le i \le d - 1}$。
\subsubsection{乘法逆元}
若 $\mathrm{ax \equiv 1 \ (\bmod\ m)}$,则称 $\mathrm{x}$ 为 $\mathrm{a}$ 模 $\mathrm{m}$ 的逆元,记作 $\mathrm{a^{-1}\ (\bmod \ m)}$。
根据费马小定理,当 $\mathrm{m}$ 为质数时,$\mathrm{b^{m - 2}}$ 即为 $\mathrm{b}$ 的乘法逆元。
只有 $\mathrm{gcd(a,m) = 1}$ 时才存在乘法逆元,可用扩展欧几里得算法求逆元(本质上还是解线性同余方程)。
\begin{lstlisting}[caption=乘法逆元]
// 扩展欧几里得算法求最小正整数逆元
void exgcd(int a,int b, int &x, int &y) {
if(!b) {
x = 1;
y = 0;
return ;
}
exgcd(b, a % b, y, x);
y -= a / b * x;
}
// 求 n mod p 的逆元
// 前提是 n 与 p 互质
int inv(int n, int p) {
int x, y;
exgcd(n, p, x, y);
return (x + p) % p;
}
\end{lstlisting}
\subsubsection{线性求 1 $\sim$ N 的逆元}
\textbf{重要前提:p 是质数。}
显然,1 对 p 的逆元是 1。
对于某个数 i,利用带余除法,有 $p = k \times i + j$,有 $k \times i + j \equiv 0 (\bmod\ p)$,有 $k \times j ^ {-1} + i ^ {-1} \equiv 0 (\bmod\ p)$,即 $i ^ {-1} \equiv -\lfloor \frac{p}{i} \rfloor \times (p \bmod\ i)^{-1} (\bmod\ p)$。
$$
i^{-1} \equiv\left\{\begin{array}{ll}
1, & \text { if } i=1, \\
-\left\lfloor\frac{p}{i}\right\rfloor(p \bmod i)^{-1}, & \text { otherwise. }
\end{array} \quad(\bmod p)\right.
$$
\begin{lstlisting}[caption=线性求逆元]
// 线性求 1 ~ n mod p 的逆元
// 前提是 p 是质数
int get_inv(int n, int p) {
inv[1] = 1;
for (int i = 2; i <= n; ++i) {
inv[i] = (long long)(p - p / i) * inv[p % i] % p;
}
}
\end{lstlisting}
\subsubsection{欧拉函数}
$\mathrm{1 \sim N}$ 中与 $\mathrm{N}$ 互质的数的个数被称为欧拉函数。
若在算数基本定理中,$\mathrm{N = p_1^{c_1}p_2^{c_2}\cdots p_m^{c_m}}$,则
$$
\mathrm{\varphi(N) = N \times \frac{p_1 - 1}{p_1} \times \frac{p_2 - 1}{p_2} \times \cdots \times \frac{p_m - 1}{p_m} = N \times \prod_{prime\ p \mid N} (1 - \frac{1}{p})}
$$
\begin{lstlisting}[caption=欧拉函数]
// 求单个数的 Euler 函数,时间复杂度 O(sqrt(n))
int phi(int n) {
int ans = n;
for(int i = 2; i <= sqrt(n); i++)
if(n % i == 0) {
ans = ans / i * (i - 1);
while(n % i == 0)
n /= i;
}
if(n > 1)
ans = ans / n * (n - 1);
return ans;
}
\end{lstlisting}
\noindent \textbf{欧拉函数的性质}
\begin{enumerate}
\item $\mathrm{\varphi}$ 是积性函数
\item 若 $\mathrm{f}$ 是积性函数,且在算数基本定理中有 $\mathrm{n = p_1^{c_1}p_2^{c_2}\cdots p_m^{c_m}}$,则 $\mathrm{f(n) = \begin{matrix} \prod_{i=1}^m f(p_i^{c_i}) \end{matrix}}$
\item $\mathrm{\forall \ n > 1}$,$\mathrm{1 \sim n}$ 中与 $\mathrm{n}$ 互质的数和为 $\mathrm{n \times \varphi(n) / 2}$
\item 设 $\mathrm{p}$ 为质数,若 $\mathrm{p \mid n}$ 且 $\mathrm{p^2 \mid n}$,则 $\mathrm{\varphi(n) = \varphi(n / p) \times p}$
\item 设 $\mathrm{p}$ 为质数,若 $\mathrm{p \mid n}$ 且 $\mathrm{p^2 \nmid n}$,则 $\mathrm{\varphi(n) = \varphi(n / p) \times (p - 1)}$
\item $\mathrm{\begin{matrix} \sum_{d \mid n} \end{matrix} \varphi(d) = n}$
\item 若 $\mathrm{n = p^k}$,其中 $\mathrm{p}$ 是质数,那么 $\mathrm{\varphi(n) = p^k - p^{k - 1}}$
\item 若 $\mathrm{n}$ 是质数,显然有 $\mathrm{\varphi(n) = n - 1}$
\end{enumerate}
第二条对所有积性函数都适用,后面几条仅对欧拉函数适用。
\subsubsection{求 $2 \sim N$ 中每个数的欧拉函数}
利用线性筛法中,每个数字都会被自己的最小质因子筛掉的性质,以及上面的性质 4、5,进行优化。时间复杂度 $\mathrm{O(n)}$。
\begin{lstlisting}[caption=求 2 $\sim$ N 中每个数的欧拉函数]
int v[maxn], prime[maxn], phi[maxn];
void euler(int n) {
// v 用来记录最小质因子
memset(v, 0, sizeof(v));
m = 0;
for(int i = 2; i <= n; i++) {
if(v[i] == 0) {
v[i] = i, prime[++m] = i;
phi[i] = i - 1;
}
for(int j = 1; j <= m; j++) {
if(prime[j] > v[i] || prime[j] > n / i)
break;
v[i * prime[j]] = prime[j];
phi[i * prime[j]] = phi[i] * (i % prime[j] ? prime[j] - 1 : prime[j]);
}
}
}
\end{lstlisting}
\subsubsection{欧拉定理}
\noindent \textbf{费马小定理}
若 $\mathrm{p}$ 是质数,则 $\mathrm{a^p \equiv a (\bmod\ p)}$,实际上是欧拉定理的特殊情况。
\noindent \textbf{欧拉定理}
若正整数 $\mathrm{a, n}$ 互质,则 $\mathrm{a^{\varphi(n)} \equiv 1(\bmod \ n)}$。
\subsubsection{扩展欧拉定理}
用于降幂:
$$
\mathrm{a^b \equiv \begin{cases}a^{b \bmod \varphi(n)}, & \operatorname{gcd}(a, n)=1, \\ a^b, & \operatorname{gcd}(a, n) \neq 1, b<\varphi(n), \quad(\bmod\ n) \\ a^{(b \bmod \varphi(n))+\varphi(n)}, & \operatorname{gcd}(a, n) \neq 1, b \geq \varphi(n) .\end{cases}}
$$
\begin{lstlisting}[caption=扩展欧拉定理]
// 其中 a 和 n 都是 int;b 是大整数
int ExEuler(int a, Big b, int n) {
int phi_n = phi(n);
if(gcd(a, n) == 1)
return quickPow(a, b % phi_n, n);
else if(b < phi_n)
return quickPow(a, b, n);
return quickPow(a, b % phi_n + phi_n, n);
}
\end{lstlisting}
\subsubsection{中国剩余定理求解线性同余方程组}
设 $\mathrm{m_1, m_2, \cdots , m_n}$ 是两两互质的整数,$\mathrm{m = \begin{matrix} \prod_{i = 1}^n m_i \end{matrix}}$,$\mathrm{M_i =m / m_i}$,$\mathrm{t_i}$ 是线性同余方程 $\mathrm{M_it_i \equiv 1\ (\bmod \ m_i)}$ 的一个解。对于任意的 $\mathrm{n}$ 个整数 $\mathrm{a_1, a_2, \cdots, a_n}$,方程组
$$
\begin{cases}
x \equiv a_1\ (\bmod \ m_1) \\
x \equiv a_2\ (\bmod \ m_2) \\
\vdots \\
x \equiv a_n\ (\bmod \ m_n)
\end{cases}
$$
有整数解,解为 $\mathrm{x = \begin{matrix} \sum_{i = 1}^n \end{matrix} a_iM_it_i}$,在模 $\mathrm{m}$ 意义下有唯一解。
\begin{lstlisting}[caption=中国剩余定理]
// 求同余方程组最小非负整数解
int intChina(int r) {
int Mi, x0, y0, d, ans = 0;
int M = 1;
for(int i = 1; i <= r; i++)
M *= m[i];
for(int i = 1; i <= r; i++) {
Mi = M / m[i];
// 求出每个 ti 的逆元
exgcd(Mi, m[i], x0, y0);
ans = (ans + Mi * x0 * a[i]) % M;
}
return (ans + M) % M;
}
\end{lstlisting}
\subsubsection{扩展中国剩余定理}
中国剩余定理要求 $\mathrm{N}$ 个模数两两互质,其实这个条件表苛刻。当 $\mathrm{N}$ 个模数不满足两两互质时,中国剩余定理不再适用。
可以考虑用数学归纳法,假设已经求出了前 $\mathrm{k-1}$ 个方程构成的方程组的一个解 $\mathrm{x}$。记 $\mathrm{m = lcm(m_1, m_2, \cdots, m_{k - 1})}$,则 $\mathrm{x + i \times m\ (i \in Z)}$ 是前 $\mathrm{k - 1}$ 个方程的通解。
考虑第 $\mathrm{k}$ 个方程,求出一个整数 $\mathrm{t}$,$\mathrm{s.t.\ x + t \times m \equiv a_k \ (\bmod \ m_k)}$,该方程等价于$\mathrm{t \times m \equiv a_k - x \ (\bmod \ m_k)}$,其中 $\mathrm{t}$ 是未知量,这就是一个线性同余方程,可以使用 exgcd 求解。若有解,则 $\mathrm{x’ = x + t \times m}$ 就是前 $\mathrm{k}$ 个方程构成的方程组的一个解。
简单来说,就是使用 $\mathrm{n}$ 次扩展欧几里得算法。
\begin{lstlisting}[caption=扩展中国剩余定理]
// EXCRT,用来解决模数不一定两两互质的情况
// 容易溢出时将数据类型都换为 __int128_t
ll excrt(int n) {
// x 保存前 k - 1 个方程的解, M 保存前 k - 1 个数的 lcm
ll x = a[1] % m[1], M = m[1];
for(int i = 2; i <= n ;i++) {
ll t, y;
ll d = exgcd(M, m[i], t, y);
t = (a[i] - x) / d * t;
x = x + t * M;
M = lcm(M, m[i]);
x = (x % M + M) % M;
}
return x;
}
\end{lstlisting}
\subsubsection{高次同余方程}
给定整数 $\mathrm{a,b,p}$,其中 $\mathrm{a,p}$ 互质,求一个非负整数 $\mathrm{x}$,$\mathrm{s.t. \ a^x \equiv b\ (\bmod \ p)}$。
\begin{lstlisting}[caption=BSGS 求解高次同余方程]
// Baby Step, Giant Step 算法
// 注意:a 与 p 必须互质
// 无解时返回 -1
int BSGS(int a, int b, int p) {
map<int, int> hash;
hash.clear();
b %= p;
int t = (int)sqrt(p) + 1;
for(int j = 0; j < t; j++) {
int val = (long long)b * pow(a, j, p) % p;
hash[val] = j;
}
a = pow(a, t, p);
if(a == 0)
return b == 0 ? 1 : -1;
for(int i = 0; i <= t; i++) {
int val = pow(a, i, p);
int j = hash.find(val) == hash.end() ? -1 : hash[val];
if(j >= 0 && i * t >= j)
return i * t - j;
}
return -1;
}
\end{lstlisting}
\subsubsection{阶乘取模问题}
由中国剩余定理,阶乘取模问题可以转化为模数为素数幂 $p^\alpha$ 的情形,可以对 $n!$ 做如下分解:
$$
n!=p^{\nu_p(n!)}(n!)_p
$$
其中,$\nu_p(n!)$ 表示阶乘 $n!$ 的素因数分解中 $p$ 的幂次,$(n!)_p$ 表示在阶乘 $n!$ 的结果中去除所有 $p$ 的幂次得到的整数。下面将讨论 $(n!)_p$ 在素数幂模下的余数以及幂次 $\nu_p(n!)$ 的具体计算方法。
\noindent \textbf{Wilson 定理}
对于自然数 $n>1$,当且仅当 $n$ 是素数时,$(n-1)!\equiv -1\pmod n$。
推广之后,可以求出 $(n!)_p$。预处理的时间复杂度为 $O(p^\alpha)$,单次询问的时间复杂度为 $O(\log_p n)$。
\begin{lstlisting}[caption=利用推广的 Wilson 定理求 $(n!)_p$]
// Calculate (n!)_p mod pa.
// pa 是某个素数的幂次
int factmod(int n, int p, int pa) {
std::vector<int> f(pa);
f[0] = 1;
for (int i = 1; i < pa; ++i)
f[i] = i % p ? (long long)f[i - 1] * i % pa : f[i - 1];
bool neg = p != 2 || pa <= 4;
int res = 1;
while (n > 1) {
if ((n / pa) & neg)
res = pa - res;
res = (long long)res * f[n % pa] % pa;
n /= p;
}
return res;
}
\end{lstlisting}
\noindent \textbf{Legendre 公式}
Legendre 公式 可以求出 $\mathrm{p}$ 的幂次 $\nu_p(n!)$。它的时间复杂度为 $O(\log n)$。
\begin{lstlisting}[caption=利用 Legendre 公式求 $\nu_p(n!)$]
// Obtain multiplicity of p in n!.
int multiplicity_factorial(int n, int p) {
int count = 0;
do {
n /= p;
count += n;
} while (n);
return count;
}
\end{lstlisting}
利用上面两个公式可以在 $O(logn)$ 时间内求出阶乘的模(只能求出模数为质数的幂次的情景,当模数不是质数幂次时可以用中国剩余定理进行求解)。
\subsection{线性代数}
\subsubsection{矩阵乘法}
矩阵乘法满足结合律、分配率,但不满足交换律。
矩阵相当于状态转移,一维空间状态转移矩阵中 $\mathrm{A_{i,j}}$ 表示状态 $\mathrm{i}$ 如何递推,得到下一时刻的状态 $\mathrm{j}$。同理可扩展到高维空间(但是建议在代码实现层面,将高维空间压缩为一维空间)。
简单来说,状态转移矩阵中每个位置的值由状态 $\mathrm{i}$ 如何转移到状态 $\mathrm{j}$ 决定。
\begin{lstlisting}[caption=矩阵快速幂]
// 矩阵乘法加速线性递推:矩阵快速幂
// 时间复杂度 O(n^3logT),其中 T 是递推总轮数
struct matrix {
ll a[sz][sz];
matrix() {memset(a, 0, sizeof(a));}
matrix operator*(const matrix& T) const {
matrix res;
int r;
for (int i = 0; i < sz; ++i)
for (int k = 0; k < sz; ++k) {
r = a[i][k];
for (int j = 0; j < sz; ++j)
res.a[i][j] += T.a[k][j] * r, res.a[i][j] %= mod;
}
return res;
}
}
matrix matQuickPow(matrix a, int b) {
matrix ret;
for(int i = 0; i < sz; i++)
ret.a[i][i] = 1;
while(b) {
if(b & 1)
ret = ret * a;
a = a * a;
b >>= 1;
}
return ret;
}
\end{lstlisting}
\subsubsection{高斯消元法求解线性方程组}
高斯消元用来求解线性方程组。
线性方程组的所有系数可以写成一个 $\mathrm{m \times n}$ 的系数矩阵,再加上每个方程等号右侧的常数,可以写成一个 $\mathrm{m \times (n + 1)}$ 的增广矩阵, 然后对增广矩阵进行初等行变换:
- 用一个非零的数乘某一行
- 把其中一行的若干倍加到另一行上
- 交换两行的位置
直到将增广矩阵化为一个上三角矩阵,依次回带,得到一个简化阶梯型矩阵,给出了方程组的解。
高斯消元完成后,若存在系数全为零,常数不为零的行,则方程组无解;若主元恰好有 $\mathrm{n}$ 个,方程组有唯一解;若主元有 $\mathrm{k < n}$ 个,方程组有无穷多个解。
高斯消元法就是通过初等行变换把增广矩阵变为简化阶梯型矩阵的线性方程组求解算法。
\begin{lstlisting}[caption=高斯消元法]
// 时间复杂度 O(n^3)
// 无解时输出 -1, 无穷多解时输出 0
void gaussElimination(int n) {
int nwline = 0;
for(int k = 0; k < n; k++) {
// 用于行主元选取
int maxRow = nwline;
for(int i = nwline + 1; i < n; i++)
if(fabs(p[i][k]) > fabs(p[maxRow][k]))
maxRow = i;
if(fabs(p[maxRow][k]) < eps)
continue;
// 交换当前行和最大元素所在行
for(int i = 0; i <= n; i++)
swap(p[nwline][i], p[maxRow][i]);
for(int i = 0; i < n; i++) {
if(i == nwline)
continue;
double mul = p[i][k] / p[nwline][k];
for(int j = k; j <= n; j++)
p[i][j] -= p[nwline][j] * mul;
}
nwline++;
}
// 存在找不到主元的情况
// 一定是无解或者无穷多解,优先判断无解
if(nwline < n) {
while(nwline < n)
if(!(fabs(p[nwline++][n]) < eps)) {
printf("-1");
return ;
}
printf("0");
return ;
}
for(int i = 0; i < n; i++)
printf("%.2lf\n", fabs(p[i][n] / p[i][i]) < eps ? 0 : p[i][n] / p[i][i]);
}
\end{lstlisting}
\subsubsection{高斯消元法求解异或方程组}
$$