目录
  1. 1. 写在前面
  2. 2. 常用C++STL(模板)合集
    1. 2.1. 万能(误)算法头文件(部分)
    2. 2.2. 动态数组(vector)、双向链表(list)
    3. 2.3. 普通队列、双端队列、优先队列
    4. 2.4.
    5. 2.5. pair(成组)、set(有序元素序列)
    6. 2.6. map(映射)
    7. 2.7. 一些C++的功能/特性
    8. 2.8. 强大的pb_ds库
  3. 3. 数据结构部分
    1. 3.1. 单调队列
    2. 3.2. 单调栈
    3. 3.3. 归并排序、逆序对数
    4. 3.4. 线段树/树状数组
    5. 3.5. 可持久化线段树(主席树)
      1. 3.5.1. 朴素版本
      2. 3.5.2. 简约版本
    6. 3.6. 归并树(划分树)
    7. 3.7. ST表(倍增)求区间最值
      1. 3.7.1. 一维
      2. 3.7.2. 二维
    8. 3.8. 字典树
    9. 3.9. 树链剖分
      1. 3.9.1. 点权
      2. 3.9.2. 边权
    10. 3.10. 并查集
      1. 3.10.1. 普通并查集
      2. 3.10.2. 带权并查集
  4. 4. 数论部分
    1. 4.1. 一些基础的公式、定理
    2. 4.2. 因数、余数、质数基础
      1. 4.2.1. 埃氏素数筛法
      2. 4.2.2. 欧拉函数
      3. 4.2.3. 欧几里得定理求a、b最大公因数、最小公倍数
      4. 4.2.4. 扩展欧几里得
      5. 4.2.5. 线性同余方程
      6. 4.2.6. 模线性方程组求解(中国剩余定理)
      7. 4.2.7. 乘法逆元
      8. 4.2.8. 二次剩余
      9. 4.2.9. 普通分解质因数
    3. 4.3. 约数的和取余
    4. 4.4. 组合数相关
      1. 4.4.1. 一些定理
      2. 4.4.2. 预处理组合数
      3. 4.4.3. 求单个组合数
      4. 4.4.4. 卢卡斯定理
      5. 4.4.5. 斯特林数
    5. 4.5. Miller_Rabin 快速大素数测试
    6. 4.6. Pollard_rho大数分解质因数
    7. 4.7. 欧拉降幂公式
    8. 4.8. BSGS大步小步算法
    9. 4.9. 莫比乌斯函数、反演
    10. 4.10. 牛顿迭代法求平方根
    11. 4.11. 高斯消元
    12. 4.12. 线性基
    13. 4.13. 斐波那契数列
  5. 5. 图论部分
    1. 5.1. 邻接表存图
    2. 5.2. Vector伪邻接表存图
    3. 5.3. Dijkstra单源最短路
    4. 5.4. SPFA单源最短路
    5. 5.5. Floyd多源最短路
    6. 5.6. 最小生成树MST
      1. 5.6.1. Kruskal
      2. 5.6.2. Prim
    7. 5.7. 最近公共祖先(LCA)
    8. 5.8. 二分图匹配
    9. 5.9. 网络流
    10. 5.10. 树分治(点分治)
  6. 6. 字符串处理部分
    1. 6.1. KMP
    2. 6.2. exKMP
  7. 7. 博弈论部分
    1. 7.1. nim博弈
    2. 7.2. Bash博弈
    3. 7.3. Wythoff博弈
    4. 7.4. 斐波那契博弈
    5. 7.5. 环形博弈
    6. 7.6. SG函数
  8. 8. 动态规划部分
    1. 8.1. 背包
    2. 8.2. 区间DP
  9. 9. 计算几何部分
    1. 9.1. Kuangbin几何模板
  10. 10. 杂项
    1. 10.1. 竞赛常用快速读入骚操作
    2. 10.2. C++高精度模板
    3. 10.3. Java大整数
    4. 10.4. Java大实数
    5. 10.5. 手动扩栈
    6. 10.6. 对拍
个人模板整理

写在前面

本篇博客主要用于本人的ICPC比赛用模板整理,方便查阅。

如有网友发现哪里写的不好、哪里可以改进,欢迎提出。

持续更新中…

常用C++STL(模板)合集

