算法竞赛板子

算法竞赛板子

目录

1. ST表_区间最值_gcd_按位与_按位或
2. 树状数组
3. 快读
4. 带权并查集
5. 欧拉筛
6. 组合数
7. lucas定理求组合数
8. 离散化
9. 线形基
10. 主席树
11. 约瑟夫环
12. tarjan 求静态LCA
13. tarjan 求无向图割点
14. tarjan 求无向图割点后的连通块
15. tarjan有向图强连通分量缩点
16. exgcd
17. bsgs
18. exbsgs
19. 无向图三元环计数
20. Trie
21. 一些杂项
22. 计算几何

1. ST表_区间最值_gcd_按位与_按位或

struct ST // ST表_区间最值_gcd_按位与_按位或
{int n = 0;vector<int> a;vector<vector<int>> mx, mn;ST(int _n, vector<int> _a){n = _n;a = _a;mx.assign(n + 1, vector<int>(32));mn.assign(n + 1, vector<int>(32, 1e18));for (int j = 0; j < 30; j++) // j 是每一层状态for (int i = 1; i <= n; i++){if (i + (1 << j) - 1 > n)continue;if (!j){mx[i][j] = a[i];mn[i][j] = a[i];}else{mx[i][j] = max(mx[i][j - 1], mx[i + (1 << j - 1)][j - 1]);mn[i][j] = min(mn[i][j - 1], mn[i + (1 << j - 1)][j - 1]);}}}// queryint askmx(int l, int r){assert(l <= r);int k = __lg(r - l + 1);return max(mx[l][k], mx[r + 1 - (1 << k)][k]);}int askmn(int l, int r){assert(l <= r);int k = __lg(r - l + 1);return min(mn[l][k], mn[r + 1 - (1 << k)][k]);}
}; // ST st(n,a);

2. 树状数组

// 树状数组  快速求前缀和 / 前缀最大值
// 维护差分数组(区间加 区间求和) 初始化要add(i,a[i]-a[i-1])
// 位置(统计个数) ...
struct Tree
{int n;vector<int> tr;Tree(int n1){n = n1 + 2;tr.assign(n + 2, 0);}void add(int x, int c) // 加点{for (int i = x; i <= n; i += i & -i)tr[i] += c;}int ask(int x) // 前缀查询{int res = 0;for (int i = x; i; i -= i & -i)res += tr[i];return res;}int ask(int l, int r) // 区间查询{assert(l <= r);if (l > r)return 0ll;return ask(r) - ask(l - 1);}
}; //    Tree tr(n);  tr.add(x,c)

3. 快读

int re()
{int s = 0, f = 1;char ch = getchar();while (ch > '9' || ch < '0'){if (ch == '-')f = -1;ch = getchar();}while ('0' <= ch && ch <= '9')s = s * 10 + ch - 48, ch = getchar();return s * f;
}
void wr(int s)
{if (s < 0)putchar('-'), s = -s;if (s > 9)wr(s / 10);putchar(s % 10 + 48);
}
void wr(int s, char ch) { wr(s), putchar(ch); }

4. 带权并查集

// 带权并查集
struct DSU
{vector<int> p, vs, es; // 集合数 点数 边数 (对一个连通块而言)DSU(int n1)            // p[x]不一定为根节点  find(x)一定是根节点{int n = n1 + 2;p.assign(n, 0);vs.assign(n, 0);es.assign(n, 0);for (int i = 1; i <= n1; i++)p[i] = i, vs[i] = 1, es[i] = 0;}int find(int x) // 找到根节点{if (p[x] == x)return x;int px = find(p[x]);return p[x] = px;}bool same(int a, int b){return find(a) == find(b);}void merge(int a, int b) // 合并集合{int pa = find(a);int pb = find(b);if (pa == pb) // pa pb 均为根节点 p[pa]==pa{es[pa]++; // 1个集合 边+1return;}p[pb] = p[pa];        // 改变b的根节点vs[pa] += vs[pb];     // 将b合并进aes[pa] += es[pb] + 1; // 2个集合}int size(int a) //  集合内的元素的个数{return vs[find(a)];}
};
//  DSU dsu(n);

5. 欧拉筛