注意:C++STL全部需要使用命名空间std

  • 万能(误)算法头文件(部分)

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
#include <algorithm>
using namespace std;
int main() {
iterator begin, end; // 指代某种数据结构首尾迭代器
T i, x, a, b;

sort(begin, end, <cmp>); // 排序函数,默认从小到大
// 遇到需要特殊排序的需要编写cmp函数 or 重载内部运算符
next_permutation(begin, end); // 下一个排列
prev_permutation(begin, end); // 前一个排列

set_union(begin(a), end(a), begin(b), end(b), begin(c));
// 取两个有序序列a、b的并集,存放到c中
set_intersection(begin(a), end(a), begin(b), end(b), begin(c));
// 取两个有序序列a、b的交集,存放到c中
set_difference(begin(a), end(a), begin(b), end(b), begin(c));
// 取两个有序序列a、b的差集,存放到c中
unique(begin, end); // 有序数据去重
merge(begin(a), end(a), begin(b), end(b), begin(c), cmp);
// 合并两个有序序列a、b,存放到c中,cmp可定义新序列排列方式

lower_bound(begin, end, x); // 返回x的前驱迭代器
// 在普通的升序序列中,x的前驱指的是第一个大于等于x的值
upper_bound(begin, end, x); // 返回x的后继迭代器
// 在普通的升序序列中,x的后继指的是第一个大于x的值
// 上述两个函数时间复杂度为O(log2(n)),内部实现是二分
// 如果找不到这样的值,会返回end

find(begin, end, x); // O(n)查找x
binary_search(begin, end, x) // 二分查找x,返回bool
min(a, b); max(a, b); // 返回a、b中的最小/最大值
fill(begin, end, x); // 往容器的[begin, end)内填充x
swap(a, b); // 交换a、b的值

return 0;
}
  • 动态数组(vector)、双向链表(list)

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
#include <vector>
#include <list>
using namespace std;
int main() {
T i;
unsigned int n, x;
bool flag;
iterator it;

// 动态数组部分
// 注意vector的空间需要预留两倍大小
vector<T> v;
v.push_back(i); // 往数组尾添加一个元素i
v[x]; // 访问第x - 1个元素
v.begin(); // 返回头元素的迭代器
v.end(); // 返回末尾迭代器(尾元素的下一个)
n = v.size(); // 数组中元素数量
v.pop_back(); // 删除最后一个元素
v.erase(it); // 删除某个的元素
v.insert(x, i); // 在x位置插入元素i
// erase、insert时间复杂度为O(n)
v.clear(); // 清空数组,不释放空间
flag = v.empty(); // 判断数组是否为空(真值)

// 链表部分
list<T> li;
li.push_front(i); // 在链头添加一个元素i
li.push_back(i); // 在链尾添加一个元素i
li.pop_front(i); // 删除链表头元素
li.pop_back(i); // 删除链表尾元素
li.erase(it); // 删除某个的元素
li.insert(x, i); // 在x位置插入元素i O(n)
li.begin(); // 返回头元素的迭代器
li.end(); // 返回末尾迭代器(尾元素的下一个)
n = li.size(); // 链表中元素数量
li.remove(i); // 删除链表中所有值为i的元素
li.unique(); // 移除所有连续相同元素,留下一个
li.reverse(); // 反转链表
li.clear(); // 情况链表,不释放空间

return 0;
}
  • 普通队列、双端队列、优先队列

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
#include <queue> // 队列头文件
#include <deque> // 双端队列头文件
using namespace std;
int main() {
T i;
unsigned int n, x;
bool flag;

// 普通队列部分,注意queue没有迭代器
queue<T> q, tmp_q; // 定义普通队列
q.push(i); // 队尾插入元素i
q.pop(); // 弹出队首元素
i = q.front(); // 访问队首元素
i = q.back(); // 访问队尾元素
n = q,size(); // 队内元素数量
flag = q.empty(); // 判断队列是否为空(真值)
q.swap(tmp_q); // 交换两个队列元素

// 优先队列部分,注意其没有迭代器
priority_queue<T> pq; // 定义优先队列
pq.push(i); // 队尾插入元素i
pq.pop(); // 弹出队首元素
i = pq.top(); // 访问队首元素
n = q,size(); // 队内元素数量
flag = q.empty(); // 判断队列是否为空(真值)
q.swap(tmp_q); // 交换两个队列元素
// 注意优先队列内部是使用<运算符,默认大根堆
// 可以采用重载运算符或加入运算符类自定义排列方式
// 例:priority_queue<T, vector<T>, greater<T> > 小根堆
/*
struct node {
int x, y;
};
bool operator < (node a, node b) {
// 这里注意是<右边的元素会放在前面
if(a.x != b.x) return a.x < b.x;
else return a.y < b.y;
}
priority_queue<node>
*/

// 双端队列部分
// 注意deque用到了map来映射,时间复杂度上常数略大
deque<T> dq; // 定义双端队列
// 可以称为vector、list、queue的结合体
// 用法类似,这里只给代码不做注释
dq.push_back(i);
dq.push_front(i);
dq.front();
dq.back();
dq.pop_front();
dq.pop_back();
dq.begin();
dq.end();
dq[x];
n = dq.size();
flag = dq.empty();
dq.insert(x, i);

return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include <stack>
using namespace std;
int main() {
T i;
unsigned int n;
bool flag;

stack<T> st; // 注意stack没有迭代器
st.push(i); // 往栈顶加入一个元素
st.pop(); // 弹出栈顶元素
i = st.top(); // 获得栈顶元素的值
flag = st.empty(); // 判断是否为空(真值)
n = st.size(); // 获得栈内元素个数

return 0;
}
  • pair(成组)、set(有序元素序列)

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
#include <set>
#include <pair>
using namespace std;
int main() {
T i;
T1 t1;
T2 t2;
iterator it;
unsigned int n;
bool flag;

// pair是将两种元素组成一对
pair<T1, T2> p;
p = make_pair(t1, t2); // 将t1、t2构造成一对
// pair支持比较,遵循字典序
p.first; // 访问第一个元素,这里是t1
p.second; // 访问第二个元素,这里是t2

// set内部是RB-tree维护
set<int> st; // 注意,set内元素不重复
st.insert(i); // 往set内插入一个元素i
// 时间复杂度O(log2(n)) 这里会返回一个<pair>迭代器
// first指向插入元素后所在的迭代器
// second指向是否插入成功(真值)
st.begin(); // 返回首迭代器
st.end(); // 范围尾迭代器
st.erase(it); st.erase(i);
// 删除某个元素
st.equal_range(i); // 返回几何中与i相等的上下限两个迭代器
flag = st.empty();
n = st.size();
st.clear();
// set内置了lower_bound和upeer_bound函数
// 用法和algorithm的一样

// 可重复元素set
multiset<int> mst;
// 用法与set大致相同
// 唯一不同只在删除函数上
mst.erase(i); // 会删除所有值为i的元素

return 0;
}
1
2
3
4
5
6
7
// 如果需要给set自定义排序顺序
struct CMP {
bool operator() (const int& a, const int& b) const {
return a > b; // 返回真值则代表左边的值优先级高
}
};
multiset<int, CMP> mst;
  • map(映射)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#include <map>
using namespace std;
int main() {
T1 t1;
T2 t2;

// map将两种元素做映射,一种指向另一种
// 内部也是RB-tree维护
map<T1, T2> mp;
mp[t1] = t2; // 直接让t1对应到t2
mp[t1]; // 访问t1对应的内容,时间复杂度O(log2(n))
// 如果t1没有指向任何内容,则会返回T2类型的初始值

return 0;
}
  • 一些C++的功能/特性

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
#include <bits/stdc++.h> // 标准库头文件
using namespace std;
int main() {

__int128 a; // 128位整数,最大值大概10^38次方
// C++11以上可用,无法用标准方法读入

cin.tie(0); cout.tie(0);
ios::sync_with_stdio(false);
// 关闭cin、cout同步流,此举后不可混用scanf/printf

auto x; // 自动变量,可以是任意属性
// 举个例子
std::set<int> st;
std::for(auto i:st); // C++版for_each
// 其中i是auto变量,也可改成set<int>::iterator

// 所有的STL容器push/insert操作都可替换为emplcace
// 速度上减小常数(不用临时变量)
// 例:
int i;
std::set<int> st;
st.emplace(i);
std::vector<int> vc;
vc.emplace_back(i);

return 0;
}
  • 强大的pb_ds库

这部分并不是标准库STL,是扩展模板库。

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
#include <bits/stdc++.h>
#include <bits/extc++.h> // 扩展库头文件
/*
* 这里如果没有bits/extc++.h的话需要
* ext/pb_ds/tree_policy.hpp
* ext/pb_ds/assoc_container.hpp
* ext/pb_ds/priority_queue_policy.hpp
* ext/pb_ds/trie_policy.hpp
* ext/rope
* ...
*/
using namespace __gnu_pbds;
using namespace __gnu_cxx;
int main() {

// 哈希表部分,用法与map一样,效率在C++11以下效率高
// 注意,这部分在namespace __gnu_pbds下
cc_hash_table<string, int> mp1; // 拉链法
gp_hash_table<string, int> mp2; // 查探法(快一些)

// 优先队列部分,比STL中高级
priority_queue<int, std::greater<int>, TAG> pq;
/*
* 第一个参数是数据类型
* 第二个是排序方式
* 第三个是堆的类型
* 其中堆的类型有下面几种
* pairing_heap_tag
* thin_heap_tag
* binomial_heap_tag
* c_binomial_heap_tag
* binary_heap_tag
* 其中pairing_heap_tag最快
* 并且这个东西是带默认参数的,只需要定义一个int
*/
// 比STL中的优先队列多了join和迭代器
// 例子:
priority_queue<int> pq1, pq2;
pq1.join(pq2); // 会将pq2合并到pq1上
pq1.begin(); pq1.end(); // 可遍历

// 红黑树(平衡树)部分,与set相似,但更快
tree <
int,
null_type,
std::less<>,
rb_tree_tag,
tree_order_statistics_node_update
> t, tre;
/*
* int 关键字类型
* null_type 无映射(低版本g++为null_mapped_type)
* less<int> 从小到大排序
* rb_tree_tag 红黑树(splay_tree_tag splay)
* tree_order_statistics_node_update 结点更新
*/
int i, k;
t.insert(i); // 插入
t.erase(i); // 删除
t.order_of_key(i);
// 询问这个tree中有多少个比i小的元素
t.find_by_order(k);
// 找第k + 1小的元素的迭代器,如果order太大会返回end()
t.join(tre); // tre合并到t上
t.split(i, tre); // 小于等于i的保留,其余的属于tre
// 基本操作有size()/empty()/begin()/end()等
// 同样内置lower_bound/upper_bound

// 可持久化平衡树部分
// 注意,这部分在namespace __gun_cxx下
rope<char> str;
// 待我学习完后再更新

return 0;
}

数据结构部分

  • 单调队列

1
// 待更新
  • 单调栈

可以O(n)O(n)地解决一些区间长度与区间最值关系问题。

例如:求nn个数中max[L,R]n{mini[L,R]{ai}(RL+1)}\displaystyle\max_{[L, R] \in n}\{\min_{i\in[L, R]}\{a_i\} \cdot (R-L+1)\}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#define ll long long
const int MAX_N = 2e6 + 10;

struct data {
ll a, id;
}a[MAX_N];

stack<data> st;

int main() {
ll n, ans = 0;
data last;
scanf("%d", &n);
for(int i = 1; i <= n; ++i) {
scanf("%d", &a[i].a);
a[i].id = i;
if(st.empty()) st.push(a[i]);
else {
while(!st.empty() && st.top().a > a[i].a) {
last = st.top();
st.pop();
if(st.empty())
ans = max((i - 1) * last.a, ans);
else ans = max((i - st.top().id - 1) * last.a, ans);
} st.push(a[i]);
}
}
while(!st.empty()) {
last = st.top();
st.pop();
if(st.empty()) ans = max(n * last.a, ans);
else ans = max((n - st.top().id) * last.a, ans);
}
printf("%lld\n", ans);
return 0;
}
  • 归并排序、逆序对数

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
// 归并排(升)序,p为逆序对数
#define mm (l + r) >> 1
#define lson l, mid
#define rson mid + 1, r
const int MAX_N = 1e5 + 10;
int a[MAX_N], tmp[MAX_N];
void m_sort(int l, int r, int &p) {
if(l == r) return ;
int mid = mm;
m_sort(lson, p), m_sort(rson, p);
int i = l, j = mid + 1, k = l;
while(i <= mid && j <= r) {
if(a[i] <= a[j])
tmp[k++] = a[i++];
else {
tmp[k++] = a[j++];
p += mid - i + 1;
}
}
while(i <= mid) tmp[k++] = a[i++];
while(j <= r) tmp[k++] = a[j++];
for(i = l; i <= r; ++i) a[i] = tmp[i];
}
/*
* 逆序数的应用
* 1.假设序列中任意相邻数可以交换,逆序对数是让序列升序的最少交换次数
* 2.设f[n][m]代表逆序对数为m的n元排列个数,则有
* f[n + 1][m] = f[n][m] + f[n][m - 1] + …… + f[n][m - n]
* 化简得f[n + 1][m] = f[n][m] + f[n + 1][m - 1]
*/
  • 线段树/树状数组

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define lowbit(x) (x & -x)
#define mm (l + r) >> 1
#define lson l, mid, p << 1 // 左儿子
#define rson mid + 1, r, p << 1|1 // 右儿子
const int MAX_N = 200000 + 10; // 数量大小

ll sum[MAX_N], tree[MAX_N << 2], lazy[MAX_N << 2], a[MAX_N];

inline void push_up(int p) {
tree[p] = tree[p << 1] + tree[p << 1|1];
}

inline void push_down(int l, int r, int p) {
if(lazy[p]) {
int mid = mm;
// 如果当前区间有lazy标记,则需要下放
tree[p << 1] += (mid - l + 1) * lazy[p];
tree[p << 1|1] += (r - mid) * lazy[p];
// 先更新两个儿子节点的值
lazy[p << 1] += lazy[p];
lazy[p << 1|1] += lazy[p];
// 标记下放到两个儿子节点
lazy[p] = 0; // 当前节点标记取消
}
}

void build(int l, int r, int p) {
if(l == r) {
tree[p] = a[l]; // 叶子节点信息
lazy[p] = 0;
return ;
}
int mid = mm;
build(lson), build(rson); // 分别往左往右
push_up(p); // 将儿子节点信息上传
}

ll query(int l, int r, int p, int L, int R) {
if(L <= l && r <= R) return tree[p];
// 区间包含则返回该节点所有信息
int mid = mm;
ll ret = 0; // ret用来保存信息
push_down(l, r, p); // 标记下放
if(L <= mid) ret += query(l, mid, p << 1, L, R);
// 取左区间的部分信息
if(R > mid) ret += query(mid + 1, r, p << 1|1, L, R);
// 取右区间的部分信息
return ret; // 将以上结果返回
}

void update(int l, int r, int p, int L, int R, int s) {
if(L <= l && r <= R) {
// 如果当前访问区间被要访问的区间所包含
tree[p] += (r - l + 1) * s;
// 首先更新当前节点的值
lazy[p] += s; // 打上lazy标签
return ;
}
int mid = mm;
push_down(l, r, p); // 标记下放
if(L <= mid) update(l, mid, p << 1, L, R, s);
// L <= mid 说明需要更新左区间
if(R > mid) update(mid + 1, r, p << 1 | 1, L, R, s);
// R > mid 说明需要更新右区间
tree[p] = tree[p << 1] + tree[p << 1|1];
// 别忘了要更新 p 节点的信息哦~
}

ll ask(int x) {
ll ret = 0;
for(; x; x -= lowbit(x)) ret += sum[x];
return ret;
}
ll get_sum(int l, int r) {
// 询问区间[l, r]的和
return ask(r) - ask(l - 1);
}
void change(int n, int x, ll s) {
// 给x位置上的值加上s
for(; x <= n; x += lowbit(x))
sum[x] += s;
}
  • 可持久化线段树(主席树)

    主席树本质上是n棵权值线段树,本质上不可修改,套上树状数组后可以支持单点修改(待更新)

  • 朴素版本

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
// 代码是求区间第k小的值
// 大部分情况需要离散化
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define mm (l + r) >> 1
#define lson l, mid
#define rson mid + 1, r
const int MAX_N = (int)1e5 + 10;

struct node {
int ls, rs; // 用来存左右儿子下标
ll sum; // 数字和
}tree[MAX_N * 40];

int rt[MAX_N], CNT, a[MAX_N], p[MAX_N];
vector<int> vt;
// re是各版本线段树根节点下标 CNT是新建节点下标
// a是原数组,p是离散化后对应下标

int disc(int n) { // 离散化
for(int i = 1; i <= n; ++ i)
vt.push_back(a[i]);
sort(vt.begin(), vt.end());
unique(vt.begin(), vt.end());
for(int i = 1; i <= n; ++ i)
p[i] = lower_bound(vt.begin(), vt.end(), a[i]) - vt.begin() + 1;
}

int build(int l, int r) {
// 建立空树
int pos = ++ CNT;
if(l == r) return pos;
int mid = mm;
tree[pos].ls = build(lson); // 往左建树
tree[pos].rs = build(rson); // 往右建树
return pos;
}

inline void push_up(int pos) {
tree[pos].sum = tree[ tree[pos].ls ].sum + tree[ tree[pos].rs ].sum;
}

int update(int last, int l, int r, int k) {
// 更新第k点的值
int pos = ++ CNT;
if(l == r) {
tree[pos].sum = tree[last].sum + 1;
return pos;
}
int mid = mm;
tree[pos].ls = tree[last].ls, tree[pos].rs = tree[last].rs;
if(k <= mid)
tree[pos].ls = update(tree[last].ls, lson, k);
else
tree[pos].rs = update(tree[last].rs, rson, k);
push_up(pos);
return pos;
}


int query(int posL, int posR, int l, int r, int k) {
if(l == r) return l;
int cnt = (int)(tree[ tree[posR].ls ].sum - tree[ tree[posL].ls ].sum);
int mid = mm;
if(k <= cnt) return query(tree[posL].ls, tree[posR].ls, lson, k);
else return query(tree[posL].rs, tree[posR].rs, rson, k - cnt);
}

void init(int n) {
CNT = 0;
build(1, n);
for(int i = 1; i <= n*40; ++ i)
tree[i].sum = tree[i].ls = tree[i].rs = 0;
vt.clear();
}

int main() {
int n, q, i, l, r, k;

while(~ scanf("%d %d", &n, &q)) {
init(n); // 初始化
for(i = 1; i <= n; ++ i) scanf("%d", &a[i]);
disc(n);
for(i = 1; i <= n; ++ i)
rt[i] = update(rt[i - 1], 1, n, p[i]);
while(q --) {
scanf("%d %d %d", &l, &r, &k);
printf("%d\n", vt[query(rt[l - 1], rt[r], 1, n, k) - 1]);
}
}

return 0;
}
  • 简约版本

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
// 简写主席树模板,没有建初始树过程
// 可以不要离散化,要注意空间
#include <bits/stdc++.h>
using namespace std;
const int MAX_N = 1e5 + 10;
#define ll long long
#define mm (l + r) >> 1
#define lson l, mid
#define rson mid + 1, r
struct node {
int ls, rs; // 用来存左右儿子下标
int sum; // 数字个数和
void copy(const node tmp) {
ls = tmp.ls;
rs = tmp.rs;
sum = tmp.sum;
}
}tree[MAX_N * 40];
int CNT, rt[MAX_N];
ll a[MAX_N];

void init(int n) {
CNT = 0;
for(int i = 1; i <= n*40; ++i)
tree[i].sum = tree[i].ls = tree[i].rs = 0;
for(int i = 0; i <= n; ++i) rt[i] = 0;
}

inline void push_up(int pos) {
tree[pos].sum = tree[ tree[pos].ls ].sum + tree[ tree[pos].rs ].sum;
}

void update(int &pos, int pre, ll l, ll r, ll val) {
pos = ++CNT;
tree[pos].copy(tree[pre]);
if(l == r) {
++tree[pos].sum;
return;
}
ll mid = mm;
if(val <= mid) update(tree[pos].ls, tree[pre].ls, lson, val);
else update(tree[pos].rs, tree[pre].rs, rson, val);
push_up(pos);
}

int query(int posL, int posR, ll l, ll r, int val) {
if(l == r) return l;
int cnt = (int)(tree[ tree[posR].ls ].sum - tree[ tree[posL].ls ].sum);
ll mid = mm;
if(val <= cnt) return query(tree[posL].ls, tree[posR].ls, lson, val);
else return query(tree[posL].rs, tree[posR].rs, rson, val - cnt);
}
  • 归并树(划分树)

    本质上是把归并排序和二叉树思想结合,可以解决部分主席树可以解决的东西,但不能回滚和修改。

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
// 求区间第k小
#include <bits/stdc++.h>
using namespace std;
#define mm (l + r) >> 1
#define lson l, mid, dep + 1
#define rson mid + 1, r, dep + 1
const int MAX_N = 1e5 + 10;

int a[MAX_N], tree[32][MAX_N];

void build(int l, int r, int dep) {
if(l > r) return ;
if(l == r) {
tree[dep][l] = a[l];
return ;
}
int mid = mm, i = l, j = mid + 1, k = l;
build(lson), build(rson);
while(i <= mid && j <= r) {
if(tree[dep + 1][i] <= tree[dep + 1][j])
tree[dep][k ++] = tree[dep + 1][i ++];
else
tree[dep][k ++] = tree[dep + 1][j ++];
}
while(i <= mid)
tree[dep][k ++] = tree[dep + 1][i ++];
while(j <= r)
tree[dep][k ++] = tree[dep + 1][j ++];
return ;
}

int query(int l, int r, int dep, int ql, int qr, int key) {
if(qr < l || r < ql) return 0;
if(ql <= l && r <= qr)
return lower_bound(&tree[dep][l], &tree[dep][r] + 1, key) - &tree[dep][l];
int mid = mm, ret = 0;
if(qr > mid) ret += query(rson, ql, qr, key);
if(ql <= mid) ret += query(lson, ql, qr, key);
return ret;
}

int main() {
int n, m, i, ql, qr, k, l, r, mid, tot;
while(~scanf("%d %d", &n, &m)) {
for(i = 1; i <= n; ++i)
scanf("%d", &a[i]);
build(1, n, 1);
while(m --) {
scanf("%d %d %d", &ql, &qr, &k);
l = 1, r = n + 1, mid = mm;
while(l + 1 < r) {
mid = mm;
tot = query(1, n, 1, ql, qr, tree[1][mid]);
if(tot <= k - 1) l = mid;
else r = mid;
}
printf("%d\n", tree[1][l]);
}
}
return 0;
}
  • ST表(倍增)求区间最值

  • 一维

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
// 一维区间求最大值
const int MAX_N = 1e5 + 10;
int dp[MAX_N][20], vm[MAX_N];
/*
需要改成存最小(大)值下标,可以增加下列函数
int Min(int x, int y) {
return val[x] <= val[y] ? x : y;
}
*/
void init_RMQ(int n, int b[]) {
vm[0] = -1;
for(int i = 1; i <= n; ++i) {
vm[i] = ((i & (i-1)) == 0) ? vm[i - 1]+1 : vm[i - 1];
dp[i][0] = b[i];
}
for(int j = 1; j <= vm[n]; ++j)
for(int i = 1; i + (1<<j) - 1 <= n; ++i)
dp[i][j] = max(dp[i][j - 1], dp[i + (1<<(j-1)) ][j - 1]);
}
int RMQ(int l, int r) {
int k = vm[r - l + 1];
return max(dp[l][k], dp[r - (1<<k) + 1][k]);
}
  • 二维

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
// 二维RMQ,预处理复杂度O(n * m * log2(n) * log2(m))
const int MAX_N = 310;
int val[MAX_N][MAX_N];
int dp[MAX_N][MAX_N][9][9];
int vm[MAX_N]; // 二进制位数-1,使用前需初始化

void init_VM() {
vm[0] = -1;
for(int i = 1; i < MAX_N; ++i)
vm[i] = ((i & (i-1)) == 0) ? vm[i - 1]+1 : vm[i - 1];
}

void init_RMQ(int n, int m) {
for(int i = 1; i <= n; ++i)
for(int j = 1; j <= m; ++j)
dp[i][j][0][0] = val[i][j];
for(int ii = 0; ii <= vm[n]; ++ii)
for(int jj = 0; jj <= vm[m]; ++jj)
if(ii + jj)
for(int i = 1; i + (1<<ii) - 1 <= n; ++i)
for(int j = 1; j + (1<<jj) - 1 <= m; ++j) {
if(ii) dp[i][j][ii][jj] = max(dp[i][j][ii - 1][jj],
dp[i + (1<<(ii-1)) ][j][ii - 1][jj]);
else dp[i][j][ii][jj] = max(dp[i][j][ii][jj - 1],
dp[i][j + (1<<(jj-1)) ][ii][jj - 1]);
}
}

int RMQ(int x1, int y1, int x2, int y2) {
// 求子矩阵最大值
int k1 = vm[x2 - x1 + 1];
int k2 = vm[y2 - y1 + 1];
x2 = x2 - (1<<k1) + 1;
y2 = y2 - (1<<k2) + 1;
return max( max(dp[x1][y1][k1][k2], dp[x1][y2][k1][k2]),
max(dp[x2][y1][k1][k2], dp[x2][y2][k1][k2]) );
}
  • 字典树

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
// k叉字典树,空间占用较大
struct Kx_Trie {
static const int MAX_N = 1010; // 节点数
static const int MAX_M = 26; // 字符数
static const int Tfirst = 'a'; // 第一个字符

int CNT;
struct node {
int son[MAX_M], sum; // sum是前缀数量
bool flag; // true代表这是个完整的单词
char c;
void init() {
sum = flag = c = 0;
memset(son, 0, sizeof(son));
}
}tree[MAX_N];

void init() {
CNT = 0;
tree[0].init();
}

void insert(char* str) {
int rt, nxt;
for(rt = 0; *str; rt = nxt, ++str) {
nxt = tree[rt].son[*str - Tfirst];
if(nxt == 0) {
nxt = tree[rt].son[*str - Tfirst] = ++CNT;
tree[CNT].init();
tree[CNT].c = *str;
}
++tree[nxt].sum;
} tree[rt].flag = true;
}

void erase(char* str) {
// 假设已经插入str
int rt = 0;
for(; *str; ++str) {
--tree[rt].sum;
rt = tree[rt].son[*str - Tfirst];
}
--tree[rt].sum;
tree[rt].flag = false;
}

bool search(char* str) {
// 查询是否有这个单词
int rt = 0;
for(; *str; ++str) {
rt = tree[rt].son[*str - Tfirst];
if(!rt) return false;
}
return tree[rt].flag;
}

int prefix(char *str) {
int rt, len;
for(rt = len = 0; *str; ++str, ++len) {
rt = tree[rt].son[*str - Tfirst];
if(!rt || !tree[rt].sum) break ;
}
return len;
}
};
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
// 左孩子右兄弟字典树,时间换空间
struct Ls_Rb_Trie {
static const int MAX_N = 1010; // 节点数
int top;
struct node {
char c;
int l, r, sum;
bool flag;
void init() {
l = r = sum = flag = c = 0;
}
}tree[MAX_N];

void init() {
top = 0;
memset(tree, 0, sizeof(tree));
}

bool search(char* str) {
int rt;
for(rt = 0; *str; ++str) {
for(rt = tree[rt].l; rt; rt = tree[rt].r) {
if(tree[rt].c == *str) break ;
}
if(!rt || !tree[rt].sum) return false;
}
return tree[rt].flag;
}

void insert(char* str) {
int i, rt;
for(rt = 0; *str; ++str, rt = i) {
for(i = tree[rt].l; i; i = tree[i].r)
if(tree[i].c == *str)
break ;
if(i == 0) {
tree[rt].l = i = ++top;
tree[top].init();
tree[top].r = tree[rt].l;
tree[top].c = *str;
}
++tree[i].sum;
}
tree[rt].flag = true;
}

void erase(char* str) {
int rt;
for(rt = 0; *str; ++str) {
for(rt = tree[rt].l; rt; rt = tree[rt].r) {
if(tree[rt].c == *str)
break;
}
--tree[rt].sum;
}
tree[rt].flag = 0;
}

int prefix(char* str) {
// 最长前缀,
int rt = 0, lv;
for(lv = 0; *str; ++str, ++lv) {
for(rt = tree[rt].l; rt; rt = tree[rt].r)
if(tree[rt].c == *str)
break ;
if(rt == 0 || !tree[rt].sum) break ;
}
return lv;
}
};
  • 树链剖分

  • 点权

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
// 查询单点值,修改路径上的点权,HDU3966
const int MAX_N = 5e4 + 10;
struct node {
int to, go;
}edge[MAX_N << 1];
int head[MAX_N], tot, pos;
int top[MAX_N]; // top[v]表示v所在重链的顶端节点
int fa[MAX_N]; // 父亲节点
int deep[MAX_N]; // 深度
int num[MAX_N]; // num[v]表示v为根的子树节点数量
int p[MAX_N]; // p[v]表示v对应的位置
int fp[MAX_N]; // fp和p数组相反
int son[MAX_N]; // 重儿子

void init() {
tot = 0;
memset(head, -1, sizeof(head));
pos = 1; // 树状数组编号
memset(son, -1, sizeof(son));
}

void add_edge(int u, int v) {
edge[tot].to = v, edge[tot].go = head[u], head[u] = tot++;
}

void dfs(int u, int pre, int dep) {
deep[u] = dep;
fa[u] = pre;
num[u] = 1;
for(int i = head[u]; i != -1; i = edge[i].go) {
int v = edge[i].to;
if(v != pre) {
dfs(v, u, dep+1);
num[u] += num[v];
if(son[u] == -1 || num[v] > num[son[u]])
son[u] = v;
}
}
}

void get_pos(int u, int sp) {
top[u] = sp;
p[u] = pos++;
fp[p[u]] = u;
if(son[u] == -1) return ;
get_pos(son[u], sp);
for(int i = head[u]; i != -1; i = edge[i].go) {
int v = edge[i].to;
if(v != son[u] && v != fa[u])
get_pos(v, v);
}
}

// 树状数组部分
int c[MAX_N];
int n;

inline int lowbit(int x) {
return x & (-x);
}

int get_sum(int i) {
int s = 0;
while(i > 0) {
s += c[i];
i -= lowbit(i);
}
return s;
}

void add(int i, int val) {
while(i <= n) {
c[i] += val;
i += lowbit(i);
}
}

// u->v路径上的点值改变为val
void update(int u, int v, int val) {
int f1 = top[u], f2 = top[v];
while(f1 != f2) {
if(deep[f1] < deep[f2]) {
swap(f1, f2);
swap(u, v);
}
add(p[f1], val);
add(p[u]+1, -val);
u = fa[f1];
f1 = top[u];
}
if(deep[u] > deep[v]) swap(u, v);
add(p[u], val);
add(p[v]+1, -val);
}

int a[MAX_N];

int main() {
int m, q;
while(~ scanf("%d %d %d", &n, &m, &q)) {
int u, v, k;
char op[10];
init();
for(int i = 1; i <= n; ++i)
scanf("%d", &a[i]);
while(m--) {
scanf("%d %d", &u, &v);
add_edge(u, v);
add_edge(v, u);
}
dfs(1, 0, 0);
get_pos(1, 1);
memset(c, 0, sizeof(c));
for(int i = 1; i <= n; ++i) {
add(p[i], a[i]);
add(p[i]+1, -a[i]);
}
while(q--) {
scanf("%s", op);
if(op[0] == 'Q') {
scanf("%d", &u);
printf("%d\n", get_sum(p[u]));
} else {
scanf("%d %d %d", &u, &v, &k);
if(op[0] == 'D') k = -k;
update(u, v, k);
}
}
}
return 0;
}
  • 边权

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
// 修改单条边,查询路径边权最大值 SPOJ QTREE
#define mm (l + r) >> 1
#define lson l, mid, p << 1
#define rson mid + 1, r, p << 1|1
const int MAX_N = 1e5 + 10;
struct node {
int to, go;
}edge[MAX_N << 1];
int head[MAX_N], tot, pos;
// tot是边表下标,pos是剖分后点个数(1开始)
int top[MAX_N]; // top[v]表示v所在重链的顶端节点
int fa[MAX_N]; // 父亲节点
int deep[MAX_N]; // 深度
int num[MAX_N]; // num[v]表示v为根的子树节点数量
int p[MAX_N]; // p[v]表示v对应的位置
int fp[MAX_N]; // fp和p数组相反
int son[MAX_N]; // 重儿子

void init() {
tot = 0;
memset(head, -1, sizeof(head));
pos = 1;
memset(son, -1, sizeof(son));
}

void add_edge(int u, int v) {
edge[tot].to = v, edge[tot].go = head[u], head[u] = tot++;
}

//第一遍 dfs 求出fa, deep, num, son
void dfs(int u, int pre, int dep) {
deep[u] = dep;
fa[u] = pre;
num[u] = 1;
for(int i = head[u]; i != -1; i = edge[i].go) {
int v = edge[i].to;
if(v != pre) {
dfs(v, u, dep+1);
num[u] += num[v];
if(son[u] == -1 || num[v] > num[son[u]])
son[u] = v;
}
}
}

//第二遍 dfs 求出top, p
void get_pos(int u, int sp) {
top[u] = sp;
p[u] = pos++;
fp[p[u]] = u;
if(son[u] == -1) return ;
get_pos(son[u], sp);
for(int i = head[u]; i != -1; i = edge[i].go) {
int v = edge[i].to;
if(v != son[u] && v != fa[u])
get_pos(v, v);
}
}

// 线段树部分
int tree[MAX_N << 2];

inline void push_up(int p) {
tree[p] = max(tree[p << 1], tree[p << 1|1]);
}

void build(int l, int r, int p) {
if(l == r) {
tree[p] = 0;
return ;
}
int mid = mm;
build(lson), build(rson);
}

void update(int l, int r, int p, int k, int val) {
if(l == r && r == k) {
tree[p] = val;
return ;
}
int mid = mm;
if(k <= mid) update(lson, k, val);
else update(rson, k, val);
push_up(p);
}

int query(int l, int r, int p, int L, int R) {
if(L <= l && r <= R) return tree[p];
int mid = mm, ret = INT_MIN;
if(L <= mid) ret = max(ret, query(lson, L, R));
if(R > mid) ret = max(ret, query(rson, L, R));
return ret;
}

int find(int u, int v) {
// 查询u->v上最大边权
int f1 = top[u], f2 = top[v];
int tmp = 0;
while(f1 != f2) {
if(deep[f1] < deep[f2]) {
swap(f1, f2);
swap(u, v);
}
tmp = max(tmp, query(1, pos, 1, p[f1], p[u]));
u = fa[f1], f1 = top[u];
}
if(u == v) return tmp;
if(deep[u] > deep[v]) swap(u, v);
return max(tmp, query(1, pos, 1, p[ son[u] ], p[v]));
}

int G[MAX_N][3];

int main() {
int T, n, u, v;
char op[10];
scanf("%d", &T);
while(T --) {
init();
scanf("%d", &n);
for(int i = 1; i < n; ++i) {
scanf("%d %d %d", &G[i][0], &G[i][1], &G[i][2]);
add_edge(G[i][0], G[i][1]);
add_edge(G[i][1], G[i][0]);
}
dfs(1, 0, 0);
get_pos(1, 1);
build(1, pos, 1);
for(int i = 1; i < n; ++i) {
if(deep[ G[i][0] ] > deep[ G[i][1] ])
swap(G[i][0], G[i][1]);
update(1, pos, 1, p[ G[i][1] ], G[i][2]);
}
while(scanf("%s", op)) {
if(op[0] == 'D') break ;
scanf("%d %d", &u, &v);
if(op[0] == 'Q')
printf("%d\n", find(u, v));
else update(1, pos, 1, p[ G[u][1] ], v);
}
}
return 0;
}
  • 并查集

  • 普通并查集

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
const int MAX_N = 1e5 + 10;
int f[MAX_N];
void init(int n) {
for(int i = 0; i <= n; ++i)
f[i] = i;
}
int find(int x) {
return f[x] == x ? x : f[x] = find(f[x]);
}
bool merge(int x, int y) {
x = find(x), y = find(y);
if(x == y) return false ;
f[y] = x;
return true; // false表示已经在同一集合
}
  • 带权并查集

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
const int MAX_N = 1e5 + 10;
const int MOD = 3;
// MOD实际上指的是关系数
int f[MAX_N], p[MAX_N];
// p数组存权值

void init(int n) {
for(int i = 0; i <= n; ++i)
f[i] = i, p[i] = 0;
}

int find(int x) {
int tmp = f[x];
if(f[x] != x) {
f[x] = find(f[x]); // 压缩路径
p[x] = (p[x] + p[tmp]) % MOD;
// 合并权值
}
return f[x];
}

void merge(int x, int y, int s) {
// 合并集合并加权,这里一般要s-1
int fx = find(x), fy = find(y);
if(fx != fy) {
f[fy] = fx;
p[fy] = (p[x] - p[y] + s + MOD) % MOD;
}
}

bool solve(int x, int y, int ask, int n) {
// 判断x、y权值关系,是否满足要求
int fx = find(x), fy = find(y);
bool flag = true;
if(x > n || y > n || (ask == 1 && x == y))
flag = false ;
else if(fx == fy)
flag = ((p[y] - p[x] + MOD) % MOD == ask);
return flag;
}

数论部分

  • 一些基础的公式、定理

  1. 组合数公式:Cnm=Cn1m1+Cn1mC_n^m = C_{n-1}^{m-1} + C_{n-1}^m,即杨辉三角,可用在O(n2)O(n^2)下求出所有组合数
  2. Lucas定理:Cnm=Cn/pm/pCn%pm%pmodpC_n^m = C_{n/p}^{m/p} * C_{n\%p}^{m\%p} \mod p,要求pp是质数
  3. 错排公式:an=(n1)an1an2a_n = (n-1) * a_{n-1} * a_{n-2}
  4. 欧拉定理:若gcd(a,n)=1gcd(a,n) = 1,则有aϕ(n)1(modn)a^{\phi(n)} \equiv 1 \pmod n,其中ϕ(n)\phi(n)是欧拉函数。
  5. 费马小定理:若gcd(a,p)=1gcd(a, p) = 1,则有ap11(modp)a^{p-1} \equiv 1 \pmod p
  6. 费马大定理:当n>2n > 2时,xn+yn=znx^n + y^n = z^n没有正整数解
  7. 威尔逊定理:当且仅当pp是素数时,有(p1)!1(modp)(p-1)! \equiv 1 \pmod p
  8. 快速幂:对于任意正整数a,na,n,满足an={(an/2)2n是偶数a(an/2)2n是奇数a^n = \begin{cases} {(a^{n/2})}^2 &\text{n是偶数} \\ a * (a^{n/2})^2 &\text{n是奇数} \end{cases}
  9. 平方和公式:x=1nx2=(2n+1)(n+1)n/6\sum_{x=1}^n x^2 = (2*n + 1) * (n + 1) * n / 6
  10. 等比数列求和:Sn=a1(1qn)/(1q)S_n = a_1 * (1-q^n) / (1-q)
  11. 素数间隔定理:对于大质数pnp_npn1p_{n-1}之间间隔不超过log2pn\log_2p_n
  12. 待更新……
  • 因数、余数、质数基础

  • 埃氏素数筛法

1
2
3
4
5
6
7
8
9
10
11
12
13
const int MAX_N = 1e5 + 10;
int CNT, prime[MAX_N];
bool vis[MAX_N];
void make_prime() {
CNT = 0;
for(int i = 2; i < MAX_N; ++i) {
if(!vis[i]) {
prime[++CNT] = i;
for(int j = 2; j * i < MAX_N; ++j)
vis[i*j] = true;
}
}
}
  • 欧拉函数

对于任意正整数nn,欧拉函数是[1,n][1,n]中与nn互质的数的个数,通常写作φ(n)\varphi(n)。特别地,φ(1)=1\varphi(1) = 1

通式为:φ(x)=xi=1n(11pi)\varphi(x) = x \displaystyle\prod_{i=1}^{n} (1 - \frac{1}{p_i}),其中p1,p2,,pnp_1, p_2, \dots, p_nnn的所有质因子。

  1. 求解单个数字其欧拉函数:
1
2
3
4
5
6
7
8
9
10
11
int phi(int n) {
int ans = n;
for(int i = 2; i * i <= n; ++i) {
if(n % i == 0) {
ans -= ans / i;
while(n % i == 0) n /= i;
}
}
if(n > 1) ans -= ans / n;
return ans;
}
  1. 打表求欧拉函数:
1
2
3
4
5
6
7
8
9
10
11
12
const int MAX_N = 1e5 + 10;
int phi[MAX_N];
void euler() {
for(int i = 2; i < MAX_N; ++i) {
if(!phi[i]) {
for(int j = i; j < MAX_N; j += i) {
if(!phi[j]) phi[j] = j;
phi[j] = phi[j] / i * (i-1);
}
}
}
}
  1. 线性筛素数、欧拉函数:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
const int MAX_N = 1e6 + 10;
int CNT, prime[MAX_N], phi[MAX_N];
bool vis[MAX_N];
void make_prime_phi() {
CNT = 0;
memset(vis, false, sizeof(vis));
for(int i = 2; i < MAX_N; ++i) {
if(!vis[i]) {
prime[++CNT] = i;
phi[i] = i-1;
}
for(int j = 1; j <= CNT && i*prime[j] < MAX_N; ++j) {
vis[ i*prime[j] ] = true;
if(i % prime[j] == 0) {
phi[ i*prime[j] ] = phi[i] * prime[j];
break ;
} else phi[ i*prime[j] ] = phi[i] * (prime[j]-1);
}
}
}

欧拉函数的一些性质:

  1. gcd(n,m)=1gcd(n, m) = 1,则φ(nm)=φ(n)φ(m)\varphi(nm) = \varphi(n) \varphi(m)
  2. n>2n>2时,所有的φ(n)\varphi(n)都是偶数
  3. nn为奇数时,有φ(2n)=φ(n)\varphi(2n) = \varphi(n)
  4. gcd(i,n)=1ni=nφ(n)/2\displaystyle\sum_{gcd(i, n) = 1}^{n}{i} = n * \varphi(n) / 2
  5. dnφ(d)=n\displaystyle\sum_{d|n} \varphi(d) = n
  • 欧几里得定理求a、b最大公因数、最小公倍数

1
2
3
4
5
6
int gcd(int a, int b) {
return b ? gcd(b, a % b) : a;
}
int lcm(int a, int b) {
return a / gcd(a, b) * b;
}
  • 扩展欧几里得

求一组整数特解(x,y)(x, y),满足ax+by=gcd(a,b)ax + by = gcd(a, b),其中aabb已知。

易知上述方程可转化为bx+(a%b)y=gcd(b,a%b)=gcd(a,b)bx' + (a \% b)y' = gcd(b, a \% b) = gcd(a, b)

则有ax+by=bx+(a%b)yax + by = bx' + (a \% b)y',整理后得

ax+by=ay+b(xaby)ax + by = ay' + b(x' - \lfloor{\frac{a}{b}}\rfloor y')

得出解:x=y,y=xabyx = y', y = x' - \lfloor{\frac{a}{b}}\rfloor y',其中当递归到b=0b = 0时,可得一组特解x=1,y=0x = 1, y = 0

1
2
3
4
5
6
7
8
9
10
int ex_gcd(int a, int b, int &x, int &y) {
if(!b) {
x = 1, y = 0;
return a;
}
int d = ex_gcd(b, a%b, x, y);
int t = x;
x = y, y = t - a / b * y;
return d;
}

得到一组特解(x1,y1)(x_1, y_1)后,方程ax+by=gcd(a,b)ax + by = gcd(a, b)的通解为:x=x1+k(b/d),y=y1+k(a/d)x = x_1 + k * (b/d), y = y_1 + k * (a / d),其中kk为任意整数。ax+by=cgcd(a,b)ax + by = c*gcd(a, b)cc为任意正整数)的通解为:x=cx1+k(b/d),y=cy1+k(a/d)x = c * x_1 + k * (b/d), y = c * y_1 + k * (a / d)

  • 线性同余方程

形如:axb(modp)ax \equiv b \pmod p的方程,称为线性同余方程,有如下两个定理:

  1. d=gcd(a,p)d = gcd(a, p),当dbd|b时,方程恰有且好有dd个模nn不同余的解,否则方程无解。
  2. x0x_0是方程任一解,则该方程通解为xi=(x0+k(p/d))%nx_i = (x_0 + k * (p / d)) \% n,其中k=0,1,2,,d1k = 0, 1, 2, \dots, d-1

注意:扩展欧几里得求出axpy=gcd(a,p)ax-py=gcd(a,p)的一组解后,并不一定是原同余方程的解。

特别地,设t=n/dt = n / d,若求出一特解xx,则该方程的最小正整数解为x%tx \% t,这里的模运算右对齐。

  • 模线性方程组求解(中国剩余定理)

对于模线性方程组如:

{xa1(modm1)xa2(modm2)xan(modmn)\begin{cases} {x \equiv a_1 \pmod{m_1}} \\ {x \equiv a_2 \pmod{m_2}} \\ \dots \\ {x \equiv a_n \pmod{m_n}}\end{cases}

假设整数m1,m2,,mnm_1, m_2, \dots, m_n两两互质,通解可以构造如下:

M=i=1nmiM = \displaystyle\prod_{i=1}^{n} m_i,并设Mi=M/miM_i = M / m_iti=Mi1t_i = M_i^{-1},其中Mi1M_i^{-1}表示MiM_imim_i意义下的倒数,即Miti1(modmi)M_it_i \equiv 1 \pmod{m_i}

方程组通解为:x=i=1naitiMi+kMx = \displaystyle\sum_{i = 1}^{n} {a_it_iM_i + kM}kk为任意整数。

求最小正整数解:

1
2
3
4
5
6
7
8
9
10
11
int crt(int r[], int m[], int n) {
int M = 1, ans = 0, x, y, Mi;
for(int i = 1; i <= n; ++i) M *= m[i];
for(int i = 1; i <= n; ++i) {
Mi = M / m[i];
ex_gcd(Mi, m[i], x, y);
ans = (ans + r[i] * x * Mi) % M;
}
while(ans < 0) ans += M;
return ans;
}

若数m1,m2,,mnm_1, m_2, \dots, m_n不一定互质,求最小正整数解:

1
2
3
4
5
6
7
8
9
10
11
12
int ex_crt(int m[], int r[], int n) {
int M = m[1], R = r[1], d, x, y;
for(int i = 2; i <= n; ++i) {
d = ex_gcd(M, m[i], x, y);
if((R - r[i]) % d) return -1;
x = (R - r[i]) / d * x % m[i];
R -= M * x, M = M / d * m[i];
R %= M;
}
R = (R % M + M) % M;
return R;
}
  • 乘法逆元

在可模运算的数据集合内,定义任意一个aa,若存在一个ee,满足:

a×ee×aa(modp)a \times e \equiv e \times a \equiv a \pmod p

则称ee为该种数据的在模pp意义下的幺元。若存在一个bb,满足:

a×bb×ae(modp)a \times b \equiv b \times a \equiv e \pmod p

则称bbaa在模pp意义下的逆元,通常写作inv(a)inv(a)a1a^{-1}

由上述可知,在正整数集内,11就是一个幺元(pp可取任意值)。

  1. 扩展欧几里得求逆元:

即求同余方程ax1(modp)ax \equiv 1 \pmod p的最小正整数解,前提是gcd(a,p)=1gcd(a, p) = 1

1
2
3
4
5
6
const int MOD = 1e9 + 7;
int inv(int a) {
int x, y;
ex_gcd(a, MOD, x, y);
return (x % MOD + MOD) % MOD;
}
  1. 费马小定理求逆元:

根据ap1modp=1a^{p-1} \mod p = 1a1ap2(modp)a^{-1} \equiv a^{p-2} \pmod p,前提是pp为质数且aa不是pp的倍数。

1
2
3
4
5
6
7
8
9
10
11
12
13
const int MOD = 1e9 + 7;
int qmod(int n, int a) {
int tmp = 1;
if(n) {
tmp = qmod(n >> 1, a);
tmp = tmp * tmp % MOD;
if(n & 1) tmp = tmp * a % MOD;
}
return tmp;
}
int inv(int a) {
return qmod(MOD - 2, a);
}
  1. 线性递推逆元:

前提是模数为一个奇质数。

1
2
3
4
5
6
7
8
const int MAX_N = 1e5 + 10;
const int MOD = 1e9 + 7;
int inv[MAX_N]; // 最大可求到min(MAX_N, MOD)
void make_inv() {
inv[1] = 1;
for(int i = 2; i < MOD && i < MAX_N; ++i)
inv[i] = (MOD - MOD / i) * inv[MOD % i] % MOD;
}
  • 二次剩余

x2n(modp)x^2 \equiv n \pmod p有解,则称nn是模pp的二次剩余。常用于模意义下的开方,即nx(modp)\sqrt{n} \equiv x \pmod p

引入一个概念,勒让德符号:

(ap)={0,a0(modp)1,x2a(modp)1,x2≢a(modp)\left(\frac{a}{p}\right) = \begin{cases} 0, &{a \equiv 0 \pmod p} \\ 1, &{\exists x^2 \equiv a \pmod p} \\ -1, & {\forall x^2 \not\equiv a \pmod p} \end{cases}

又有欧拉准则(前提为pp是奇质数):

(ap)=±1ap12(modp)\left(\frac{a}{p}\right) = \pm1 \equiv a^{\frac{p-1}{2}} \pmod p

可以利用上述欧拉准则进行求解:

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
#define ll long long

struct Point {
// 二次域 x + y*sqrt(w)
ll x, y, w;
Point(int a = 0, int b = 0, int c = 0) {
x = a, y = b, w = c;
}
};

Point mul(Point a, Point b, ll MOD) {
// 二次域乘法
Point ans;
ans.x = (a.x * b.x % MOD + a.y * b.y % MOD * a.w % MOD) % MOD;
ans.y = (a.x * b.y % MOD + a.y * b.x % MOD) % MOD;
ans.w = a.w;
return ans;
}

Point qmod(ll n, Point a, ll MOD) {
Point ans = Point(1, 0, a.w);
for(; n; n >>= 1, a = mul(a, a, MOD))
if(n & 1) ans = mul(ans, a, MOD);
return ans;
}

ll qmod(ll n, ll a, ll MOD) {
ll ans = 1LL;
for(; n; n >>= 1, a = a * a % MOD)
if(n & 1) ans = ans * a % MOD;
return ans;
}

ll lgd(ll a, ll p) {
// 勒让德符号
ll ans = qmod((p - 1) >> 1, a, p);
if(ans + 1 == p) return -1;
return ans;
}

ll solve(ll b, ll p) {
// 求x^2 = b(mod p)
if(!b || p == 2) return b;
if(lgd(b, p) == -1) return -1;
ll w, a;
do {
a = rand() % p;
w = a * a - b;
w = (w % p + p) % p;
} while(lgd(w, p) != -1);
Point ans, tmp = Point(a, 1, w);
ans = qmod((p + 1) >> 1, tmp, p);
return ans.x;
}
  • 普通分解质因数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
const int MAX_N = 100 + 10;
int factor[MAX_N][2], fatcnt;
int get_factors(int x) {
// 需要先make_prime
fatcnt = 0;
for(int i = 1; prime[i] * prime[i] <= x; ++i) {
factor[fatcnt][1] = 0;
if(x % prime[i] == 0) {
factor[++fatcnt][0] = prime[i];
while(x % prime[i] == 0) {
++factor[fatcnt][1];
x /= prime[i];
}
}
}
if(x > 1) {
factor[++fatcnt][0] = x;
factor[fatcnt][1] = 1;
}
return fatcnt;
}
  • 约数的和取余

ABA^B所有因子的和取模MOD,这里有一个计算1+p+p2++pn1 + p + p^2 + \dots + p^n的算法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
ll qmod(ll a, ll n, ll MOD) {
ll ans = 1;
for(; n; n >>= 1, a = a*a % MOD)
if(n & 1) ans = ans*a % MOD;
return ans;
}

ll sum(ll p, ll n, ll MOD) {
if(p == 0) return 0;
if(n == 0) return 1;
if(n & 1)
return (1 + qmod(p, n / 2 + 1, MOD)) % MOD * sum(p, n / 2, MOD) % MOD;
return ((1 + qmod(p, n / 2 + 1, MOD)) % MOD * sum(p, n / 2, MOD) % MOD + qmod(p, n / 2, MOD) % MOD) % MOD;
}

ll solve(ll A, ll B, ll MOD) {
get_factors(A);
ll ans = 1;
for(int i = 1; i <= fatcnt; ++i)
ans = ans * sum(factor[i][0], factor[i][1]*B, MOD) % MOD;
return ans;
}
  • 组合数相关

  • 一些定理

  1. m=0nCnm=2n\displaystyle\sum_{m=0}^{n} C_n^m = 2^n
  2. 奇数项和 = 偶数项和 = 2n12^{n-1}
  3. kCnk=nCn1k1k * C_n^k = n * C_{n-1}^{k-1}
  • 预处理组合数

1
2
3
4
5
6
7
8
9
10
const int MAX_N = 1010;
const int MOD = 1e9 + 7;
int C[MAX_N][MAX_N];
void make_C() {
for(int i = 0; i < MAX_N; ++i)
C[i][0] = C[i][i] = 1;
for(int i = 2; i < MAX_N; ++i)
for(int j = 1; j < i; ++j)
C[i][j] = (C[i - 1][j - 1] + C[i-1][j]) % MOD;
}
  • 求单个组合数

1
2
3
4
5
6
7
8
9
10
const int MAX_N = 1e5 + 10;
const int MOD = 1e9 + 7;
int fac[MAX_N], inv[MAX_N];
// fac是阶乘,inv是阶乘逆元
void cal(int n, int m) {
// 需要make_inv
if(n < m) return 0;
if(n - m < m) m = n - m;
return fac[n] * inv[m] % MOD * inv[ n-m ] % MOD;
}
  • 卢卡斯定理

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#define ll long long
ll inv(ll a, ll p) {
// a^(p-2) % p
return qmod(p-2, a, p);
}

ll cal(ll a, ll b, ll p) {
if(a < b) return 0;
if(b > a-b) b = a-b;
ll fz = 1, fm = 1;
for(ll i = 0; i < b; ++i) {
fz = fz * (a-i) % p;
fm = fm * (i+1) % p;
}
return fz * inv(fm, p) % p;
}

ll lucas(ll a, ll b, ll p) {
if(b == 0 || a == b) return 1;
return cal(a % p, b % p, p) * lucas(a / p, b / p, p) % p;
}
  • 斯特林数

  1. 第一类斯特林数:形如[nm]\left[\begin{matrix} {n} \\ {m} \end{matrix}\right]或者s(n,m)s(n, m),表示将nn个不同元素分成mm组,每组是一个的方案数。由于是一个环,所以(1,2,3,4,5)(1, 2, 3, 4, 5)(5,4,3,2,1)(5, 4, 3, 2, 1)其实是同一种方案。其中又分无符号斯特林数sus_{u}和有符号斯特林数sss_{s}
    第一类斯特林数又有如下性质:
    • su(n,0)=0ns_{u}(n, 0) = 0^n
    • su(n,n)=su(n,1)=1s_{u}(n, n) = s_{u}(n, 1) = 1
    • su(n,n1)=Cn2s_{u}(n, n-1) = C_n^2
    • su(n,2)=(n1)i=1n11is_{u}(n, 2) = (n-1) \cdot \displaystyle\sum_{i=1}^{n-1} {\frac{1}{i}}
    • su(n,n2)=2Cn3+3Cn4s_{u}(n, n-2) = 2 \cdot C_n^3 + 3 \cdot C_n^4
    • k=0nsu(n,k)=n!\displaystyle\sum_{k=0}^{n}{s_{u}(n, k)} = n!
    • su(n+1,m)=su(n,m1)+nsu(n,m)s_{u}(n+1, m) = s_{u}(n, m-1) + n \cdot s_{u}(n, m)
    • ss(n,m)=(1)n+msu(n,m)s_{s}(n, m) = {(-1)}^{n+m}s_{u}(n, m)
    • k=0nss(n,k)=0n\displaystyle\sum_{k=0}^{n}{s_{s}(n, k)} = 0^n,注意00=10^0 = 1
    • 代码实现第一类斯特林数:
1
2
3
4
5
6
7
8
9
10
11
#define ll long long
const int MAX_N = 2000 + 10;
const int MOD = 1e9 + 7;
ll S[MAX_N][MAX_N];
void make_Stirling_1() {
memset(S, 0LL, sizeof(S));
S[0][0] = S[1][1] = 1LL;
for(int i = 2; i < MAX_N; ++i)
for(int j = 1; j <= i; ++j)
S[i][j] = (S[i-1][j-1] + (i-1) * S[i-1][j] % MOD) % MOD;
}
  1. 第二类斯特林数:形如{nm}\left\{\begin{matrix} {n} \\ {m}\end{matrix}\right\}或者S(n,m)S(n, m),表示将nn个不同元素拆分成mm个集合的方案数。
    则有如下性质:
    • 前三个性质与第一类相同
    • S(n,2)=2n11S(n, 2) = 2^{n-1} - 1
    • S(n,n1)=Cn2S(n, n-1) = C_n^2
    • S(n,n2)=Cn3+3Cn4S(n, n-2) = C_n^3 + 3 \cdot C_n^4
    • S(n,n3)=Cn4+10Cn5+15Cn6S(n, n-3) = C_n^4 + 10 \cdot C_n^5 + 15 \cdot C_n^6
    • S(n,3)=12(3n1+1)2n1S(n, 3) = \frac{1}{2}(3^{n-1} + 1) - 2^{n-1}
    • k=0nS(n,k)=Bn\displaystyle\sum_{k=0}^{n}{S(n, k)} = B_nBnB_n指的是Bell数
    • S(n,m)=1m!k=0m(1)kCmk(mk)nS(n, m) = \frac{1}{m!} \displaystyle_{k=0}^{m} {(-1)}^{k}C_m^k {(m-k)}^n
    • S(n+1,m)=S(n,m1)+mS(n,m)S(n+1, m) = S(n, m-1) + m \cdot S(n, m)
    • 代码实现第二类斯特林数:
1
2
3
4
5
6
7
8
9
10
11
#define ll long long
const int MAX_N = 2000 + 10;
const int MOD = 1e9 + 7;
ll S[MAX_N][MAX_N];
void make_Stirling_2() {
memset(S, 0LL, sizeof(S));
S[0][0] = S[1][1] = 1LL;
for(int i = 2; i < MAX_N; ++i)
for(int j = 1; j <= i; ++j)
S[i][j] = (S[i-1][j-1] + j * S[i-1][j] % MOD) % MOD;
}
  • Miller_Rabin 快速大素数测试

可以测试2632^{63}(甚至更高)以下的数是否是素数(但可能是伪素数)

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
#define ull unsigned long long
const int S = 4; // 测试次数,一般4~8次

ull mulit_mod(ull n, ull a, ull MOD) {
a %= MOD, n %= MOD;
ull ans = 0;
for(; n; n >>= 1, a = (a + a) % MOD)
if(n & 1) ans = (ans + a) % MOD;
return ans;
}

ull pow_mod(ull n, ull a, ull MOD) {
ull ans = 1;
for(; n; n >>= 1, a = mulit_mod(a, a, MOD))
if(n & 1) ans = mulit_mod(ans, a, MOD);
return ans;
}

bool Miller_Rabin(ull n, ull a) {
ull d = n-1, s = 0, t;
while(!(d&1)) d >>= 1, ++s;
t = pow_mod(d, a, n);
if(t == 1 || t == -1) return true;
for(int i = 0; i < s; ++i) {
if(t == n-1) return true;
t = mulit_mod(t, t, n);
}
return false;
}

bool is_prime(ull n) {
if(n < 2) return false;
if(n == 2) return true;
if((n&1) == 0) return false;
srand(time(NULL));
for(int i = 0; i < S; ++i) {
ull a = rand() % (n-1) + 1;
if(a != 1 && n % a == 0) return false;
if(!Miller_Rabin(n, a))
return false;
}
return true;
}
  • Pollard_rho大数分解质因数

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
ll factor[100];
int fatcnt; //质因子个数
ll gcd(ll a, ll b) {
ll t;
while(b)
t = a, a = b, b = t%b;
return a < 0 ? -a : a;
}

ll Pollard_rho(ll x, ll c) {
ll i = 1, k = 2;
srand(time(NULL));
ll x0 = rand() % x;
ll y = x0;
while(true) {
++i;
x0 = (mulit_mod(x0, x0, x) + c) % x;
ll d = gcd(y - x0, x);
if(d != 1 && d != x) return d;
if(y == x0) return x;
if(i == k) {y = x0, k += k;}
}
}
// k一般设置107左右
void find_fac(ll n, ll k) {
if(n == 1) return ;
if(is_prime(n)) { // 调用Miller_Rabin
factor[++fatcnt] = n;
return ;
}
ll p = n;
int c = k;
while(p >= n) p = Pollard_rho(p, c--);
find_fac(p, k);
find_fac(n / p, k);
}
  • 欧拉降幂公式

ab{ab%φ(p),gcd(a,p)=1ab,gcd(a,p)1,b<φ(p)ab%φ(p)+φ(p),gcd(a,p)1,bφ(p)(modp)a^b \equiv \begin{cases} a^{b \% \varphi(p)}, &{\gcd(a, p) = 1} \\ a^b, &{\gcd(a, p) \neq 1, b < \varphi(p)} \\ a^{b \% \varphi(p) + \varphi(p)}, &{\gcd(a, p) \neq 1, b \geq \varphi(p)} \end{cases} \pmod p

利用欧拉降幂求

aaaanmodp{\underbrace{a^{a^{a^{\sdot^{\sdot^{\sdot^{a}}}}}}}_n} \mod p

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
#define ll long long
//需要make_phi
ll qmod(ll a, ll n, ll MOD)
{
ll ans = 1;
for(; n; n >>= 1, a = a*a % MOD)
if(n & 1) ans = ans*a % MOD;
return ans;
}

ll get_log(ll x, ll a) {
if(x < 1) return -1;
return 1 + get_log(log(x)/log(a), a);
}

ll gcd(ll a, ll b) {
return b ? gcd(b, a%b) : a;
}

ll f(ll a, ll n, ll p) {
if(a == 1 || n == 0) return 1%p;
if(p == 1) return 0;
if(n == 1) return a % p;
if(gcd(a, p) == 1) {
return qmod(a, f(a, n-1, phi[p]), p);
} else {
ll k = get_log(phi[p], a);
if(k <= n-1)
return qmod(a, f(a, n-1, phi[p]) + phi[p], p);
else return qmod(a, f(a, n-1, phi[p]), p);
}
}
  • BSGS大步小步算法

Baby-Step Giant-Step算法求解axb(modn)a^x \equiv b \pmod n,求一个x[0,n)x \in {[0, n)},对于a,b,na, b, n没有特殊要求。

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
#define ll long long
const int MOD = 76543; // 用于Hash
int head[MOD], top;
struct node {
int hs, id, go;
node() {}
node(int a, int b, int c) {
hs = a, id = b, go = c;
}
}a[MOD];

void insert(int x, int y) {
int k = x%MOD;
a[++top] = node(x, y, head[k]);
head[k] = top;
}

int find(int x) {
int k = x%MOD;
for(int i = head[k]; i != -1; i = a[i].go)
if(a[i].hs == x)
return a[i].id;
return -1;
}

int BSGS(int a, int b, int n) {
// a^x = b (mod n)
memset(head, -1, sizeof(head));
top = 0;
if(b == 1) return 0;
int m = sqrt(1.0 * n), j;
ll x = 1, p = 1;
for(int i = 0; i < m; ++i, p = p*a % n)
insert(p*b % n, i);
for(ll i = m; ; i += m) {
if((j = find(x = x*p % n)) != -1) return i-j;
if(i > n) break ;
}
return -1;
}
  • 莫比乌斯函数、反演

莫比乌斯函数:

μ(d)={1,d=0(1)k,d=p1p2pk,pipj0,others\mu(d) = \begin{cases} 1, & {d = 0} \\ {(-1)}^k, & {d = p_1 \cdot p_2 \dots p_k}, \forall p_i \not = p_j \\ 0, & {others} \end{cases}

有一些这样的性质:

dnμ(d)={1,n=10,n>0\displaystyle\sum_{d|n} {\mu(d)} = \begin{cases} 1, & n = 1 \\ 0, & n > 0 \end{cases}

dnμ(d)d=φ(n)n,nN+\displaystyle\sum_{d|n} {\frac{\mu(d)}{d}} = {\frac{\varphi(n)}{n}} , {n \in N_+}

根据这样的性质便可以线性筛莫比乌斯函数:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
const int MAX_N = 1e6 + 10;
bool vis[MAX_N];
int CNT, prime[MAX_N], mu[MAX_N];

void make_mobius() {
memset(vis, false, sizeof(vis));
mu[1] = 1;
CNT = 0;
for(int i = 2; i < MAX_N; ++i) {
if(!vis[i]) {
prime[++CNT] = i;
mu[i] = -1;
}
for(int j = 1; j <= CNT && i*prime[j] < MAX_N; ++j) {
vis[i*prime[j]] = true;
if(i % prime[j] == 0) {
mu[i * prime[j]] = 0;
break ;
} else mu[i * prime[j]] = -mu[i];
}
}
}

莫比乌斯反演:

形式一:

F(n)=dnf(d)    f(n)=dnμ(d)F(nd)F(n) = \displaystyle\sum_{d|n} {f(d)} \iff f(n) = \displaystyle\sum_{d|n} {\mu(d)F(\frac{n}{d})}

形式二:

F(n)=ndf(d)    f(n)=ndμ(dn)F(d)F(n) = \displaystyle\sum_{n|d} {f(d)} \iff f(n) = \displaystyle\sum_{n|d} {\mu(\frac{d}{n})F(d)}

  • 牛顿迭代法求平方根

    1. C++ 版本,精度有限
1
2
3
4
5
6
7
8
9
10
double sqrt_newton(double n) {
const double eps = 1E-15;
double x = 1;
while (1) {
double nx = (x + n / x) / 2;
if (abs(x - nx) < eps) break;
x = nx;
}
return x;
}
  1. JAVA版本,精度较高
1
2
3
4
5
6
7
8
9
10
11
12
public static BigInteger isqrtNewton(BigInteger n) {
BigInteger a = BigInteger.ONE.shiftLeft(n.bitLength() / 2);
boolean p_dec = false;
for(;;) {
BigInteger b = n.divide(a).add(a).shiftRight(1);
if(a.compareTo(b) == 0 || a.compareTo(b) < 0 && p_dec)
break;
p_dec = a.compareTo(b) > 0;
a = b;
}
return a;
}
  • 高斯消元

1
// 待更新
  • 线性基

  • 斐波那契数列

    1. 通项公式:

Fn=(1+52)n(152)n5F_n = \frac{ {\left({ \frac{1 + \sqrt{5}}{2} } \right)}^n - {\left({ \frac{1 - \sqrt{5}}{2} } \right)}^n} {\sqrt{5}}

1
// 待更新

图论部分

  • 邻接表存图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#define ll long long
// 使用前需要将tot清零,head初始化为-1
const ll INF = 1LL << 60;
const int MAX_N = 1e5 + 10; // 点数
const int MAX_M = 1e6 + 10; // 边数

struct Edge {
int to, dis, go;
Edge(int a = 0, int b = 0, int c = 0) {
to = a, dis = b, go = c;
}
}edge[MAX_M];
int head[MAX_N], tot;
void add_edge(int u, int v, int dis) {
edge[++tot] = Edge(v, dis, head[u]), head[u] = tot;
}
  • Vector伪邻接表存图

1
2
3
4
5
6
7
8
9
10
11
12
const int MAX_N = 1e5 + 10;
struct node {
int to, dis;
node() {}
node(int a = 0, int b = 0) {
to = a, dis = b;
}
};
vector<node> G[MAX_N];
void add_edge(int u, int v, int dis) {
G[u].push_back(node(v, dis));
}
  • Dijkstra单源最短路

时间复杂度O(nlog2n)O(n \cdot \log_2n),不可以有负权边

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
//需要邻接表存图
#define ll long long
struct node {
int u; ll dis;
node(int a = 0, ll b = 0) {
u = a, dis = b;
}
};
bool operator < (node a, node b) {
// 这里注意是<右边的元素会放在前面
if(a.dis != b.dis) return a.dis > b.dis;
else return a.u > b.u;
}

bool vis[MAX_N];
ll dist[MAX_N];

void Dijkstra(int s, int n) {
for(int i = 1; i <= n; ++i)
dist[i] = INF, vis[i] = false;
priority_queue<node> q;
dist[s] = 0;
q.push(node(s, 0));
node t;
int u, v, i;
ll tmp;
while(!q.empty()) {
t = q.top();
q.pop();
u = t.u;
if(vis[u]) continue ;
vis[u] = true;
for(i = head[u]; i != -1; i = edge[i].go) {
v = edge[i].to;
tmp = dist[u] + edge[i].dis;
if(dist[v] > tmp) {
dist[v] = tmp;
q.push(node(v, dist[v]));
}
}
}
}
  • SPFA单源最短路

时间复杂度O(kE)O(kE)kk是点入队次数,EE是边数,其时间复杂度不固定,最坏可达到O(nE)O(nE),但是可以判断是否有负权环(某个点入队nn次),也可以解决负权边问题。

有两种优化:用双端队列、栈代替普通队列

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
//需要邻接表存图
#define ll long long
bool vis[MAX_N];
ll dist[MAX_N];

void SPFA(int s, int n) {
for(int i = 1; i <= n; ++i)
dist[i] = INF, vis[i] = false;
queue<int> q;
dist[s] = 0;
q.push(s);
vis[s] = true;
int u, v, i;
ll tmp;
while(!q.empty()) {
u = q.front();
q.pop();
vis[u] = false;
for(i = head[u]; i != -1; i = edge[i].go) {
v = edge[i].to;
tmp = dist[u] + edge[i].dis;
if(dist[v] > tmp) {
dist[v] = tmp;
if(!vis[v]) {
vis[v] = true;
q.push(v);
}
}
}
}
}
  • Floyd多源最短路

时间复杂度O(n3)O(n^3)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
const int MAX_N = 110; // 点数
const int INF = 0x3f3f3f3f;
int G[MAX_N][MAX_N], dist[MAX_N][MAX_N];

void floyd(int n) {
for(int i = 1; i <= n; ++i)
for(int j = 1; j <= n; ++j)
dist[i][j] = INF;
for(int k = 1; k <= n; ++k)
for(int i = 1; i <= n; ++i)
for(int j = 1; j <= n; ++j) {
if(G[i][k] && G[k][j]) // 表示两点连通
dist[i][j] = min(dist[i][j], G[i][k] + G[k][j]);
if(G[i][j]) dist[i][j] = min(dist[i][j], G[i][j]);
}
}
  • 最小生成树MST

  • Kruskal

利用并查集思想,用于解决大部分无向图的最小生成树问题。

时间复杂度大概在O(mlog2m+nα(n)))O(m * log_2m + n * \alpha(n)))

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
//需要邻接表
bool cmp(Edge x, Edge y) {
if(x.dis != y.dis)
return x.dis < y.dis;
if(x.from != y.from)
return x.from < y.from;
return x.to < y.to;
}

int f[MAX_N];
int find(int x) {
return x == f[x] ? x : f[x] = find(f[x]);
}
void merge(int x, int y) {
x = find(x), y = find(y);
if(x != y) f[y] = x;
}
bool is_merge(int x, int y) {
x = find(x), y = find(y);
return x == y;
}

int kruskal(int n, int m) {
int k = 0, ans = 0;
for(int i = 1; i <= n; ++i) f[i] = i;
sort(edge + 1, edge + 1 + m, cmp);
for(int i = 1; i <= m && k < n-1; ++i) {
int u = edge[i].from, v = edge[i].to;
if(!is_merge(u, v)) {
merge(u, v);
ans += edge[i].dis;
++k;
}
}
return k == n-1 ? ans : -1;
}
  • Prim

将点集分为已访问点VV和未访问点EE,初始全是未访问点,先将任一点分到已访问点。选择VVEE最短的一条边,将EE中该点加入VV,直到所有点都在VV

用堆优化后时间复杂度为O((n+m)log2n)O((n+m) * log_2n),看起来比Kruskal优秀,但实际跑起来不一定。

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
// 需要邻接表
struct node {
int u, dis;
node(int a = 0, int b = 0) {
u = a, dis = b;
}
};
bool operator< (node x, node y) {
if(x.dis != y.dis) return x.dis > y.dis;
if(x.u != y.u) return x.u > y.u;
}
bool vis[MAX_N];
int dist[MAX_N];

void init(int n) {
tot = 0;
for(int i = 0; i <= n; ++i)
vis[i] = false, head[i] = -1, dist[i] = inf;
}

int prim(int n) {
// 默认选择点1为初始点
int i, u, v, k = 0, ans = 0;
node t;
priority_queue<node> q;
dist[1] = 0;
q.push(node(1, 0));
//vis[1] = true;
while(k < n && !q.empty()) {
t = q.top(); q.pop();
u = t.u;
if(vis[u]) continue ;
vis[u] = true, ans += t.dis, ++k;
for(i = head[u]; i != -1; i = edge[i].go) {
v = edge[i].to;
if(!vis[v] && dist[v] > edge[i].dis) {
dist[v] = edge[i].dis;
q.push(node(v, dist[v]));
}
}
}
return k == n ? ans : -1;
}
  • 最近公共祖先(LCA)

1
// 待更新
  • 二分图匹配

1
// 待更新
  • 网络流

1
// 待更新
  • 树分治(点分治)

1
// 待更新

字符串处理部分

  • KMP

1
// 待更新
  • exKMP

1
// 待更新

博弈论部分

  • nim博弈

nn堆物品,每堆aia_i个,两个玩家轮流取走任意堆的任意个物品,不可不取,取走最后一个的获胜。

定义NimSum=i=1naiNimSum = \displaystyle\bigoplus_{i=1}^{n} a_i,当且仅当NimSumNimSum00时,先手必败。

  • Bash博弈

只有一堆nn个物品,两个人轮流从中取物,最多取mm个,不可不取,最后取光者为胜。

当且仅当nmod(m+1)=0n \mod (m+1) = 0时,先手必败。

  • Wythoff博弈

有两堆各若干个物品,两人轮流从中取物,至少取一堆中若干件物品,或从两堆中同时取相同件物品,最后取完者胜利。

设两堆物品分别为x,y(x<y)x,y (x < y)个,令z=yxz = y-x,当且仅当5+12z=x\left\lfloor\frac{\sqrt{5} + 1}{2} \cdot z \right\rfloor = x时,先手必败。

  • 斐波那契博弈

有一堆若干个物品,两人轮流取,先手第一次取若干个,但不能把物品取完,之后每次去的物品数不能超过上一次取的物品数的两倍且至少一件,取走最后一件物品的人获胜。

当且仅当物品总数是斐波那契数时,先手必胜。

  • 环形博弈

若干个物品连成环,每次取11个或是相邻的22个,取到最后一个的胜利。

当且仅当物品数量小于等于22时,先手必胜。

  • SG函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
//f[N]:可改变当前状态的方式,N为方式的种类,f[N]要在getSG之前先预处理
//SG[]:0~n的SG函数值
//S[]:为x后继状态的集合
int f[N],SG[MAXN],S[MAXN];
void getSG(int n){
int i,j;
memset(SG,0,sizeof(SG));
//因为SG[0]始终等于0,所以i从1开始
for(i = 1; i <= n; i++){
//每一次都要将上一状态 的 后继集合 重置
memset(S,0,sizeof(S));
for(j = 0; f[j] <= i && j <= N; j++)
S[SG[i-f[j]]] = 1; //将后继状态的SG函数值进行标记
for(j = 0;; j++) if(!S[j]){ //查询当前后继状态SG值中最小的非零值
SG[i] = j;
break;
}
}
}

动态规划部分

  • 背包

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
const int MAX_N = 1e5 + 10; // 物品件数
const int MAX_M = 1e5 + 10; // 最大cost
int N, M, dp[MAX_M]; // 件数、cost、答案

void ZeroOnePack(int cost, int weight) {
// 0-1背包,对于一件物品的代价cost,价值weight
for(int i = M; i >= cost; --i)
dp[i] = max(dp[i], dp[i-cost] + weight);
}

void CompletePack(int cost, int weight) {
// 完全背包
for(int i = cost; i <= M; ++i)
dp[i] = max(dp[i], dp[i-cost] + weight);
}

void MultiplePack(int cost, int weight, int cnt) {
// 多重背包,将件数拆成二进制表示
if(cost * cnt >= M) CompletePack(cost, weight);
else {
for(int k = 1; k < cnt; cnt -= k, k <<= 1)
ZeroOnePack(k * cost, k * weight);
if(cnt > 0) ZeroOnePack(cnt * cost, cnt * weight);
}
}
  • 区间DP

针对可以将两个、多个部分整合(分解)得到答案的数据,建立转移方程。

f(i,j)=merge{f(i,k),f(k)}+costf(i, j) = merge\{f(i, k), f(k)\} + cost

costcost指代这段区间合并(分解)所要花费的代价。

计算几何部分

  • Kuangbin几何模板

    1. 二维计算几何
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
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
const double eps = 1e-8;
const double inf = 1e20;
const double pi = acos(-1.0);
const int maxp = 1010;
//`Compares a double to zero`
int sgn(double x){
if(fabs(x) < eps)return 0;
if(x < 0)return -1;
else return 1;
}
//square of a double
inline double sqr(double x){return x*x;}
/*
* Point
* Point() - Empty constructor
* Point(double _x,double _y) - constructor
* input() - double input
* output() - %.2f output
* operator == - compares x and y
* operator < - compares first by x, then by y
* operator - - return new Point after subtracting curresponging x and y
* operator ^ - cross product of 2d points
* operator * - dot product
* len() - gives length from origin
* len2() - gives square of length from origin
* distance(Point p) - gives distance from p
* operator + Point b - returns new Point after adding curresponging x and y
* operator * double k - returns new Point after multiplieing x and y by k
* operator / double k - returns new Point after divideing x and y by k
* rad(Point a,Point b)- returns the angle of Point a and Point b from this Point
* trunc(double r) - return Point that if truncated the distance from center to r
* rotleft() - returns 90 degree ccw rotated point
* rotright() - returns 90 degree cw rotated point
* rotate(Point p,double angle) - returns Point after rotateing the Point centering at p by angle radian ccw
*/
struct Point{
double x,y;
Point(){}
Point(double _x,double _y){
x = _x;
y = _y;
}
void input(){
scanf("%lf%lf",&x,&y);
}
void output(){
printf("%.2f %.2f\n",x,y);
}
bool operator == (Point b)const{
return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;
}
bool operator < (Point b)const{
return sgn(x-b.x)== 0?sgn(y-b.y)<0:x<b.x;
}
Point operator -(const Point &b)const{
return Point(x-b.x,y-b.y);
}
//叉积
double operator ^(const Point &b)const{
return x*b.y - y*b.x;
}
//点积
double operator *(const Point &b)const{
return x*b.x + y*b.y;
}
//返回长度
double len(){
return hypot(x,y);//库函数
}
//返回长度的平方
double len2(){
return x*x + y*y;
}
//返回两点的距离
double distance(Point p){
return hypot(x-p.x,y-p.y);
}
Point operator +(const Point &b)const{
return Point(x+b.x,y+b.y);
}
Point operator *(const double &k)const{
return Point(x*k,y*k);
}
Point operator /(const double &k)const{
return Point(x/k,y/k);
}
//`计算pa 和 pb 的夹角`
//`就是求这个点看a,b 所成的夹角`
//`测试 LightOJ1203`
double rad(Point a,Point b){
Point p = *this;
return fabs(atan2( fabs((a-p)^(b-p)),(a-p)*(b-p) ));
}
//`化为长度为r的向量`
Point trunc(double r){
double l = len();
if(!sgn(l))return *this;
r /= l;
return Point(x*r,y*r);
}
//`逆时针旋转90度`
Point rotleft(){
return Point(-y,x);
}
//`顺时针旋转90度`
Point rotright(){
return Point(y,-x);
}
//`绕着p点逆时针旋转angle`
Point rotate(Point p,double angle){
Point v = (*this) - p;
double c = cos(angle), s = sin(angle);
return Point(p.x + v.x*c - v.y*s,p.y + v.x*s + v.y*c);
}
};
/*
* Stores two points
* Line() - Empty constructor
* Line(Point _s,Point _e) - Line through _s and _e
* operator == - checks if two points are same
* Line(Point p,double angle) - one end p , another end at angle degree
* Line(double a,double b,double c) - Line of equation ax + by + c = 0
* input() - inputs s and e
* adjust() - orders in such a way that s < e
* length() - distance of se
* angle() - return 0 <= angle < pi
* relation(Point p) - 3 if point is on line
* 1 if point on the left of line
* 2 if point on the right of line
* pointonseg(double p) - return true if point on segment
* parallel(Line v) - return true if they are parallel
* segcrossseg(Line v) - returns 0 if does not intersect
* returns 1 if non-standard intersection
* returns 2 if intersects
* linecrossseg(Line v) - line and seg
* linecrossline(Line v) - 0 if parallel
* 1 if coincides
* 2 if intersects
* crosspoint(Line v) - returns intersection point
* dispointtoline(Point p) - distance from point p to the line
* dispointtoseg(Point p) - distance from p to the segment
* dissegtoseg(Line v) - distance of two segment
* lineprog(Point p) - returns projected point p on se line
* symmetrypoint(Point p) - returns reflection point of p over se
*
*/
struct Line{
Point s,e;
Line(){}
Line(Point _s,Point _e){
s = _s;
e = _e;
}
bool operator ==(Line v){
return (s == v.s)&&(e == v.e);
}
//`根据一个点和倾斜角angle确定直线,0<=angle<pi`
Line(Point p,double angle){
s = p;
if(sgn(angle-pi/2) == 0){
e = (s + Point(0,1));
}
else{
e = (s + Point(1,tan(angle)));
}
}
//ax+by+c=0
Line(double a,double b,double c){
if(sgn(a) == 0){
s = Point(0,-c/b);
e = Point(1,-c/b);
}
else if(sgn(b) == 0){
s = Point(-c/a,0);
e = Point(-c/a,1);
}
else{
s = Point(0,-c/b);
e = Point(1,(-c-a)/b);
}
}
void input(){
s.input();
e.input();
}
void adjust(){
if(e < s)swap(s,e);
}
//求线段长度
double length(){
return s.distance(e);
}
//`返回直线倾斜角 0<=angle<pi`
double angle(){
double k = atan2(e.y-s.y,e.x-s.x);
if(sgn(k) < 0)k += pi;
if(sgn(k-pi) == 0)k -= pi;
return k;
}
//`点和直线关系`
//`1 在左侧`
//`2 在右侧`
//`3 在直线上`
int relation(Point p){
int c = sgn((p-s)^(e-s));
if(c < 0)return 1;
else if(c > 0)return 2;
else return 3;
}
// 点在线段上的判断
bool pointonseg(Point p){
return sgn((p-s)^(e-s)) == 0 && sgn((p-s)*(p-e)) <= 0;
}
//`两向量平行(对应直线平行或重合)`
bool parallel(Line v){
return sgn((e-s)^(v.e-v.s)) == 0;
}
//`两线段相交判断`
//`2 规范相交`
//`1 非规范相交`
//`0 不相交`
int segcrossseg(Line v){
int d1 = sgn((e-s)^(v.s-s));
int d2 = sgn((e-s)^(v.e-s));
int d3 = sgn((v.e-v.s)^(s-v.s));
int d4 = sgn((v.e-v.s)^(e-v.s));
if( (d1^d2)==-2 && (d3^d4)==-2 )return 2;
return (d1==0 && sgn((v.s-s)*(v.s-e))<=0) ||
(d2==0 && sgn((v.e-s)*(v.e-e))<=0) ||
(d3==0 && sgn((s-v.s)*(s-v.e))<=0) ||
(d4==0 && sgn((e-v.s)*(e-v.e))<=0);
}
//`直线和线段相交判断`
//`-*this line -v seg`
//`2 规范相交`
//`1 非规范相交`
//`0 不相交`
int linecrossseg(Line v){
int d1 = sgn((e-s)^(v.s-s));
int d2 = sgn((e-s)^(v.e-s));
if((d1^d2)==-2) return 2;
return (d1==0||d2==0);
}
//`两直线关系`
//`0 平行`
//`1 重合`
//`2 相交`
int linecrossline(Line v){
if((*this).parallel(v))
return v.relation(s)==3;
return 2;
}
//`求两直线的交点`
//`要保证两直线不平行或重合`
Point crosspoint(Line v){
double a1 = (v.e-v.s)^(s-v.s);
double a2 = (v.e-v.s)^(e-v.s);
return Point((s.x*a2-e.x*a1)/(a2-a1),(s.y*a2-e.y*a1)/(a2-a1));
}
//点到直线的距离
double dispointtoline(Point p){
return fabs((p-s)^(e-s))/length();
}
//点到线段的距离
double dispointtoseg(Point p){
if(sgn((p-s)*(e-s))<0 || sgn((p-e)*(s-e))<0)
return min(p.distance(s),p.distance(e));
return dispointtoline(p);
}
//`返回线段到线段的距离`
//`前提是两线段不相交,相交距离就是0了`
double dissegtoseg(Line v){
return min(min(dispointtoseg(v.s),dispointtoseg(v.e)),min(v.dispointtoseg(s),v.dispointtoseg(e)));
}
//`返回点p在直线上的投影`
Point lineprog(Point p){
return s + ( ((e-s)*((e-s)*(p-s)))/((e-s).len2()) );
}
//`返回点p关于直线的对称点`
Point symmetrypoint(Point p){
Point q = lineprog(p);
return Point(2*q.x-p.x,2*q.y-p.y);
}
};
//圆
struct circle{
Point p;//圆心
double r;//半径
circle(){}
circle(Point _p,double _r){
p = _p;
r = _r;
}
circle(double x,double y,double _r){
p = Point(x,y);
r = _r;
}
//`三角形的外接圆`
//`需要Point的+ / rotate() 以及Line的crosspoint()`
//`利用两条边的中垂线得到圆心`
//`测试:UVA12304`
circle(Point a,Point b,Point c){
Line u = Line((a+b)/2,((a+b)/2)+((b-a).rotleft()));
Line v = Line((b+c)/2,((b+c)/2)+((c-b).rotleft()));
p = u.crosspoint(v);
r = p.distance(a);
}
//`三角形的内切圆`
//`参数bool t没有作用,只是为了和上面外接圆函数区别`
//`测试:UVA12304`
circle(Point a,Point b,Point c,bool t){
Line u,v;
double m = atan2(b.y-a.y,b.x-a.x), n = atan2(c.y-a.y,c.x-a.x);
u.s = a;
u.e = u.s + Point(cos((n+m)/2),sin((n+m)/2));
v.s = b;
m = atan2(a.y-b.y,a.x-b.x) , n = atan2(c.y-b.y,c.x-b.x);
v.e = v.s + Point(cos((n+m)/2),sin((n+m)/2));
p = u.crosspoint(v);
r = Line(a,b).dispointtoseg(p);
}
//输入
void input(){
p.input();
scanf("%lf",&r);
}
//输出
void output(){
printf("%.2lf %.2lf %.2lf\n",p.x,p.y,r);
}
bool operator == (circle v){
return (p==v.p) && sgn(r-v.r)==0;
}
bool operator < (circle v)const{
return ((p<v.p)||((p==v.p)&&sgn(r-v.r)<0));
}
//面积
double area(){
return pi*r*r;
}
//周长
double circumference(){
return 2*pi*r;
}
//`点和圆的关系`
//`0 圆外`
//`1 圆上`
//`2 圆内`
int relation(Point b){
double dst = b.distance(p);
if(sgn(dst-r) < 0)return 2;
else if(sgn(dst-r)==0)return 1;
return 0;
}
//`线段和圆的关系`
//`比较的是圆心到线段的距离和半径的关系`
int relationseg(Line v){
double dst = v.dispointtoseg(p);
if(sgn(dst-r) < 0)return 2;
else if(sgn(dst-r) == 0)return 1;
return 0;
}
//`直线和圆的关系`
//`比较的是圆心到直线的距离和半径的关系`
int relationline(Line v){
double dst = v.dispointtoline(p);
if(sgn(dst-r) < 0)return 2;
else if(sgn(dst-r) == 0)return 1;
return 0;
}
//`两圆的关系`
//`5 相离`
//`4 外切`
//`3 相交`
//`2 内切`
//`1 内含`
//`需要Point的distance`
//`测试:UVA12304`
int relationcircle(circle v){
double d = p.distance(v.p);
if(sgn(d-r-v.r) > 0)return 5;
if(sgn(d-r-v.r) == 0)return 4;
double l = fabs(r-v.r);
if(sgn(d-r-v.r)<0 && sgn(d-l)>0)return 3;
if(sgn(d-l)==0)return 2;
if(sgn(d-l)<0)return 1;
}
//`求两个圆的交点,返回0表示没有交点,返回1是一个交点,2是两个交点`
//`需要relationcircle`
//`测试:UVA12304`
int pointcrosscircle(circle v,Point &p1,Point &p2){
int rel = relationcircle(v);
if(rel == 1 || rel == 5)return 0;
double d = p.distance(v.p);
double l = (d*d+r*r-v.r*v.r)/(2*d);
double h = sqrt(r*r-l*l);
Point tmp = p + (v.p-p).trunc(l);
p1 = tmp + ((v.p-p).rotleft().trunc(h));
p2 = tmp + ((v.p-p).rotright().trunc(h));
if(rel == 2 || rel == 4)
return 1;
return 2;
}
//`求直线和圆的交点,返回交点个数`
int pointcrossline(Line v,Point &p1,Point &p2){
if(!(*this).relationline(v))return 0;
Point a = v.lineprog(p);
double d = v.dispointtoline(p);
d = sqrt(r*r-d*d);
if(sgn(d) == 0){
p1 = a;
p2 = a;
return 1;
}
p1 = a + (v.e-v.s).trunc(d);
p2 = a - (v.e-v.s).trunc(d);
return 2;
}
//`得到过a,b两点,半径为r1的两个圆`
int gercircle(Point a,Point b,double r1,circle &c1,circle &c2){
circle x(a,r1),y(b,r1);
int t = x.pointcrosscircle(y,c1.p,c2.p);
if(!t)return 0;
c1.r = c2.r = r;
return t;
}
//`得到与直线u相切,过点q,半径为r1的圆`
//`测试:UVA12304`
int getcircle(Line u,Point q,double r1,circle &c1,circle &c2){
double dis = u.dispointtoline(q);
if(sgn(dis-r1*2)>0)return 0;
if(sgn(dis) == 0){
c1.p = q + ((u.e-u.s).rotleft().trunc(r1));
c2.p = q + ((u.e-u.s).rotright().trunc(r1));
c1.r = c2.r = r1;
return 2;
}
Line u1 = Line((u.s + (u.e-u.s).rotleft().trunc(r1)),(u.e + (u.e-u.s).rotleft().trunc(r1)));
Line u2 = Line((u.s + (u.e-u.s).rotright().trunc(r1)),(u.e + (u.e-u.s).rotright().trunc(r1)));
circle cc = circle(q,r1);
Point p1,p2;
if(!cc.pointcrossline(u1,p1,p2))cc.pointcrossline(u2,p1,p2);
c1 = circle(p1,r1);
if(p1 == p2){
c2 = c1;
return 1;
}
c2 = circle(p2,r1);
return 2;
}
//`同时与直线u,v相切,半径为r1的圆`
//`测试:UVA12304`
int getcircle(Line u,Line v,double r1,circle &c1,circle &c2,circle &c3,circle &c4){
if(u.parallel(v))return 0;//两直线平行
Line u1 = Line(u.s + (u.e-u.s).rotleft().trunc(r1),u.e + (u.e-u.s).rotleft().trunc(r1));
Line u2 = Line(u.s + (u.e-u.s).rotright().trunc(r1),u.e + (u.e-u.s).rotright().trunc(r1));
Line v1 = Line(v.s + (v.e-v.s).rotleft().trunc(r1),v.e + (v.e-v.s).rotleft().trunc(r1));
Line v2 = Line(v.s + (v.e-v.s).rotright().trunc(r1),v.e + (v.e-v.s).rotright().trunc(r1));
c1.r = c2.r = c3.r = c4.r = r1;
c1.p = u1.crosspoint(v1);
c2.p = u1.crosspoint(v2);
c3.p = u2.crosspoint(v1);
c4.p = u2.crosspoint(v2);
return 4;
}
//`同时与不相交圆cx,cy相切,半径为r1的圆`
//`测试:UVA12304`
int getcircle(circle cx,circle cy,double r1,circle &c1,circle &c2){
circle x(cx.p,r1+cx.r),y(cy.p,r1+cy.r);
int t = x.pointcrosscircle(y,c1.p,c2.p);
if(!t)return 0;
c1.r = c2.r = r1;
return t;
}

//`过一点作圆的切线(先判断点和圆的关系)`
//`测试:UVA12304`
int tangentline(Point q,Line &u,Line &v){
int x = relation(q);
if(x == 2)return 0;
if(x == 1){
u = Line(q,q + (q-p).rotleft());
v = u;
return 1;
}
double d = p.distance(q);
double l = r*r/d;
double h = sqrt(r*r-l*l);
u = Line(q,p + ((q-p).trunc(l) + (q-p).rotleft().trunc(h)));
v = Line(q,p + ((q-p).trunc(l) + (q-p).rotright().trunc(h)));
return 2;
}
//`求两圆相交的面积`
double areacircle(circle v){
int rel = relationcircle(v);
if(rel >= 4)return 0.0;
if(rel <= 2)return min(area(),v.area());
double d = p.distance(v.p);
double hf = (r+v.r+d)/2.0;
double ss = 2*sqrt(hf*(hf-r)*(hf-v.r)*(hf-d));
double a1 = acos((r*r+d*d-v.r*v.r)/(2.0*r*d));
a1 = a1*r*r;
double a2 = acos((v.r*v.r+d*d-r*r)/(2.0*v.r*d));
a2 = a2*v.r*v.r;
return a1+a2-ss;
}
//`求圆和三角形pab的相交面积`
//`测试:POJ3675 HDU3982 HDU2892`
double areatriangle(Point a,Point b){
if(sgn((p-a)^(p-b)) == 0)return 0.0;
Point q[5];
int len = 0;
q[len++] = a;
Line l(a,b);
Point p1,p2;
if(pointcrossline(l,q[1],q[2])==2){
if(sgn((a-q[1])*(b-q[1]))<0)q[len++] = q[1];
if(sgn((a-q[2])*(b-q[2]))<0)q[len++] = q[2];
}
q[len++] = b;
if(len == 4 && sgn((q[0]-q[1])*(q[2]-q[1]))>0)swap(q[1],q[2]);
double res = 0;
for(int i = 0;i < len-1;i++){
if(relation(q[i])==0||relation(q[i+1])==0){
double arg = p.rad(q[i],q[i+1]);
res += r*r*arg/2.0;
}
else{
res += fabs((q[i]-p)^(q[i+1]-p))/2.0;
}
}
return res;
}
};

/*
* n,p Line l for each side
* input(int _n) - inputs _n size polygon
* add(Point q) - adds a point at end of the list
* getline() - populates line array
* cmp - comparision in convex_hull order
* norm() - sorting in convex_hull order
* getconvex(polygon &convex) - returns convex hull in convex
* Graham(polygon &convex) - returns convex hull in convex
* isconvex() - checks if convex
* relationpoint(Point q) - returns 3 if q is a vertex
* 2 if on a side
* 1 if inside
* 0 if outside
* convexcut(Line u,polygon &po) - left side of u in po
* gercircumference() - returns side length
* getarea() - returns area
* getdir() - returns 0 for cw, 1 for ccw
* getbarycentre() - returns barycenter
*
*/
struct polygon{
int n;
Point p[maxp];
Line l[maxp];
void input(int _n){
n = _n;
for(int i = 0;i < n;i++)
p[i].input();
}
void add(Point q){
p[n++] = q;
}
void getline(){
for(int i = 0;i < n;i++){
l[i] = Line(p[i],p[(i+1)%n]);
}
}
struct cmp{
Point p;
cmp(const Point &p0){p = p0;}
bool operator()(const Point &aa,const Point &bb){
Point a = aa, b = bb;
int d = sgn((a-p)^(b-p));
if(d == 0){
return sgn(a.distance(p)-b.distance(p)) < 0;
}
return d > 0;
}
};
//`进行极角排序`
//`首先需要找到最左下角的点`
//`需要重载号好Point的 < 操作符(min函数要用) `
void norm(){
Point mi = p[0];
for(int i = 1;i < n;i++)mi = min(mi,p[i]);
sort(p,p+n,cmp(mi));
}
//`得到凸包`
//`得到的凸包里面的点编号是0$\sim$n-1的`
//`两种凸包的方法`
//`注意如果有影响,要特判下所有点共点,或者共线的特殊情况`
//`测试 LightOJ1203 LightOJ1239`
void getconvex(polygon &convex){
sort(p,p+n);
convex.n = n;
for(int i = 0;i < min(n,2);i++){
convex.p[i] = p[i];
}
if(convex.n == 2 && (convex.p[0] == convex.p[1]))convex.n--;//特判
if(n <= 2)return;
int &top = convex.n;
top = 1;
for(int i = 2;i < n;i++){
while(top && sgn((convex.p[top]-p[i])^(convex.p[top-1]-p[i])) <= 0)
top--;
convex.p[++top] = p[i];
}
int temp = top;
convex.p[++top] = p[n-2];
for(int i = n-3;i >= 0;i--){
while(top != temp && sgn((convex.p[top]-p[i])^(convex.p[top-1]-p[i])) <= 0)
top--;
convex.p[++top] = p[i];
}
if(convex.n == 2 && (convex.p[0] == convex.p[1]))convex.n--;//特判
convex.norm();//`原来得到的是顺时针的点,排序后逆时针`
}
//`得到凸包的另外一种方法`
//`测试 LightOJ1203 LightOJ1239`
void Graham(polygon &convex){
norm();
int &top = convex.n;
top = 0;
if(n == 1){
top = 1;
convex.p[0] = p[0];
return;
}
if(n == 2){
top = 2;
convex.p[0] = p[0];
convex.p[1] = p[1];
if(convex.p[0] == convex.p[1])top--;
return;
}
convex.p[0] = p[0];
convex.p[1] = p[1];
top = 2;
for(int i = 2;i < n;i++){
while( top > 1 && sgn((convex.p[top-1]-convex.p[top-2])^(p[i]-convex.p[top-2])) <= 0 )
top--;
convex.p[top++] = p[i];
}
if(convex.n == 2 && (convex.p[0] == convex.p[1]))convex.n--;//特判
}
//`判断是不是凸的`
bool isconvex(){
bool s[2];
memset(s,false,sizeof(s));
for(int i = 0;i < n;i++){
int j = (i+1)%n;
int k = (j+1)%n;
s[sgn((p[j]-p[i])^(p[k]-p[i]))+1] = true;
if(s[0] && s[2])return false;
}
return true;
}
//`判断点和任意多边形的关系`
//` 3 点上`
//` 2 边上`
//` 1 内部`
//` 0 外部`
int relationpoint(Point q){
for(int i = 0;i < n;i++){
if(p[i] == q)return 3;
}
getline();
for(int i = 0;i < n;i++){
if(l[i].pointonseg(q))return 2;
}
int cnt = 0;
for(int i = 0;i < n;i++){
int j = (i+1)%n;
int k = sgn((q-p[j])^(p[i]-p[j]));
int u = sgn(p[i].y-q.y);
int v = sgn(p[j].y-q.y);
if(k > 0 && u < 0 && v >= 0)cnt++;
if(k < 0 && v < 0 && u >= 0)cnt--;
}
return cnt != 0;
}
//`直线u切割凸多边形左侧`
//`注意直线方向`
//`测试:HDU3982`
void convexcut(Line u,polygon &po){
int &top = po.n;//注意引用
top = 0;
for(int i = 0;i < n;i++){
int d1 = sgn((u.e-u.s)^(p[i]-u.s));
int d2 = sgn((u.e-u.s)^(p[(i+1)%n]-u.s));
if(d1 >= 0)po.p[top++] = p[i];
if(d1*d2 < 0)po.p[top++] = u.crosspoint(Line(p[i],p[(i+1)%n]));
}
}
//`得到周长`
//`测试 LightOJ1239`
double getcircumference(){
double sum = 0;
for(int i = 0;i < n;i++){
sum += p[i].distance(p[(i+1)%n]);
}
return sum;
}
//`得到面积`
double getarea(){
double sum = 0;
for(int i = 0;i < n;i++){
sum += (p[i]^p[(i+1)%n]);
}
return fabs(sum)/2;
}
//`得到方向`
//` 1 表示逆时针,0表示顺时针`
bool getdir(){
double sum = 0;
for(int i = 0;i < n;i++)
sum += (p[i]^p[(i+1)%n]);
if(sgn(sum) > 0)return 1;
return 0;
}
//`得到重心`
Point getbarycentre(){
Point ret(0,0);
double area = 0;
for(int i = 1;i < n-1;i++){
double tmp = (p[i]-p[0])^(p[i+1]-p[0]);
if(sgn(tmp) == 0)continue;
area += tmp;
ret.x += (p[0].x+p[i].x+p[i+1].x)/3*tmp;
ret.y += (p[0].y+p[i].y+p[i+1].y)/3*tmp;
}
if(sgn(area)) ret = ret/area;
return ret;
}
//`多边形和圆交的面积`
//`测试:POJ3675 HDU3982 HDU2892`
double areacircle(circle c){
double ans = 0;
for(int i = 0;i < n;i++){
int j = (i+1)%n;
if(sgn( (p[j]-c.p)^(p[i]-c.p) ) >= 0)
ans += c.areatriangle(p[i],p[j]);
else ans -= c.areatriangle(p[i],p[j]);
}
return fabs(ans);
}
//`多边形和圆关系`
//` 2 圆完全在多边形内`
//` 1 圆在多边形里面,碰到了多边形边界`
//` 0 其它`
int relationcircle(circle c){
getline();
int x = 2;
if(relationpoint(c.p) != 1)return 0;//圆心不在内部
for(int i = 0;i < n;i++){
if(c.relationseg(l[i])==2)return 0;
if(c.relationseg(l[i])==1)x = 1;
}
return x;
}
};
//`AB X AC`
double cross(Point A,Point B,Point C){
return (B-A)^(C-A);
}
//`AB*AC`
double dot(Point A,Point B,Point C){
return (B-A)*(C-A);
}
//`最小矩形面积覆盖`
//` A 必须是凸包(而且是逆时针顺序)`
//` 测试 UVA 10173`
double minRectangleCover(polygon A){
//`要特判A.n < 3的情况`
if(A.n < 3)return 0.0;
A.p[A.n] = A.p[0];
double ans = -1;
int r = 1, p = 1, q;
for(int i = 0;i < A.n;i++){
//`卡出离边A.p[i] - A.p[i+1]最远的点`
while( sgn( cross(A.p[i],A.p[i+1],A.p[r+1]) - cross(A.p[i],A.p[i+1],A.p[r]) ) >= 0 )
r = (r+1)%A.n;
//`卡出A.p[i] - A.p[i+1]方向上正向n最远的点`
while(sgn( dot(A.p[i],A.p[i+1],A.p[p+1]) - dot(A.p[i],A.p[i+1],A.p[p]) ) >= 0 )
p = (p+1)%A.n;
if(i == 0)q = p;
//`卡出A.p[i] - A.p[i+1]方向上负向最远的点`
while(sgn(dot(A.p[i],A.p[i+1],A.p[q+1]) - dot(A.p[i],A.p[i+1],A.p[q])) <= 0)
q = (q+1)%A.n;
double d = (A.p[i] - A.p[i+1]).len2();
double tmp = cross(A.p[i],A.p[i+1],A.p[r]) *
(dot(A.p[i],A.p[i+1],A.p[p]) - dot(A.p[i],A.p[i+1],A.p[q]))/d;
if(ans < 0 || ans > tmp)ans = tmp;
}
return ans;
}

//`直线切凸多边形`
//`多边形是逆时针的,在q1q2的左侧`
//`测试:HDU3982`
vector<Point> convexCut(const vector<Point> &ps,Point q1,Point q2){
vector<Point>qs;
int n = ps.size();
for(int i = 0;i < n;i++){
Point p1 = ps[i], p2 = ps[(i+1)%n];
int d1 = sgn((q2-q1)^(p1-q1)), d2 = sgn((q2-q1)^(p2-q1));
if(d1 >= 0)
qs.push_back(p1);
if(d1 * d2 < 0)
qs.push_back(Line(p1,p2).crosspoint(Line(q1,q2)));
}
return qs;
}
//`半平面交`
//`测试 POJ3335 POJ1474 POJ1279`
//***************************
struct halfplane:public Line{
double angle;
halfplane(){}
//`表示向量s->e逆时针(左侧)的半平面`
halfplane(Point _s,Point _e){
s = _s;
e = _e;
}
halfplane(Line v){
s = v.s;
e = v.e;
}
void calcangle(){
angle = atan2(e.y-s.y,e.x-s.x);
}
bool operator <(const halfplane &b)const{
return angle < b.angle;
}
};
struct halfplanes{
int n;
halfplane hp[maxp];
Point p[maxp];
int que[maxp];
int st,ed;
void push(halfplane tmp){
hp[n++] = tmp;
}
//去重
void unique(){
int m = 1;
for(int i = 1;i < n;i++){
if(sgn(hp[i].angle-hp[i-1].angle) != 0)
hp[m++] = hp[i];
else if(sgn( (hp[m-1].e-hp[m-1].s)^(hp[i].s-hp[m-1].s) ) > 0)
hp[m-1] = hp[i];
}
n = m;
}
bool halfplaneinsert(){
for(int i = 0;i < n;i++)hp[i].calcangle();
sort(hp,hp+n);
unique();
que[st=0] = 0;
que[ed=1] = 1;
p[1] = hp[0].crosspoint(hp[1]);
for(int i = 2;i < n;i++){
while(st<ed && sgn((hp[i].e-hp[i].s)^(p[ed]-hp[i].s))<0)ed--;
while(st<ed && sgn((hp[i].e-hp[i].s)^(p[st+1]-hp[i].s))<0)st++;
que[++ed] = i;
if(hp[i].parallel(hp[que[ed-1]]))return false;
p[ed]=hp[i].crosspoint(hp[que[ed-1]]);
}
while(st<ed && sgn((hp[que[st]].e-hp[que[st]].s)^(p[ed]-hp[que[st]].s))<0)ed--;
while(st<ed && sgn((hp[que[ed]].e-hp[que[ed]].s)^(p[st+1]-hp[que[ed]].s))<0)st++;
if(st+1>=ed)return false;
return true;
}
//`得到最后半平面交得到的凸多边形`
//`需要先调用halfplaneinsert() 且返回true`
void getconvex(polygon &con){
p[st] = hp[que[st]].crosspoint(hp[que[ed]]);
con.n = ed-st+1;
for(int j = st,i = 0;j <= ed;i++,j++)
con.p[i] = p[j];
}
};
//***************************

const int maxn = 1010;
struct circles{
circle c[maxn];
double ans[maxn];//`ans[i]表示被覆盖了i次的面积`
double pre[maxn];
int n;
circles(){}
void add(circle cc){
c[n++] = cc;
}
//`x包含在y中`
bool inner(circle x,circle y){
if(x.relationcircle(y) != 1)return 0;
return sgn(x.r-y.r)<=0?1:0;
}
//圆的面积并去掉内含的圆
void init_or(){
bool mark[maxn] = {0};
int i,j,k=0;
for(i = 0;i < n;i++){
for(j = 0;j < n;j++)
if(i != j && !mark[j]){
if( (c[i]==c[j])||inner(c[i],c[j]) )break;
}
if(j < n)mark[i] = 1;
}
for(i = 0;i < n;i++)
if(!mark[i])
c[k++] = c[i];
n = k;
}
//`圆的面积交去掉内含的圆`
void init_add(){
int i,j,k;
bool mark[maxn] = {0};
for(i = 0;i < n;i++){
for(j = 0;j < n;j++)
if(i != j && !mark[j]){
if( (c[i]==c[j])||inner(c[j],c[i]) )break;
}
if(j < n)mark[i] = 1;
}
for(i = 0;i < n;i++)
if(!mark[i])
c[k++] = c[i];
n = k;
}
//`半径为r的圆,弧度为th对应的弓形的面积`
double areaarc(double th,double r){
return 0.5*r*r*(th-sin(th));
}
//`测试SPOJVCIRCLES SPOJCIRUT`
//`SPOJVCIRCLES求n个圆并的面积,需要加上init\_or()去掉重复圆(否则WA)`
//`SPOJCIRUT 是求被覆盖k次的面积,不能加init\_or()`
//`对于求覆盖多少次面积的问题,不能解决相同圆,而且不能init\_or()`
//`求多圆面积并,需要init\_or,其中一个目的就是去掉相同圆`
void getarea(){
memset(ans,0,sizeof(ans));
vector<pair<double,int> >v;
for(int i = 0;i < n;i++){
v.clear();
v.push_back(make_pair(-pi,1));
v.push_back(make_pair(pi,-1));
for(int j = 0;j < n;j++)
if(i != j){
Point q = (c[j].p - c[i].p);
double ab = q.len(),ac = c[i].r, bc = c[j].r;
if(sgn(ab+ac-bc)<=0){
v.push_back(make_pair(-pi,1));
v.push_back(make_pair(pi,-1));
continue;
}
if(sgn(ab+bc-ac)<=0)continue;
if(sgn(ab-ac-bc)>0)continue;
double th = atan2(q.y,q.x), fai = acos((ac*ac+ab*ab-bc*bc)/(2.0*ac*ab));
double a0 = th-fai;
if(sgn(a0+pi)<0)a0+=2*pi;
double a1 = th+fai;
if(sgn(a1-pi)>0)a1-=2*pi;
if(sgn(a0-a1)>0){
v.push_back(make_pair(a0,1));
v.push_back(make_pair(pi,-1));
v.push_back(make_pair(-pi,1));
v.push_back(make_pair(a1,-1));
}
else{
v.push_back(make_pair(a0,1));
v.push_back(make_pair(a1,-1));
}
}
sort(v.begin(),v.end());
int cur = 0;
for(int j = 0;j < v.size();j++){
if(cur && sgn(v[j].first-pre[cur])){
ans[cur] += areaarc(v[j].first-pre[cur],c[i].r);
ans[cur] += 0.5*(Point(c[i].p.x+c[i].r*cos(pre[cur]),c[i].p.y+c[i].r*sin(pre[cur]))^Point(c[i].p.x+c[i].r*cos(v[j].first),c[i].p.y+c[i].r*sin(v[j].first)));
}
cur += v[j].second;
pre[cur] = v[j].first;
}
}
for(int i = 1;i < n;i++)
ans[i] -= ans[i+1];
}
};
  1. 三维计算几何
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
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
const double eps = 1e-8;
int sgn(double x){
if(fabs(x) < eps)return 0;
if(x < 0)return -1;
else return 1;
}
struct Point3{
double x,y,z;
Point3(double _x = 0,double _y = 0,double _z = 0){
x = _x;
y = _y;
z = _z;
}
void input(){
scanf("%lf%lf%lf",&x,&y,&z);
}
void output(){
scanf("%.2lf %.2lf %.2lf\n",x,y,z);
}
bool operator ==(const Point3 &b)const{
return sgn(x-b.x) == 0 && sgn(y-b.y) == 0 && sgn(z-b.z) == 0;
}
bool operator <(const Point3 &b)const{
return sgn(x-b.x)==0?(sgn(y-b.y)==0?sgn(z-b.z)<0:y<b.y):x<b.x;
}
double len(){
return sqrt(x*x+y*y+z*z);
}
double len2(){
return x*x+y*y+z*z;
}
double distance(const Point3 &b)const{
return sqrt((x-b.x)*(x-b.x)+(y-b.y)*(y-b.y)+(z-b.z)*(z-b.z));
}
Point3 operator -(const Point3 &b)const{
return Point3(x-b.x,y-b.y,z-b.z);
}
Point3 operator +(const Point3 &b)const{
return Point3(x+b.x,y+b.y,z+b.z);
}
Point3 operator *(const double &k)const{
return Point3(x*k,y*k,z*k);
}
Point3 operator /(const double &k)const{
return Point3(x/k,y/k,z/k);
}
//点乘
double operator *(const Point3 &b)const{
return x*b.x+y*b.y+z*b.z;
}
//叉乘
Point3 operator ^(const Point3 &b)const{
return Point3(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);
}
double rad(Point3 a,Point3 b){
Point3 p = (*this);
return acos( ( (a-p)*(b-p) )/ (a.distance(p)*b.distance(p)) );
}
//变换长度
Point3 trunc(double r){
double l = len();
if(!sgn(l))return *this;
r /= l;
return Point3(x*r,y*r,z*r);
}
};
struct Line3
{
Point3 s,e;
Line3(){}
Line3(Point3 _s,Point3 _e)
{
s = _s;
e = _e;
}
bool operator ==(const Line3 v)
{
return (s==v.s)&&(e==v.e);
}
void input()
{
s.input();
e.input();
}
double length()
{
return s.distance(e);
}
//点到直线距离
double dispointtoline(Point3 p)
{
return ((e-s)^(p-s)).len()/s.distance(e);
}
//点到线段距离
double dispointtoseg(Point3 p)
{
if(sgn((p-s)*(e-s)) < 0 || sgn((p-e)*(s-e)) < 0)
return min(p.distance(s),e.distance(p));
return dispointtoline(p);
}
//`返回点p在直线上的投影`
Point3 lineprog(Point3 p)
{
return s + ( ((e-s)*((e-s)*(p-s)))/((e-s).len2()) );
}
//`p绕此向量逆时针arg角度`
Point3 rotate(Point3 p,double ang)
{
if(sgn(((s-p)^(e-p)).len()) == 0)return p;
Point3 f1 = (e-s)^(p-s);
Point3 f2 = (e-s)^(f1);
double len = ((s-p)^(e-p)).len()/s.distance(e);
f1 = f1.trunc(len); f2 = f2.trunc(len);
Point3 h = p+f2;
Point3 pp = h+f1;
return h + ((p-h)*cos(ang)) + ((pp-h)*sin(ang));
}
//`点在直线上`
bool pointonseg(Point3 p)
{
return sgn( ((s-p)^(e-p)).len() ) == 0 && sgn((s-p)*(e-p)) == 0;
}
};
struct Plane
{
Point3 a,b,c,o;//`平面上的三个点,以及法向量`
Plane(){}
Plane(Point3 _a,Point3 _b,Point3 _c)
{
a = _a;
b = _b;
c = _c;
o = pvec();
}
Point3 pvec()
{
return (b-a)^(c-a);
}
//`ax+by+cz+d = 0`
Plane(double _a,double _b,double _c,double _d)
{
o = Point3(_a,_b,_c);
if(sgn(_a) != 0)
a = Point3((-_d-_c-_b)/_a,1,1);
else if(sgn(_b) != 0)
a = Point3(1,(-_d-_c-_a)/_b,1);
else if(sgn(_c) != 0)
a = Point3(1,1,(-_d-_a-_b)/_c);
}
//`点在平面上的判断`
bool pointonplane(Point3 p)
{
return sgn((p-a)*o) == 0;
}
//`两平面夹角`
double angleplane(Plane f)
{
return acos(o*f.o)/(o.len()*f.o.len());
}
//`平面和直线的交点,返回值是交点个数`
int crossline(Line3 u,Point3 &p)
{
double x = o*(u.e-a);
double y = o*(u.s-a);
double d = x-y;
if(sgn(d) == 0)return 0;
p = ((u.s*x)-(u.e*y))/d;
return 1;
}
//`点到平面最近点(也就是投影)`
Point3 pointtoplane(Point3 p)
{
Line3 u = Line3(p,p+o);
crossline(u,p);
return p;
}
//`平面和平面的交线`
int crossplane(Plane f,Line3 &u)
{
Point3 oo = o^f.o;
Point3 v = o^oo;
double d = fabs(f.o*v);
if(sgn(d) == 0)return 0;
Point3 q = a + (v*(f.o*(f.a-a))/d);
u = Line3(q,q+oo);
return 1;
}
};

杂项

  • 竞赛常用快速读入骚操作

在读入/输出大量(一般10710^7量级以上)的数据,cin/cout会耗费大量的时间去同步scanf/printf的缓冲区,当然,有时候甚至scanf/printf也不够用。

  1. 关闭iostream同步流
1
2
3
4
std::ios::sync_with_stdio(false);
std::cin.tie(0);
// 如果编译开启了 C++11 或更高版本,建议使用 std::cin.tie(nullptr);
// 注意,此后不可和scanf/printf混用
  1. 普通快速读写
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
inline void read_int(int &X) {
X = 0; int w = 0; char ch = 0;
while(!isdigit(ch)) w |= ch=='-', ch = getchar();
while( isdigit(ch)) X = (X<<3)+ (X<<1) + (ch-48), ch = getchar();
X = w ? -X : X;
}

inline void write_int(int x) {
static int sta[35];
int top = 0;
do {
sta[top++] = x % 10, x /= 10;
} while(x);
while(top) putchar(sta[--top] + 48);
}
  1. 利用fread快读
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
namespace fastIO {
#define BUF_SIZE 100000
//fread -> read
bool IOerror = 0;
inline char nc() {
static char buf[BUF_SIZE], *p1 = buf + BUF_SIZE, *pend = buf + BUF_SIZE;
if (p1 == pend) {
p1 = buf;
pend = buf + fread(buf, 1, BUF_SIZE, stdin);
if (pend == p1) {
IOerror = 1;
return -1;
}
}
return *p1++;
}
inline bool blank(char ch) {
return ch == ' ' || ch == '\n' || ch == '\r' || ch == '\t';
}
inline void read(int &x) {
char ch;
while (blank(ch = nc()));
if (IOerror) return;
for (x = ch - '0'; (ch = nc()) >= '0' && ch <= '9'; x = x * 10 + ch - '0');
}
#undef BUF_SIZE
};
using namespace fastIO;
  1. Real FastIO(抄来的)
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
namespace fastIO{
#define BUF_SIZE 100000
#define OUT_SIZE 100000
#define ll long long
//fread->read
bool IOerror=0;
inline char nc(){
static char buf[BUF_SIZE],*p1=buf+BUF_SIZE,*pend=buf+BUF_SIZE;
if (p1==pend){
p1=buf; pend=buf+fread(buf,1,BUF_SIZE,stdin);
if (pend==p1){IOerror=1;return -1;}
//{printf("IO error!\n");system("pause");for (;;);exit(0);}
}
return *p1++;
}
inline bool blank(char ch){return ch==' '||ch=='\n'||ch=='\r'||ch=='\t';}
inline void read(int &x){
bool sign=0; char ch=nc(); x=0;
for (;blank(ch);ch=nc());
if (IOerror)return;
if (ch=='-')sign=1,ch=nc();
for (;ch>='0'&&ch<='9';ch=nc())x=x*10+ch-'0';
if (sign)x=-x;
}
inline void read(ll &x){
bool sign=0; char ch=nc(); x=0;
for (;blank(ch);ch=nc());
if (IOerror)return;
if (ch=='-')sign=1,ch=nc();
for (;ch>='0'&&ch<='9';ch=nc())x=x*10+ch-'0';
if (sign)x=-x;
}
inline void read(double &x){
bool sign=0; char ch=nc(); x=0;
for (;blank(ch);ch=nc());
if (IOerror)return;
if (ch=='-')sign=1,ch=nc();
for (;ch>='0'&&ch<='9';ch=nc())x=x*10+ch-'0';
if (ch=='.'){
double tmp=1; ch=nc();
for (;ch>='0'&&ch<='9';ch=nc())tmp/=10.0,x+=tmp*(ch-'0');
}
if (sign)x=-x;
}
inline void read(char *s){
char ch=nc();
for (;blank(ch);ch=nc());
if (IOerror)return;
for (;!blank(ch)&&!IOerror;ch=nc())*s++=ch;
*s=0;
}
inline void read(char &c){
for (c=nc();blank(c);c=nc());
if (IOerror){c=-1;return;}
}
//fwrite->write
struct Ostream_fwrite{
char *buf,*p1,*pend;
Ostream_fwrite(){buf=new char[BUF_SIZE];p1=buf;pend=buf+BUF_SIZE;}
void out(char ch){
if (p1==pend){
fwrite(buf,1,BUF_SIZE,stdout);p1=buf;
}
*p1++=ch;
}
void print(int x){
static char s[15],*s1;s1=s;
if (!x)*s1++='0';if (x<0)out('-'),x=-x;
while(x)*s1++=x%10+'0',x/=10;
while(s1--!=s)out(*s1);
}
void println(int x){
static char s[15],*s1;s1=s;
if (!x)*s1++='0';if (x<0)out('-'),x=-x;
while(x)*s1++=x%10+'0',x/=10;
while(s1--!=s)out(*s1); out('\n');
}
void print(ll x){
static char s[25],*s1;s1=s;
if (!x)*s1++='0';if (x<0)out('-'),x=-x;
while(x)*s1++=x%10+'0',x/=10;
while(s1--!=s)out(*s1);
}
void println(ll x){
static char s[25],*s1;s1=s;
if (!x)*s1++='0';if (x<0)out('-'),x=-x;
while(x)*s1++=x%10+'0',x/=10;
while(s1--!=s)out(*s1); out('\n');
}
void print(double x,int y){
static ll mul[]={1,10,100,1000,10000,100000,1000000,10000000,100000000,
1000000000,10000000000LL,100000000000LL,1000000000000LL,10000000000000LL,
100000000000000LL,1000000000000000LL,10000000000000000LL,100000000000000000LL};
if (x<-1e-12)out('-'),x=-x;x*=mul[y];
ll x1=(ll)floor(x); if (x-floor(x)>=0.5)++x1;
ll x2=x1/mul[y],x3=x1-x2*mul[y]; print(x2);
if (y>0){out('.'); for (size_t i=1;i<y&&x3*mul[i]<mul[y];out('0'),++i); print(x3);}
}
void println(double x,int y){print(x,y);out('\n');}
void print(char *s){while (*s)out(*s++);}
void println(char *s){while (*s)out(*s++);out('\n');}
void flush(){if (p1!=buf){fwrite(buf,1,p1-buf,stdout);p1=buf;}}
~Ostream_fwrite(){flush();}
}Ostream;
inline void print(int x){Ostream.print(x);}
inline void println(int x){Ostream.println(x);}
inline void print(char x){Ostream.out(x);}
inline void println(char x){Ostream.out(x);Ostream.out('\n');}
inline void print(ll x){Ostream.print(x);}
inline void println(ll x){Ostream.println(x);}
inline void print(double x,int y){Ostream.print(x,y);} //y为小数点后几位
inline void println(double x,int y){Ostream.println(x,y);}
inline void print(char *s){Ostream.print(s);}
inline void println(char *s){Ostream.println(s);}
inline void println(){Ostream.out('\n');}
inline void flush(){Ostream.flush();} //清空
#undef ll
#undef OUT_SIZE
#undef BUF_SIZE
};
  • C++高精度模板

    1. 封装类版
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
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
#define MAXN 9999
// MAXN 是一位中最大的数字
#define MAXSIZE 10024
// MAXSIZE 是位数
#define DLEN 4
// DLEN 记录压几位
struct Big {
int a[MAXSIZE], len;
bool flag; //标记符号'-'
Big() {
len = 1;
memset(a, 0, sizeof a);
flag = 0;
}
Big(const int);
Big(const char*);
Big(const Big&);
Big& operator=(const Big&); //注意这里operator有&,因为赋值有修改……
//由于OI中要求效率
//此处不使用泛型函数
//故不重载
// istream& operator>>(istream&, BigNum&); //重载输入运算符
// ostream& operator<<(ostream&, BigNum&); //重载输出运算符
Big operator+(const Big&) const;
Big operator-(const Big&) const;
Big operator*(const Big&)const;
Big operator/(const int&) const;
// TODO: Big / Big;
Big operator^(const int&) const;
// TODO: Big ^ Big;

// TODO: Big 位运算;

int operator%(const int&) const;
// TODO: Big ^ Big;
bool operator<(const Big&) const;
bool operator<(const int& t) const;
inline void print();
};
// README::不要随随便便把参数都变成引用,那样没办法传值
Big::Big(const int b) {
int c, d = b;
len = 0;
// memset(a,0,sizeof a);
CLR(a);
while (d > MAXN) {
c = d - (d / (MAXN + 1) * (MAXN + 1));
d = d / (MAXN + 1);
a[len++] = c;
}
a[len++] = d;
}
Big::Big(const char* s) {
int t, k, index, l;
CLR(a);
l = strlen(s);
len = l / DLEN;
if (l % DLEN) ++len;
index = 0;
for (int i = l - 1; i >= 0; i -= DLEN) {
t = 0;
k = i - DLEN + 1;
if (k < 0) k = 0;
g(j, k, i) t = t * 10 + s[j] - '0';
a[index++] = t;
}
}
Big::Big(const Big& T) : len(T.len) {
CLR(a);
f(i, 0, len) a[i] = T.a[i];
// TODO:重载此处?
}
Big& Big::operator=(const Big& T) {
CLR(a);
len = T.len;
f(i, 0, len) a[i] = T.a[i];
return *this;
}
Big Big::operator+(const Big& T) const {
Big t(*this);
int big = len;
if (T.len > len) big = T.len;
f(i, 0, big) {
t.a[i] += T.a[i];
if (t.a[i] > MAXN) {
++t.a[i + 1];
t.a[i] -= MAXN + 1;
}
}
if (t.a[big])
t.len = big + 1;
else
t.len = big;
return t;
}
Big Big::operator-(const Big& T) const {
int big;
bool ctf;
Big t1, t2;
if (*this < T) {
t1 = T;
t2 = *this;
ctf = 1;
} else {
t1 = *this;
t2 = T;
ctf = 0;
}
big = t1.len;
int j = 0;
f(i, 0, big) {
if (t1.a[i] < t2.a[i]) {
j = i + 1;
while (t1.a[j] == 0) ++j;
--t1.a[j--];
// WTF?
while (j > i) t1.a[j--] += MAXN;
t1.a[i] += MAXN + 1 - t2.a[i];
} else
t1.a[i] -= t2.a[i];
}
t1.len = big;
while (t1.len > 1 && t1.a[t1.len - 1] == 0) {
--t1.len;
--big;
}
if (ctf) t1.a[big - 1] = -t1.a[big - 1];
return t1;
}
Big Big::operator*(const Big& T) const {
Big res;
int up;
int te, tee;
f(i, 0, len) {
up = 0;
f(j, 0, T.len) {
te = a[i] * T.a[j] + res.a[i + j] + up;
if (te > MAXN) {
tee = te - te / (MAXN + 1) * (MAXN + 1);
up = te / (MAXN + 1);
res.a[i + j] = tee;
} else {
up = 0;
res.a[i + j] = te;
}
}
if (up) res.a[i + T.len] = up;
}
res.len = len + T.len;
while (res.len > 1 && res.a[res.len - 1] == 0) --res.len;
return res;
}
Big Big::operator/(const int& b) const {
Big res;
int down = 0;
gd(i, len - 1, 0) {
res.a[i] = (a[i] + down * (MAXN + 1) / b);
down = a[i] + down * (MAXN + 1) - res.a[i] * b;
}
res.len = len;
while (res.len > 1 && res.a[res.len - 1] == 0) --res.len;
return res;
}
int Big::operator%(const int& b) const {
int d = 0;
gd(i, len - 1, 0) d = (d * (MAXN + 1) % b + a[i]) % b;
return d;
}
Big Big::operator^(const int& n) const {
Big t(n), res(1);
int y = n;
while (y) {
if (y & 1) res = res * t;
t = t * t;
y >>= 1;
}
return res;
}
bool Big::operator<(const Big& T) const {
int ln;
if (len < T.len) return 233;
if (len == T.len) {
ln = len - 1;
while (ln >= 0 && a[ln] == T.a[ln]) --ln;
if (ln >= 0 && a[ln] < T.a[ln]) return 233;
return 0;
}
return 0;
}
inline bool Big::operator<(const int& t) const {
Big tee(t);
return *this < tee;
}
inline void Big::print() {
printf("%d", a[len - 1]);
gd(i, len - 2, 0) { printf("%04d", a[i]); }
}

inline void print(Big s) {
// s不要是引用,要不然你怎么print(a * b);
int len = s.len;
printf("%d", s.a[len - 1]);
gd(i, len - 2, 0) { printf("%04d", s.a[i]); }
}
char s[100024];
  1. 简约版
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
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
static const int LEN = 1004;

int a[LEN], b[LEN], c[LEN], d[LEN];

void clear(int a[LEN]) {
for (int i = 0; i < LEN; ++i) a[i] = 0;
}

void read(int a[LEN]) {
static char s[LEN + 1];
scanf("%s", s);

clear(a);

int len = strlen(s);
for (int i = 0; i < len; ++i) a[len - i - 1] = s[i] - '0';
}

void print(int a[LEN]) {
int i;
for (i = LEN - 1; i >= 1; --i)
if (a[i] != 0) break;
for (; i >= 0; --i) putchar(a[i] + '0');
putchar('\n');
}

void add(int a[LEN], int b[LEN], int c[LEN]) {
clear(c);

for (int i = 0; i < LEN - 1; ++i) {
c[i] += a[i] + b[i];
if (c[i] >= 10) {
c[i + 1] += 1;
c[i] -= 10;
}
}
}

void sub(int a[LEN], int b[LEN], int c[LEN]) {
clear(c);

for (int i = 0; i < LEN - 1; ++i) {
c[i] += a[i] - b[i];
if (c[i] < 0) {
c[i + 1] -= 1;
c[i] += 10;
}
}
}

void mul(int a[LEN], int b[LEN], int c[LEN]) {
clear(c);

for (int i = 0; i < LEN - 1; ++i) {
for (int j = 0; j <= i; ++j) c[i] += a[j] * b[i - j];

if (c[i] >= 10) {
c[i + 1] += c[i] / 10;
c[i] %= 10;
}
}
}

inline bool greater_eq(int a[LEN], int b[LEN], int last_dg, int len) {
if (a[last_dg + len] != 0) return true;
for (int i = len - 1; i >= 0; --i) {
if (a[last_dg + i] > b[i]) return true;
if (a[last_dg + i] < b[i]) return false;
}
return true;
}

void div(int a[LEN], int b[LEN], int c[LEN], int d[LEN]) {
clear(c);
clear(d);

int la, lb;
for (la = LEN - 1; la > 0; --la)
if (a[la - 1] != 0) break;
for (lb = LEN - 1; lb > 0; --lb)
if (b[lb - 1] != 0) break;
if (lb == 0) {
puts("> <");
return;
}

for (int i = 0; i < la; ++i) d[i] = a[i];
for (int i = la - lb; i >= 0; --i) {
while (greater_eq(d, b, i, lb)) {
for (int j = 0; j < lb; ++j) {
d[i + j] -= b[j];
if (d[i + j] < 0) {
d[i + j + 1] -= 1;
d[i + j] += 10;
}
}
c[i] += 1;
}
}
}
  • Java大整数

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
import java.util.*;
import java.math.*;

public class BigInt {
static Scanner in = new Scanner(System.in); // 定义输入对象
public static void main(String[] args) {
BigInteger bigInt_1 = new BigInteger("100");
BigInteger bigInt_2 = BigInteger.valueOf(123);
//两种定义方式,建议使用第一种
bigInt_1.add(bigInt_2); // 加法
bigInt_1.subtract(bigInt_2); // 减法
bigInt_1.multiply(bigInt_2); // 乘法
bigInt_1.divide(bigInt_2);// 除法,向下取整
bigInt_1.divideAndRemainder(bigInt_2);
// 返回一个BigInteger[],包含商和余数
bigInt_1.remainder(bigInt_2); // 取余数,与this同符号
bigInt_1.mod(bigInt_2); // 取模,模数只能为正整数
bigInt_1.pow(10); // a^b
bigInt_1.gcd(bigInt_2); // 最大公约数
bigInt_1.compareTo(bigInt_2);
// 比较大小,<0表示this小,0表示相等,>0表示this大
bigInt_1.equals(bigInt_2); // 真值为相等
bigInt_1.negate(); // -a
bigInt_1.abs(); // 绝对值
bigInt_1.min(bigInt_2); // 最小值
bigInt_1.max(bigInt_2); // 最大值
BigInteger a = in.nextBigInteger(); // 读入
System.out.println(bigInt_1); // 输出
}
}
  • Java大实数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import java.util.*;
import java.math.*;

public class Big {
static Scanner in = new Scanner(System.in); // 定义输入对象
public static void main(String[] args) {
BigDecimal bigDec_1 = new BigDecimal("123.1");
BigDecimal bigDec_2 = BigDecimal.valueOf(233.213);
// 唯有除法与BigInteger不同,无限小数需要规定保留位数
int scale = 3; // 保留位数
bigDec_1.divide(bigDec_2, scale, RoundingMode.HALF_UP);
/*
第三个参数为保留方式,有以下几种:
CEILING 正无穷大方向取整
FLOOR 负无穷大方向取整
DOWN 向 0 的方向取整
UP 正数向正无穷大取整,负数向负无穷大取整
HALF_UP 常用的4舍5入
HALF_DOWN 6,7,8,9 向上取整
HALF_EVEN 小数位是5时,判断整数部分是奇数就进位,6,7,8,9 进位
*/
}
}
  • 手动扩栈

1
#pragma comment(linker, "/STACK:1024000000,1024000000")
  • 对拍

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#include <stdio.h>
#include <stdlib.h>
int main() {
// For Windows
while(true) {
system("data.exe > data.in");
system("std.exe < data.in > std.out");
system("test.exe < data.in > test.out");
if(system("fc std.out test.out")) {
system("pause");
return 0;
}
}
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
#include <stdio.h>
#include <stdlib.h>
int main() {
// For Linux
while(true) {
system("./data > data.in");
system("./std < data.in > std.out");
system("./test < data.in > test.out");
if(system("diff std.out test.out")) {
system("pause"); //方便查看不同处
return 0;
}
}
}
文章作者: Jeson Vendetta Xie
文章链接: http://jvxie.com/2019/08/19/个人模板整理/
版权声明: 本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 JVxie's Blog
打赏
  • 微信
  • 支付宝

评论