        vector<int> is(N), del(N);vector<int> pri;vector<int> yin(N); // 记录最小质因子auto oula = [&](int n){for (int i = 2; i < n; i++){if (!del[i])is[i] = true, pri.push_back(i);for (int j = 0; i * pri[j] <= n; j++){del[i * pri[j]] = true;yin[i * pri[j]] = pri[j]; // 每个数只会被最小质因子删去if (i % pri[j] == 0)break;}}};oula(N);

6. 组合数

constexpr int mod = 998244353;
const int N = 5e5 + 10;
//  组合数
int qp(int a, int b)
{if (!a)return 1ll;int res = 1;while (b){if (b & 1)res = res * a % mod; // 不要取模时记得取消a = a * a % mod;b >>= 1;}return res;
}int fac[N], inv[N]; // fac inv 阶乘与阶乘的逆元
void init()
{fac[0] = inv[0] = 1;for (int i = 1; i < N; ++i)fac[i] = fac[i - 1] * i % mod;inv[N - 1] = qp(fac[N - 1], mod - 2);for (int i = N - 2; i >= 1; --i)inv[i] = inv[i + 1] * (i + 1) % mod;
}
int C(int n, int m)
{return fac[n] * inv[m] % mod * inv[n - m] % mod;
}

7. lucas定理求组合数

int qp(int a, int k, int p) 
{int res = 1 % p;while (k){if (k & 1)res = res * a % p;a = a * a % p;k >>= 1;}return res;
}int C(int a, int b, int p) // 通过定理求组合数C(a, b)
{if (a < b)return 0;int x = 1, y = 1; // x是分子，y是分母for (int i = a, j = 1; j <= b; i--, j++){x = x * i % p;y = y * j % p;}return x * qp(y, p - 2, p) % p;
}int lucas(int a, int b, int p)
{if (a < p && b < p)return C(a, b, p);return C(a % p, b % p, p) * lucas(a / p, b / p, p) % p;
}

8. 离散化

 //  离散化  nums[find(x)] == x  find(x):代替x使用 为x在nums中的指针vector<int> nums;for(int i=0;i<n;i++) {int x;cin>>x;nums.push_back(x);}auto find=[&](int x){return lower_bound(nums.begin(),nums.end(),x)-nums.begin();};sort(nums.begin(),nums.end());nums.erase(unique(nums.begin(),nums.end()),nums.end());

9. 线形基

// 线性基能使用异或运算来表示原数集使用异或运算能表示的所有数。
// 可以极大地缩小异或操作所需的查询次数。
// 1. 判断一个数能否表示成某数集子集的异或和
// 2. 求一个数表示成某数集子集异或和的方案数
// 3. 求某数集子集的最大/最小/第k大/第k小异或和
// 4. 求一个数在某数集子集异或和中的排名
// 5. 所有异或值
// 6. 任意一条 1 到 n 的路径的异或和，都可以由任意一条 1 到 n 路径的异或和与图中的一些环的异或和来组合得到。
struct T
{int bit[64];bool zero = 0;int cnt, top;int d[64];T(){zero = cnt = top = 0;memset(bit, 0, sizeof bit);memset(d, 0, sizeof d);}bool insert(int x){for (int i = 63; i >= 0; i--){if ((x >> i) & 1){if (!bit[i]){bit[i] = x;return 1;}elsex ^= bit[i];}}zero = 1;return false;}int ask(int x) //  询问是否能被异或出{for (int i = 63; i >= 0; i--)if (x >> i & 1)x ^= bit[i];return x == 0;}int askmx() // 异或最大值{int ans = 0;for (int i = 63; i >= 0; i--)if ((ans ^ bit[i]) > ans)ans ^= bit[i];return ans;}//  rebuildvoid rebuild(){cnt = 0;top = 0;for (int i = 63; i >= 0; i--)for (int j = i - 1; j >= 0; j--)if (bit[i] & (1LL << j))bit[i] ^= bit[j];for (int i = 0; i <= 63; i++)if (bit[i])d[cnt++] = bit[i];}int kth(int k) // 异或第k小 需要rebuild{k -= zero;if (!k)return 0;if (k >= (1LL << cnt))return -1;int ans = 0;for (int i = 63; i >= 0; i--)if (k & (1LL << i))ans ^= d[i];return ans;}int rak(int x) // 查询排名{int ans = 0;for (int i = cnt - 1; i >= 0; i--)if (x >= d[i])ans += (1 << i), x ^= d[i];return ans + zero;}
};

10. 主席树

//  主席树  维护离散化后的坐标
//  找区间第k小值
struct Node
{int lc, rc;int cnt;
};
struct Tree
{int n, idx = 0;vector<Node> tr;vector<int> root;Tree(int n1){n = n1 + 2;tr.assign(n * 23, Node{}); // n*4+mlognroot.assign(n, 0);}int build(int l, int r){int p = ++idx;if (l == r)return p;int mid = l + r >> 1;tr[p].lc = build(l, mid);tr[p].rc = build(mid + 1, r);tr[p].cnt = tr[tr[p].lc].cnt + tr[tr[p].rc].cnt;return p;}int insert(int last, int l, int r, int x){int p = ++idx;tr[p] = tr[last];if (l == r){tr[p].cnt++;return p;}int mid = l + r >> 1;if (x <= mid)tr[p].lc = insert(tr[last].lc, l, mid, x);elsetr[p].rc = insert(tr[last].rc, mid + 1, r, x);tr[p].cnt = tr[tr[p].lc].cnt + tr[tr[p].rc].cnt;// pushup(p,tr[p].lc,tr[p].rc);return p;}int ask(int i, int j, int l, int r, int k){if (l == r)return l;int mid = l + r >> 1;int cnt = tr[tr[j].lc].cnt - tr[tr[i].lc].cnt;if (k <= cnt)return ask(tr[i].lc, tr[j].lc, l, mid, k);elsereturn ask(tr[i].rc, tr[j].rc, mid + 1, r, k - cnt);}
}; //    Tree tr(n);void _()
{int n, m;cin >> n >> m;vector<int> a(n + 1);//  离散化  nums[find(x)] == x  find(x):为x在nums中的指针vector<int> nums;for (int i = 1; i <= n; i++){cin >> a[i];nums.push_back(a[i]);}auto find = [&](int x){return lower_bound(nums.begin(), nums.end(), x) - nums.begin();};sort(nums.begin(), nums.end());nums.erase(unique(nums.begin(), nums.end()), nums.end());Tree tr(n); // 版本数tr.root[0] = tr.build(0, nums.size() - 1);for (int i = 1; i <= n; i++)tr.root[i] = tr.insert(tr.root[i - 1], 0, nums.size() - 1, find(a[i]));while (m--){int l, r, k;cin >> l >> r;k = (r - l + 1) / 2 + 1;  // 查询第k小值cout << nums[tr.ask(tr.root[l - 1], tr.root[r], 0, nums.size() - 1, k)] << endl;}
}

11. 约瑟夫环

//  约瑟夫环  0~n-1
//  O(n)
int josephus(int n, int k)
{int res = 0;for (int i = 1; i <= n; ++i)res = (res + k) % i;return res;
}
//  O(klogn)
int josephus(int n, int k)
{if (n == 1)return 0;if (k == 1)return n - 1;if (k > n)return (josephus(n - 1, k) + k) % n; // 线性算法int res = josephus(n - n / k, k);res -= n % k;if (res < 0)res += n; // mod nelseres += res / (k - 1); // 还原位置return res;
}

12. tarjan 求静态LCA

void _() // tarjan 求静态LCA
{int n, q, root;cin >> n >> q >> root;vector<vector<int>> e(n + 1);vector<vector<pair<int, int>>> query(n + 1);vector<int> p(n + 1), vis(n + 1), ans(n + 1);for (int i = 1; i <= n; i++)p[i] = i;auto find = [&](auto &&self, int u) -> int{return p[u] = p[u] == u ? u : self(self, p[u]);};auto tarjan = [&](auto &&self, int u) -> void{vis[u] = 1;for (auto v : e[u]){if (!vis[v]){self(self, v);p[v] = u;}}for (auto [v, idx] : query[u]){if (vis[v])ans[idx] = find(find, v);}};for (int i = 1; i < n; i++){int u, v;cin >> u >> v;e[u].push_back(v);e[v].push_back(u);}for (int i = 1; i <= q; i++){int u, v;cin >> u >> v;query[u].push_back({v, i});query[v].push_back({u, i});}tarjan(tarjan, root);for (int i = 1; i <= q; i++)cout << ans[i] << '\n';
}

13. tarjan 求无向图割点

void _() // tarjan 求无向图割点
{int n, m;cin >> n >> m;vector<vector<int>> e(n + 1);vector<int> dfn(n + 1), low(n + 1), cut(n + 1);int tot = 0, root = 1, cuts = 0;for (int i = 0; i < m; i++){int u, v;cin >> u >> v;e[u].push_back(v);e[v].push_back(u);}auto tarjan = [&](auto &&self, int u) -> void{dfn[u] = low[u] = ++tot;int chd = 0;for (auto v : e[u]){if (!dfn[v]){self(self, v);low[u] = min(low[u], low[v]);if (low[v] >= dfn[u]){chd++;if (u - root || chd > 1)cut[u] = 1;}}elselow[u] = min(low[u], dfn[v]);}};for (int i = 1; i <= n; i++) // 图可不连通if (!dfn[i]){root = i;tarjan(tarjan, i);}for (int i = 1; i <= n; i++)cuts += cut[i];cout << cuts << '\n';for (int i = 1; i <= n; i++)if (cut[i])cout << i << ' ';
}

14. tarjan 求无向图割点后的连通块

void _() // tarjan 求无向图割点后的连通块
{int n, m;cin >> n >> m;vector<vector<int>> e(n + 1);vector<int> dfn(n + 1), low(n + 1), cut(n + 1);int tot = 0, root = 1, cuts = 0;vector<int> res(n + 1), sz(n + 1);for (int i = 0; i < m; i++){int u, v;cin >> u >> v;e[u].push_back(v);e[v].push_back(u);}auto tarjan = [&](auto &&self, int u) -> void{dfn[u] = low[u] = ++tot;int chd = 0;sz[u] = 1;int ss = 1;for (auto v : e[u]){if (!dfn[v]){self(self, v);sz[u] += sz[v]; // dfs搜索树low[u] = min(low[u], low[v]);if (low[v] >= dfn[u]) // 不含返祖边的子树才能作为真子树 其余作为上子树{chd++;ss += sz[v]; // ss 才是割点为子树的res[u] += sz[v] * (n - sz[v]);if (u - root || chd > 1)cut[u] = 1;}}elselow[u] = min(low[u], dfn[v]);}res[u] += ss * (n - ss) + n - 1;if (!cut[u])res[u] = n - 1 << 1;};tarjan(tarjan, root);// for (int i = 1; i <= n; i++)//     bug({i, sz[i]});for (int i = 1; i <= n; i++)cout << res[i] << '\n';
}

15. tarjan有向图强连通分量缩点

void _() // tarjan有向图强连通分量缩点
{int n, m;cin >> n >> m;vector<vector<int>> e(n + 1);vector<pair<int, int>> edges;vector<int> w(n + 1);for (int i = 1; i <= n; i++)cin >> w[i];for (int i = 1; i <= m; i++){int u, v;cin >> u >> v;e[u].push_back(v);edges.push_back({u, v});}int tot = 0;vector<int> dfn(n + 1), low(n + 1), stk(n + 1), instk(n + 1), belong(n + 1);// 时间戳 追溯值 栈维护返祖边和横插边能到达的点（祖先节点和左子树节点）belong属于哪个点代表的强连通分量vector<int> scc(n + 1), sz(n + 1); // 记录点在哪个强连通分量中  强连通分量大小int cnt = 0, top = 0;              // 强连通分量编号auto tarjan = [&](auto &&self, int u) -> void{// 进入u，盖戳、入栈dfn[u] = low[u] = ++tot;stk[++top] = u, instk[u] = 1;for (auto v : e[u]){if (!dfn[v]){self(self, v);low[u] = min(low[u], low[v]);}else if (instk[v]) // 已访问且在栈中low[u] = min(low[u], low[v]);}// 离u，记录SCCif (dfn[u] == low[u]){int v;++cnt;do{v = stk[top--];instk[v] = 0;scc[v] = cnt;++sz[cnt];belong[v] = u;if (u - v)w[u] += w[v];} while (u - v);}};for (int i = 1; i <= n; i++){if (!dfn[i])tarjan(tarjan, i);}// 缩点vector<int> in(n + 1);vector<vector<int>> ne(n + 1);for (auto [u, v] : edges){u = belong[u], v = belong[v];if (u - v){ne[u].push_back(v);in[v]++;}}auto topo = [&](){vector<int> f(n + 1);queue<int> q;for (int i = 1; i <= n; i++){if (belong[i] == i && !in[i]) // dfn[i] == low[i]{q.push(i);f[i] = w[i];}}while (q.size()){auto u = q.front();q.pop();for (auto v : ne[u]){f[v] = max(f[v], w[v] + f[u]);in[v]--;if (!in[v])q.push(v);}}return f;};auto dp = topo();int res = 0;for (int i = 1; i <= n; i++){res = max(res, dp[i]);}cout << res << '\n';
}

16. exgcd

//  exgcd  解同余方程ax+by==c  ax==c(mod b)
int ex_gcd(int a, int b, int &x, int &y)
{if (b == 0){x = 1;y = 0;return a;}int d = ex_gcd(b, a % b, x, y);int temp = x;x = y;y = temp - a / b * y;return d;
}bool liEu(int a, int b, int c, int &x, int &y)
{int d = ex_gcd(a, b, x, y);if (c % d != 0)return false;int k = c / d;x *= k;y *= k;return true;
} // 最小正整数解

17. bsgs

int exgcd(int a, int b, int& x, int& y) {if (!b) {x = 1, y = 0;return a;}int d = exgcd(b, a % b, y, x);y -= a / b * x;return d;
}ll bsgs(int a, int b, int p) {if (1 % p == b % p) return 0;  // 特判t=0,因p可能为1,故写1%pint k = sqrt(p) + 1;  // 分块数// 预处理出所有b * a^y (mod p)umap<int, int> hash;for (int i = 0, j = b % p; i < k; i++) {hash[j] = i;  // 记录值的同时记录对应的y,较大的y会覆盖较小的yj = (ll)j * a % p;}// 预处理出a^kint ak = 1;for (int i = 0; i < k; i++) ak = (ll)ak * a % p;// 遍历x∈[1,k],在哈希表中查找是否存在满足的yfor (int i = 1, j = ak; i <= k; i++) {  // j记录(a^k)^xif (hash.count(j)) return (ll)i * k - hash[j];  // 返回t值j = (ll)j * ak % p;}return -1;  // 无解
}void solve() {int a, b, m, x0, x; cin >> a >> b >> m >> x0 >> x;if (x0 == x) {cout << "YES";return;}int x1 = ((ll)a * x0 + b) % m;if (a == 0) {  // 特判if (x1 == x || b == x) cout << "YES";else cout << "NO";}else if (a == 1) {  // 特判if (!b) cout << (x == x1 ? "YES" : "NO");else {// 求不定方程(n-1)b + mp = t - X0的最小正整数解int x, y;exgcd(b, m, x, y);x = ((ll)x * (x - x1) % m + m) % m;  // 保证x是正数cout << (~x ? "YES" : "NO");}}else {int C = (ll)b * qpow(a - 1, m - 2, m) % m;  // b/(a-1)int A = (x1 + C) % m;  // 分母if (!A) {  // 特判A=0int ans = (-C + m) % m;cout << (ans == x ? "YES" : "NO");}else {int B = (x + C) % m;  // 分子int ans = bsgs(a, (ll)B * qpow(A, m - 2, m) % m, m);cout << (~ans ? "YES" : "NO");}}
}

18. exbsgs

//  bsgs  exbsgs
//  解同余方程 a^x==b(mod p) x的最小解
map<int, int> mp;
int gcd(int a, int b) { return b ? gcd(b, a % b) : a; }
int BSGS(int a, int n, int p, int ad)
{mp.clear();int m = std::ceil(std::sqrt(p));int s = 1;for (int i = 0; i < m; i++, s = 1ll * s * a % p)mp[1ll * s * n % p] = i;for (int i = 0, tmp = s, s = ad; i <= m; i++, s = 1ll * s * tmp % p)if (mp.find(s) != mp.end())if (1ll * i * m - mp[s] >= 0)return 1ll * i * m - mp[s];return -1;
}int exBSGS(int a, int b, int p)
{int n = b;a %= p;n %= p;if (n == 1 || p == 1)return 0;int cnt = 0;int d, ad = 1;while ((d = gcd(a, p)) ^ 1){if (n % d)return -1;cnt++;n /= d;p /= d;ad = (1ll * ad * a / d) % p;if (ad == n)return cnt;}//		std::printf("a: %d n: %d p: %d ad:%d\n",a,n,p,ad);int ans = BSGS(a, n, p, ad);if (ans == -1)return -1;return ans + cnt;
}

19. 无向图三元环计数

// 无向图三元环计数
// 记录每个点的度数  建由度数小到大(相同则按编号)的有向图
// u-->v-->v1 判断u,v1是否连接
void _()
{int n,m;cin>>n>>m;vi d(n+1);map<P,int> mp;For(i,m){int a,b;cin>>a>>b;d[a]++,d[b]++;mp[{a,b}]=1;}vector<int>h(n+1,-1),e(m<<1|1),ne(m<<1|1);int idx=0;auto add=[&](int a,int b){e[idx]=b,ne[idx]=h[a],h[a]=idx++;};for(auto [V,_]:mp){auto [u,v]=V;if(d[u]>d[v]) swap(u,v);if(d[u]==d[v]&&v<u) swap(u,v);add(u,v);}int ans=0;vi f(n+1);//记录For(u,n){for(int i=h[u];i+1;i=ne[i]) {int v=e[i];f[v]=u;}for(int i=h[u];i+1;i=ne[i]){int v=e[i];for(int j=h[v];j+1;j=ne[j]){int v1=e[j];if(f[v1]==u) ans++;}}}cout<<ans<<endl;
}

20. Trie

const int N = 100010;
int idx;        // 各个节点的编号，根节点编号为0
int son[N][26]; // Trie 树本身
// cnt[x] 表示：以 编号为 x 为结尾的字符串的个数
int cnt[N];void insert(string s)
{int p = 0; // 指向根节点for (int i = 0; i < s.size(); i++){// 将当前字符转换成数字（a->0, b->1,...）int u = s[i] - 'a';// 如果数中不能走到当前字符// 为当前字符创建新的节点，保存该字符if (!son[p][u])// 新节点编号为 idx + 1son[p][u] = ++idx;p = son[p][u];}// 这个时候，p 等于字符串 s 的尾字符所对应的 idx// cnt[p] 保存的是字符串 s 出现的次数// 故 cnt[p] ++cnt[p]++;
}int query(string s)
{int p = 0; // 指向根节点for (int i = 0; i < s.size(); i++){// 将当前字符转换成数字（a->0, b->1,...）int u = s[i] - 'a';// 如果走不通了，即树中没有保存当前字符// 则说明树中不存在该字符串if (!son[p][u])return 0;// 指向下一个节点p = son[p][u];}// 循环结束的时候，p 等于字符串 s 的尾字符所对应的 idx//  cnt[p] 就是字符串 s 出现的次数return cnt[p];
}

21. 一些杂项

 l~r中在bit位1的个数auto calc = [&](int N, int bit) // 1~n中在bit位1的个数{++N;return (N >> (bit + 1) << bit) + min(1ll << bit, N % (1ll << (bit + 1)));};auto cal = [&](int l, int r, int bit) // l~r中在bit位1的个数{return l == 1 ? calc(r, bit) : calc(r, bit) - calc(l - 1, bit);};

快速求因数个数for (int i = 1; i < N; i++)for (int j = i; j < N; j += i)cnt[j]++;

 区间异或和
int xor_0_n(int n)
{int res[] = {n, 1, n + 1, 0};return res[n % 4];
}
//  O(1) 区间异或和
int xor_range(int l, int r)
{return xor_0_n(r) ^ xor_0_n(l - 1);
}

随机数rng()
std::mt19937 rng(std::chrono::steady_clock::now().time_since_epoch().count());

22. 计算几何

#include <bits/stdc++.h>
#define int long long
using namespace std;using ld = long double;
#define double long double
const ld eps = 1e-9;
const ld PI = acos(-1);
#define Pt Point<T> // 简写
#define Pd Point<ld>
#define Lt Line<T>
#define Ld Line<ld>int sign(ld x) // 判符号
{if (fabs(x) < eps)return 0;elsereturn x < 0 ? -1 : 1;
}
bool equal(ld x, ld y) { return !sign(x - y); }//  点与向量
template <class T>
struct Point
{T x, y;Point() : x(0), y(0) {}         // 默认构造Point(T x, T y) : x(x), y(y) {} // 带参构造// 运算符重载Point operator+(const Point &b) const { return Point(x + b.x, y + b.y); }Point operator-(const Point &b) const { return Point(x - b.x, y - b.y); }Point operator*(T k) const { return Point(x * k, y * k); }Point operator/(T k) const { return Point(x / k, y / k); }friend bool operator<(Point a, Point b) { return equal(a.x, b.x) ? a.y < b.y - eps : a.x < b.x - eps; }friend bool operator>(Point a, Point b) { return b < a; }friend bool operator==(Point a, Point b) { return !(a < b) && !(b < a); }friend bool operator!=(Point a, Point b) { return a < b || b < a; }friend auto &operator>>(istream &is, Point &p) { return is >> p.x >> p.y; }friend auto &operator<<(ostream &os, Point p) { return os << "(" << p.x << ", " << p.y << ")"; }
};ld len(Point<ld> a) { return sqrtl(a.x * a.x + a.y * a.y); }; // 模长
ld dis(Point<ld> a, Point<ld> b)                              // 两点距离
{ld dx = a.x - b.x;ld dy = a.y - b.y;return sqrtl(dx * dx + dy * dy);
}
ld cross(Point<ld> a, Point<ld> b) { return a.x * b.y - a.y * b.x; } // 叉乘 ab*sin
ld cross(Point<ld> p1, Point<ld> p2, Point<ld> p0) { return cross(p1 - p0, p2 - p0); }
ld dot(Point<ld> a, Point<ld> b) { return a.x * b.x + a.y * b.y; } // 点乘 ab*cos
ld dot(Point<ld> p1, Point<ld> p2, Point<ld> p0) { return dot(p1 - p0, p2 - p0); }
ld angle(Point<ld> a, Point<ld> b) { return abs(atan2(abs(cross(a, b)), a.x * b.x + a.y * b.y)); } // 向量夹角Point<ld> standardize(Point<ld> vec)
{ // 转换为单位向量return vec / sqrtl(vec.x * vec.x + vec.y * vec.y);
}
Point<ld> rotate(Point<ld> p, ld rad) // 向量旋转任意角度
{return {p.x * cos(rad) - p.y * sin(rad), p.x * sin(rad) + p.y * cos(rad)};
}
ld toDeg(ld x) { return x * 180 / PI; }                                              // 弧度转角度
ld toArc(ld x) { return PI / 180 * x; }                                              // 角度转弧度
ld angle(ld a, ld b, ld c) { return acos((a * a + b * b - c * c) / (2.0 * a * b)); } // 余弦定理 三边求角C弧度// 直线
template <class T>
struct Line
{Point<T> a, b;Line(Point<T> a_ = Point<T>(), Point<T> b_ = Point<T>()) : a(a_), b(b_) {}template <class U>operator Line<U>(){ // 自动类型匹配return Line<U>(Point<U>(a), Point<U>(b));}friend auto &operator<<(ostream &os, Line l){return os << "<" << l.a << ", " << l.b << ">";}
};
// 点是否在直线上（三点是否共线）
template <class T>
bool onLine(Point<T> a, Point<T> b, Point<T> c) { return sign(cross(b, a, c)) == 0; }
template <class T>
bool onLine(Point<T> p, Line<T> l) { return onLine(p, l.a, l.b); }template <class T>
ld disPointToLine(Pt p, Lt l) // 点到直线的最近距离
{ld ans = cross(p, l.a, l.b);return abs(ans) / dis(l.a, l.b); // 面积除以底边长
}
template <class T>
bool pointOnLineSide(Pt p1, Pt p2, Lt vec) // 两点是否在直线同侧
{T val = cross(p1, vec.a, vec.b) * cross(p2, vec.a, vec.b);return sign(val) == 1;
}
template <class T>
bool pointNotOnLineSide(Pt p1, Pt p2, Lt vec) // 两点是否在线段异侧
{T val = cross(p1, vec.a, vec.b) * cross(p2, vec.a, vec.b);return sign(val) == -1;
}
Pd lineIntersection(Ld l1, Ld l2) // 两直线相交交点 先判平行
{ld val = cross(l2.b - l2.a, l1.a - l2.a) / cross(l2.b - l2.a, l1.a - l1.b);return l1.a + (l1.b - l1.a) * val;
}
template <class T>
bool lineParallel(Lt p1, Lt p2) { return sign(cross(p1.a - p1.b, p2.a - p2.b)) == 0; } // 两直线是否平行
template <class T>
bool lineVertical(Lt p1, Lt p2) { return sign(dot(p1.a - p1.b, p2.a - p2.b)) == 0; } // 两直线是否垂直
Pd project(Pd p, Ld l)                                                               // 点在直线上的投影点（垂足）
{Pd vec = l.b - l.a;ld r = dot(vec, p - l.a) / (vec.x * vec.x + vec.y * vec.y);return l.a + vec * r;
}
template <class T>
Point<T> rotate(Point<T> a, Point<T> b) // 旋转
{Point<T> vec = a - b;return {-vec.y, vec.x};
}
template <class T>
Lt midSegment(Lt l) // 线段的中垂线
{Pt mid = (l.a + l.b) / 2;return {mid, mid + rotate(l.a, l.b)};
}
ld area(Point<ld> a, Point<ld> b, Point<ld> c) { return abs(cross(b, c, a)) / 2; } // 三角形面积// 凸包
template <class T> // flag 用于判定凸包边上的点、重复的顶点是否要加入到凸包中，为 0 时代表加入凸包；为 1 时不加入凸包。
// 时间复杂度为 O(NlogN)
vector<Point<T>> staticConvexHull(vector<Point<T>> A, int flag = 1)
{int n = A.size();if (n <= 2) // 特判{return A;}vector<Point<T>> ans(n * 2);sort(A.begin(), A.end());int now = -1;for (int i = 0; i < n; i++) // 维护下凸包{while (now > 0 && cross(A[i], ans[now], ans[now - 1]) < flag){now--;}ans[++now] = A[i];}int pre = now;for (int i = n - 2; i >= 0; i--) // 维护上凸包{while (now > pre && cross(A[i], ans[now], ans[now - 1]) < flag){now--;}ans[++now] = A[i];}ans.resize(now);return ans;
}template <class T> // 点与凸包的位置关系 0代表点在凸包外面 1代表在凸壳上 2代表在凸包内部
int contains(Point<T> p, vector<Point<T>> A)
{int n = A.size();bool in = false;for (int i = 0; i < n; i++){Point<T> a = A[i] - p, b = A[(i + 1) % n] - p;if (a.y > b.y){swap(a, b);}if (a.y <= 0 && 0 < b.y && cross(a, b) < 0){in = !in;}if (cross(a, b) == 0 && dot(a, b) <= 0){return 1;}}return in ? 2 : 0;
}
template <class T>
ld area(vector<Point<T>> P) // 任意多边形的面积
{int n = P.size();ld ans = 0;for (int i = 0; i < n; i++){ans += cross(P[i], P[(i + 1) % n]);}return ans / 2.0;
}
template <class T>
vector<Point<T>> mincowski(vector<Point<T>> P1, vector<Point<T>> P2) // 计算两个凸包合成的大凸包。
{int n = P1.size(), m = P2.size();vector<Point<T>> V1(n), V2(m);for (int i = 0; i < n; i++){V1[i] = P1[(i + 1) % n] - P1[i];}for (int i = 0; i < m; i++){V2[i] = P2[(i + 1) % m] - P2[i];}vector<Point<T>> ans = {P1.front() + P2.front()};int t = 0, i = 0, j = 0;while (i < n && j < m){Point<T> val = sign(cross(V1[i], V2[j])) > 0 ? V1[i++] : V2[j++];ans.push_back(ans.back() + val);}while (i < n)ans.push_back(ans.back() + V1[i++]);while (j < m)ans.push_back(ans.back() + V2[j++]);return ans;
}ld RoatingCalipers(vector<Point<ld>> poly) // 旋转卡壳求凸包直径
{if (poly.size() == 1)return 0;if (poly.size() == 2)return len(poly[1] - poly[0]);if (poly.size() == 3)return max({len(poly[1] - poly[0]), len(poly[2] - poly[1]), len(poly[0] - poly[2])});size_t cur = 0;ld ans = 0;for (size_t i = 0; i < poly.size(); i++){size_t j = (i + 1) % poly.size();Line<ld> line(poly[i], poly[j]);while (disPointToLine(poly[cur], line) <= disPointToLine(poly[(cur + 1) % poly.size()], line)){cur = (cur + 1) % poly.size();}ans = max(ans, max(len(poly[i] - poly[cur]), len(poly[j] - poly[cur])));}return ans;
}// 圆
pair<Pd, ld> pointToCircle(Pd p, Pd o, ld r) // 点到圆的最近点
{Pd U = o, V = o;ld d = dis(p, o);if (sign(d) == 0){ // p 为圆心时返回圆心本身return {o, 0};}ld val1 = r * abs(o.x - p.x) / d;ld val2 = r * abs(o.y - p.y) / d * ((o.x - p.x) * (o.y - p.y) < 0 ? -1 : 1);U.x += val1, U.y += val2;V.x -= val1, V.y -= val2;if (dis(U, p) < dis(V, p)){return {U, dis(U, p)};}else{return {V, dis(V, p)};}
}
// 根据圆心角获取圆上某点  将圆上最右侧的点以圆心为旋转中心，逆时针旋转 rad 度
Point<ld> getPoint(Point<ld> p, ld r, ld rad) { return {p.x + cos(rad) * r, p.y + sin(rad) * r}; }// 直线是否与圆相交及交点  同时返回相交状态和交点，分为三种情况：0代表不相交；1代表相切；2代表相交。
tuple<int, Pd, Pd> lineCircleCross(Ld l, Pd o, ld r)
{Pd P = project(o, l);ld d = dis(P, o), tmp = r * r - d * d;if (sign(tmp) == -1){return {0, {}, {}};}else if (sign(tmp) == 0){return {1, P, {}};}Pd vec = standardize(l.b - l.a) * sqrtl(tmp);return {2, P + vec, P - vec};
}
// 两圆是否相交及交点
// 同时返回相交状态和交点，分为四种情况： 0-内含, 1-相离, 2-相切, 3-相交
tuple<int, Pd, Pd> circleIntersection(Pd p1, ld r1, Pd p2, ld r2)
{ld x1 = p1.x, x2 = p2.x, y1 = p1.y, y2 = p2.y, d = dis(p1, p2);if (sign(abs(r1 - r2) - d) == 1){return {0, {}, {}};}else if (sign(r1 + r2 - d) == -1){return {1, {}, {}};}ld a = r1 * (x1 - x2) * 2, b = r1 * (y1 - y2) * 2, c = r2 * r2 - r1 * r1 - d * d;ld p = a * a + b * b, q = -a * c * 2, r = c * c - b * b;ld cosa, sina, cosb, sinb;if (sign(d - (r1 + r2)) == 0 || sign(d - abs(r1 - r2)) == 0){cosa = -q / p / 2;sina = sqrt(1 - cosa * cosa);Point<ld> p0 = {x1 + r1 * cosa, y1 + r1 * sina};if (sign(dis(p0, p2) - r2)){p0.y = y1 - r1 * sina;}return {2, p0, p0};}else{ld delta = sqrt(q * q - p * r * 4);cosa = (delta - q) / p / 2;cosb = (-delta - q) / p / 2;sina = sqrt(1 - cosa * cosa);sinb = sqrt(1 - cosb * cosb);Pd ans1 = {x1 + r1 * cosa, y1 + r1 * sina};Pd ans2 = {x1 + r1 * cosb, y1 + r1 * sinb};if (sign(dis(ans1, p1) - r2))ans1.y = y1 - r1 * sina;if (sign(dis(ans2, p2) - r2))ans2.y = y1 - r1 * sinb;if (ans1 == ans2)ans1.y = y1 - r1 * sina;return {3, ans1, ans2};}
}// 两圆相交面积
ld circleIntersectionArea(Pd p1, ld r1, Pd p2, ld r2)
{ld x1 = p1.x, x2 = p2.x, y1 = p1.y, y2 = p2.y, d = dis(p1, p2);if (sign(abs(r1 - r2) - d) >= 0){return PI * min(r1 * r1, r2 * r2);}else if (sign(r1 + r2 - d) == -1){return 0;}ld theta1 = angle(r1, dis(p1, p2), r2);ld area1 = r1 * r1 * (theta1 - sin(theta1 * 2) / 2);ld theta2 = angle(r2, dis(p1, p2), r1);ld area2 = r2 * r2 * (theta2 - sin(theta2 * 2) / 2);return area1 + area2;
}// 三点确定一圆
tuple<int, Pd, ld> getCircle(Pd A, Pd B, Pd C)
{if (onLine(A, B, C)){ // 特判三点共线return {0, {}, 0};}Ld l1 = midSegment(Line<ld>{A, B});Ld l2 = midSegment(Line<ld>{A, C});Pd O = lineIntersection(l1, l2);return {1, O, dis(A, O)};
}
// --------------------------// void _() // 计算凸包周长
// {
//     int n;
//     cin >> n;
//     vector<Point<ld>> s(n);
//     for (int i = 0; i < n; i++)
//         cin >> s[i];
//     auto vec = staticConvexHull(s);
//     ld res = 0;
//     if (vec.size() > 1)
//     {
//         for (int i = 1; i < vec.size(); i++)
//             res += dis(vec[i - 1], vec[i]);
//         // 加上最后一条边（闭合）
//         res += dis(vec.back(), vec.front());
//     }
//     cout << res << '\n';
// }// void _() // 计算凸包直径
// {
//     int n;
//     cin >> n;
//     vector<Point<ld>> pionts(n);
//     for (auto &v : pionts)
//         cin >> v;
//     auto vec = staticConvexHull(pionts);
//     if (vec.size() < 2)
//     {
//         cout << "0\n";
//         return;
//     }
//     ld res = RoatingCalipers(vec); // 直径
//     cout << (int)round(res * res) << '\n';
// }