結果
| 問題 | No.2829 GCD Divination |
| ユーザー |
 rrrriki
|
| 提出日時 | 2024-08-03 19:32:24 |
| 言語 | C++17 (gcc 12.3.0 + boost 1.83.0) |
| 結果 |
AC
|
| 実行時間 | 30 ms / 2,000 ms |
| コード長 | 58,673 bytes |
| コンパイル時間 | 3,083 ms |
| コンパイル使用メモリ | 251,044 KB |
| 実行使用メモリ | 6,944 KB |
| 最終ジャッジ日時 | 2024-08-03 19:32:28 |
| 合計ジャッジ時間 | 4,458 ms |
|
ジャッジサーバーID （参考情報） |
judge3 / judge2 |
（要ログイン）
テストケース
テストケース表示| 入力 | 結果 | 実行時間 実行使用メモリ |
|---|---|---|
| testcase_00 | AC | 1 ms
6,812 KB |
| testcase_01 | AC | 2 ms
6,940 KB |
| testcase_02 | AC | 1 ms
6,940 KB |
| testcase_03 | AC | 2 ms
6,940 KB |
| testcase_04 | AC | 2 ms
6,944 KB |
| testcase_05 | AC | 30 ms
6,944 KB |
| testcase_06 | AC | 27 ms
6,940 KB |
| testcase_07 | AC | 22 ms
6,940 KB |
| testcase_08 | AC | 21 ms
6,944 KB |
| testcase_09 | AC | 17 ms
6,940 KB |
| testcase_10 | AC | 2 ms
6,940 KB |
| testcase_11 | AC | 2 ms
6,944 KB |
| testcase_12 | AC | 2 ms
6,944 KB |
| testcase_13 | AC | 1 ms
6,940 KB |
| testcase_14 | AC | 2 ms
6,944 KB |
| testcase_15 | AC | 2 ms
6,944 KB |
| testcase_16 | AC | 2 ms
6,944 KB |
| testcase_17 | AC | 1 ms
6,944 KB |
| testcase_18 | AC | 2 ms
6,940 KB |
| testcase_19 | AC | 2 ms
6,944 KB |
| testcase_20 | AC | 2 ms
6,944 KB |
| testcase_21 | AC | 2 ms
6,940 KB |
| testcase_22 | AC | 2 ms
6,944 KB |
| testcase_23 | AC | 2 ms
6,940 KB |
| testcase_24 | AC | 2 ms
6,940 KB |
| testcase_25 | AC | 2 ms
6,940 KB |
| testcase_26 | AC | 1 ms
6,944 KB |
| testcase_27 | AC | 2 ms
6,944 KB |
| testcase_28 | AC | 1 ms
6,940 KB |
| testcase_29 | AC | 2 ms
6,940 KB |
| testcase_30 | AC | 2 ms
6,944 KB |
| testcase_31 | AC | 2 ms
6,944 KB |
| testcase_32 | AC | 1 ms
6,940 KB |
| testcase_33 | AC | 2 ms
6,940 KB |
| testcase_34 | AC | 2 ms
6,944 KB |
ソースコード
/**
* author: rrrriki
* created: 02.08.2024 21:23:22
*/
//#define USE_ACL
//#define USE_BOOST
#if !__INCLUDE_LEVEL__
#include <bits/stdc++.h>
#include __FILE__
unordered_map<int, ll> calc_nums(int N) {
auto divs = divisor(N);
reverse(all(divs));
unordered_map<int, ll> nums;
for (auto d : divs) {
if (N < d) continue;
ll num = N / d; // dの倍数の数
for (auto d2 : divs) {
if (d2 <= d) break;
if (d2 % d != 0) continue; //
//ll num2 = nums[d2]; // d2の倍数の数
num -= nums[d2];
}
nums[d] = num;
}
return nums;
}
int main() {
cin.tie(0);
ios_base::sync_with_stdio(false);
ll N;
cin >> N;
unordered_map<int, long double> dp;
dp[1] = 0;
auto calc = [&](auto calc, int n) -> long double {
if (dp.count(n)) return dp[n];
auto nums = calc_nums(n);
long double res = 1;
for (auto [d, num] : nums) {
if (d == n) continue;
//cout << d << " " << num << endl;
res += (long double)num / n * calc(calc, d);
}
//cout << res << endl;
//cout << n << " " << res << endl;
res *= (long double)n / (n - 1);
//cout << n << " " << res << endl;
return dp[n] = res;
};
cout << fixed << setprecision(15) << calc(calc, N) << endl;
return 0;
}
#else
// clang-format off
#ifdef USE_ACL
#include <atcoder/all>
using namespace atcoder;
using mint = modint998244353;
//using mint = modint1000000007;
#endif
#ifdef USE_BOOST
#include <boost/algorithm/string/classification.hpp>
#include <boost/algorithm/string/compare.hpp>
#include <boost/algorithm/string/join.hpp>
#include <boost/algorithm/string/replace.hpp>
#include <boost/algorithm/string/split.hpp>
#include <boost/algorithm/string/trim.hpp>
#include <boost/dynamic_bitset.hpp>
#include <boost/integer/extended_euclidean.hpp>
#include <boost/math/tools/minima.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/cpp_int.hpp>
using namespace boost::multiprecision;
#endif
#define all(x) x.begin(), x.end()
#define YES cout << "Yes\n"
#define NO cout << "No\n"
using namespace std;
#ifdef LOCAL
#include "debug.h"
#else
#define dbg(...) 42
#endif
using ll = long long;
using vl = vector<ll>;
/// コンテナの全出力 @tparam T コンテナの型 @param A コンテナ @param gap 区切り文字
template <class T> void out_c(T &A, string gap=" ") {auto itr = A.begin(); if (itr != A.end()) {cout << *itr; itr++;} while (itr != A.end()) {cout << gap << *itr; itr++;}cout << "\n"; return;}
template <class T> void out_c_pairs(T &A, string gap_inside=" ", string gap_outside = " ") {auto itr = A.begin();if (itr != A.end()) {cout << itr->first << gap_inside << itr->second;itr++;}while (itr != A.end()) {cout << gap_outside << itr->first << gap_inside << itr->second;itr++;}cout << "\n";return;}
/// べき乗を誤差なく計算する @param x 底 @param n 指数 @return x^n
ll _pow(ll x, ll n) {if (n == 0) return 1; ll val = _pow(x, n / 2); val *= val; if (n & 1) val *= x; return val;}
// マンハッタン距離
template <class T> T mnht(T a, T b, T c, T d) {return abs(a - c) + abs(b - d);}
/// ランレングス圧縮 @param s 圧縮する文字列 @return 圧縮した文字列を格納したvector<pair<char, int>>
vector<pair<char, int>> rle(const string &s){vector<pair<char, int>> vec;int cnt = 1; for(int i = 1; i < (int)s.size(); i++) {if(s[i] != s[i-1]){vec.emplace_back(s[i-1], cnt); cnt = 0;} cnt++;} vec.emplace_back(s.back(), cnt);return vec;}
/// ランレングス圧縮 @tparam T 圧縮するvectorの型 @param v 圧縮するvector @return 圧縮したvectorを格納したvector<pair<T, int>>
template <class T> vector<pair<T, int>> rle(const vector<T> &v) {vector<pair<T, int>> vec;int cnt = 1; for(int i = 1; i < (int)v.size(); i++) {if(v[i] != v[i-1]){vec.emplace_back(v[i-1], cnt); cnt = 0;} cnt++;} vec.emplace_back(v.back(), cnt);return vec;}
// 素数
bool is_prime(ll x){for (ll i=2; i*i<=x; i++){if(x%i==0)return false;}return true;}
map<ll,int> prime_factor(ll n) {map<ll,int> ret; for(ll i=2; i*i <= n; i++) {while(n%i == 0) {ret[i]++; n /= i;}} if(n != 1) ret[n]=1;return ret;}
vector<bool> sieve_of_era(ll n) {vector<bool> ret(n+1,true); ret[0]=false; ret[1]=false; for(ll i=2; i*i<=n; i++) {if(ret[i]) {for(ll j=i*2; j<=n; j+=i) {ret[j]=false;}}} return ret;}
// 約数全列挙
vector<ll> divisor(ll n) {vector<ll> ret; for(ll i=1; i*i <= n; i++) {if(n%i == 0) {ret.push_back(i); if(i*i != n) ret.emplace_back(n/i);}} sort(all(ret)); return ret;}
// 切り捨て、切り上げ、外側
inline constexpr ll ceil_div(const ll a, const ll b) {return (a + b - 1) / b - ((a + b - 1) % b < 0);}
inline constexpr ll floor_div(const ll a, const ll b) {return a / b - (a % b < 0);}
inline constexpr ll out_div(ll x, ll y) {ll d = x / y; return d * y == x ? d : ((x > 0) == (y > 0)) ? d + 1 : d - 1;}
/// 組み合わせの全探索 @param k 組み合わせの要素数
template <typename T> bool next_combination(const T first, const T last, int k) {
const T subset = first + k;
// empty container | k = 0 | k == n
if (first == last || first == subset || last == subset) {
return false;
}
T src = subset;
while (first != src) {
src--;
if (*src < *(last - 1)) {
T dest = subset;
while (*src >= *dest) {
dest++;
}
iter_swap(src, dest);
rotate(src + 1, dest + 1, last);
rotate(subset, subset + (last - dest) - 1, last);
return true;
}
}
// restore
rotate(first, subset, last);
return false;
}
// グラフ
/**
* @brief Edgeクラスはグラフのエッジ（辺）を表します。
*
* @tparam T エッジの重みの型（デフォルトはint）
*/
template <typename T = int> struct Edge {
int from, to; // エッジの始点と終点
T cost; // エッジの重み
int idx; // エッジのインデックス（オプション）
// デフォルトコンストラクタ
Edge() = default;
// エッジをコストに基づいて比較するための演算子オーバーロード
bool operator<(const Edge &other) const { return cost < other.cost; }
bool operator>(const Edge& other) const { return cost > other.cost; }
friend std::ostream& operator<<(std::ostream& os, const Edge& edge) { os << edge.to; return os; }
// コンストラクタ
Edge(int from, int to, T cost = 1, int idx = -1)
: from(from), to(to), cost(cost), idx(idx) {}
// エッジの終点をintとして取得するためのキャスト演算子
operator int() const { return to; }
};
/**
* @brief Graphクラスはグラフのデータ構造を表します。
* @tparam T エッジの重みの型（デフォルトはint）
*/
template <typename T = int> struct Graph {
vector<vector<Edge<T>>> g; // 各ノードから出ているエッジのリスト
int es; // エッジの数
// デフォルトコンストラクタ
Graph() = default;
// ノード数nを指定するコンストラクタ
explicit Graph(int n) : g(n), es(0) {}
// グラフのサイズ（ノードの数）を返す
size_t size() const { return g.size(); }
// 有向エッジを追加する関数
void add_directed_edge(int from, int to, T cost = 1) {
g[from].emplace_back(from, to, cost, es++);
}
// 無向エッジを追加する関数
void add_edge(int from, int to, T cost = 1) {
g[from].emplace_back(from, to, cost, es);
g[to].emplace_back(to, from, cost, es++);
}
/// @brief エッジを読み込む関数 @param M エッジの数 @param padding インデックスのオフセット @param weighted 重み付きかどうか @param directed 有向かどうか
void read(int M, int padding = -1, bool weighted = false,
bool directed = false) {
for (int i = 0; i < M; i++) {
int a, b;
cin >> a >> b;
a += padding;
b += padding;
T c = T(1);
if (weighted) cin >> c;
if (directed)
add_directed_edge(a, b, c);
else
add_edge(a, b, c);
}
}
// 演算子オーバーロード：インデックスによるエッジのリストへのアクセス
inline vector<Edge<T>> &operator[](const int &k) { return g[k]; }
// 演算子オーバーロード（const版）：インデックスによるエッジのリストへのアクセス
inline const vector<Edge<T>> &operator[](const int &k) const { return g[k]; }
};
/// @brief エッジのリスト @tparam T エッジの重みの型
template <typename T = int> using Edges = vector<Edge<T>>;
// ダイクストラ法
/**
* @brief dijkstra関数はダイクストラ法を用いて最短経路を求める関数です。
* @tparam T エッジの重みの型
* @param g グラフ
* @param s 始点
* @return vector<T> 始点から各頂点への最短経路の長さ
* @note 計算量はO((E+V)logV)
*/
template <typename T> vector<T> dijkstra(Graph<T> &g, int s) {
vector<T> dist(g.size(), numeric_limits<T>::max());
dist[s] = 0;
priority_queue<pair<T, int>, vector<pair<T, int>>, greater<pair<T, int>>> pq;
pq.emplace(0, s);
while (!pq.empty()) {
auto [d, v] = pq.top();
pq.pop();
if (dist[v] < d) continue;
for (auto e : g[v]) {
if (dist[e.to] > dist[v] + e.cost) {
dist[e.to] = dist[v] + e.cost;
pq.emplace(dist[e.to], e.to);
}
}
}
return dist;
}
#ifndef USE_ACL
struct dsu {
public:
dsu() : _n(0) {}
explicit dsu(int n) : _n(n), parent_or_size(n, -1) {}
int merge(int a, int b) {
assert(0 <= a && a < _n);
assert(0 <= b && b < _n);
int x = leader(a), y = leader(b);
if (x == y) return x;
if (-parent_or_size[x] < -parent_or_size[y]) std::swap(x, y);
parent_or_size[x] += parent_or_size[y];
parent_or_size[y] = x;
return x;
}
bool same(int a, int b) {
assert(0 <= a && a < _n);
assert(0 <= b && b < _n);
return leader(a) == leader(b);
}
int leader(int a) {
assert(0 <= a && a < _n);
if (parent_or_size[a] < 0) return a;
return parent_or_size[a] = leader(parent_or_size[a]);
}
int size(int a) {
assert(0 <= a && a < _n);
return -parent_or_size[leader(a)];
}
std::vector<std::vector<int>> groups() {
std::vector<int> leader_buf(_n), group_size(_n);
for (int i = 0; i < _n; i++) {
leader_buf[i] = leader(i);
group_size[leader_buf[i]]++;
}
std::vector<std::vector<int>> result(_n);
for (int i = 0; i < _n; i++) {
result[i].reserve(group_size[i]);
}
for (int i = 0; i < _n; i++) {
result[leader_buf[i]].push_back(i);
}
result.erase(
std::remove_if(result.begin(), result.end(),
[&](const std::vector<int>& v) { return v.empty(); }),
result.end());
return result;
}
private:
int _n;
// root node: -1 * component size
// otherwise: parent
std::vector<int> parent_or_size;
};
#endif
/// @brief 重み付きUnionFind @tparam T 重みの型
template <class T> struct UnionFindWithPotential {
vector<int> dat; // 親の番号 根の場合は-1
vector<T> pot; // 親との差分
UnionFindWithPotential(int N) : dat(N, -1), pot(N, T()) {}
/// @brief xの根を求める @return 根
int root(int x) {
if (dat[x] < 0) return x;
int r = root(dat[x]);
pot[x] += pot[dat[x]];
return dat[x] = r;
}
/// @brief xのポテンシャルを求める @return P(x) - P(root(x))
T potential(int x) {
root(x);
return pot[x];
}
bool same(int x, int y) { return root(x) == root(y); }
/// @brief xとyのポテンシャルの差を求める @return P(x) - P(y)
T diff(int x, int y) { return potential(x) - potential(y); }
/// @brief P(x) = P(y) + p でマージ @param p ポテンシャルの差 @return マージできたかどうか
bool merge(int x, int y, T p) {
p += potential(y) - potential(x);
x = root(x), y = root(y);
if (x == y) return p == T();
if (dat[x] < dat[y]) swap(x, y), p = -p;
dat[y] += dat[x];
dat[x] = y;
pot[x] = p;
return true;
}
/// @brief xが属する連結成分のサイズを求める @return xが属する連結成分のサイズ
int size(int x) { return -dat[root(x)]; }
};
/**
* @brief krsukal関数はクラスカル法を用いて最小/最大全域木を求める関数です。
* @tparam T エッジの重みの型
* @param g グラフ
* @param s 最小全域木を求める場合は"min"、最大全域木を求める場合は"max"を指定
* @return T 最小/最大全域木の重み
* @note 計算量はO(ElogV)
*/
template <typename T> T kruskal(Graph<T> &g, string s = "min") {
T res = 0;
int n = g.size();
dsu dsu(n);
Edges<T> edges;
for (int i = 0; i < n; i++) {
for (auto e : g[i]) {
edges.emplace_back(e);
}
}
if (s == "max") sort(all(edges), greater<Edge<T>>());
else sort(all(edges));
for (auto e : edges) {
if (dsu.same(e.from, e.to)) continue;
dsu.merge(e.from, e.to);
res += e.cost;
}
return res;
}
/**
* @brief CumulativeSum2Dクラスは二次元累積和を計算するためのクラスです。
* @tparam T 累積和の型
*/
template <class T> struct CumulativeSum2D {
vector<vector<T> > data;
/// @brief W x H の二次元累積和を宣言 @param W 幅 @param H 高さ
CumulativeSum2D(int W, int H) : data(W + 1, vector<T>(H + 1, 0)) {}
/// @brief x, y に z を加算 @param x x座標 @param y y座標 @param z 加算する値
void add(int x, int y, T z) {
++x, ++y;
if (x >= (int)data.size() || y >= (int)data[0].size()) return;
data[x][y] += z;
}
/// @brief 二次元累積和を構築
void build() {
for (int i = 1; i < (int)data.size(); i++) {
for (int j = 1; j < (int)data[i].size(); j++) {
data[i][j] += data[i][j - 1] + data[i - 1][j] - data[i - 1][j - 1];
}
}
}
/// @brief (sx, sy) から (gx, gy) までの和を計算 [sx, gx), [sy, gy)
/// @param sx x座標の始点 @param sy y座標の始点 @param gx x座標の終点 @param gy y座標の終点 @note gxとgyは含まれない
T query(int sx, int sy, int gx, int gy) const {
return (data[gx][gy] - data[sx][gy] - data[gx][sy] + data[sx][sy]);
}
};
/// @brief 座標圧縮 @tparam T 座標の型 @param compressed 座標圧縮用のmap @param reverse_compressed 圧縮後から圧縮前の座標を取得するmap @param vectors 圧縮する座標 @return tuple<圧縮された座標>
template <typename T, typename... Vectors> auto shrink_coord(
unordered_map<T, int>& compressed,
unordered_map<int, T>& reverse_compressed,
Vectors&... vectors) {
// 座標のsetを作成
set<ll> coords;
(coords.insert(vectors.begin(), vectors.end()), ...);
// 座標圧縮用のmapを作成
int idx = 0;
for (const auto& coord : coords) {
compressed[coord] = idx;
reverse_compressed[idx] = coord;
++idx;
}
// 各vectorを圧縮
auto compress_vector = [&](auto& vec) {
vector<int> compressedVec(vec.size());
transform(vec.begin(), vec.end(), compressedVec.begin(),
[&](const ll& val) { return compressed.at(val); });
return compressedVec;
};
// 圧縮されたvectorのtupleを返す
return make_tuple(compress_vector(vectors)...);
};
// ヘルパーテンプレートを定義して、テンプレートパラメータパックの要素型を取得
template<typename T>
struct element_type {
using type = typename decay<decltype(*begin(declval<T>()))>::type;
};
template<typename... Vectors>
using element_type_t = typename element_type<Vectors...>::type;
// オーバーロード: compressed map を省略可能にする
template <typename... Vectors>
auto shrink_coord(Vectors&... vectors) {
using T = element_type_t<decay_t<decltype(vectors)>...>;
unordered_map<T, int> compressed;
unordered_map<int, T> reverse_compressed;
return shrink_coord(compressed, reverse_compressed, vectors...);
}
/// @brief 座標圧縮 @tparam T 座標の型 @param vectors 圧縮する座標 @param compressed 座標圧縮用のmap @param reverse_compressed 圧縮後から圧縮前の座標を取得するmap @return 圧縮された座標
template <typename T>
vector<vector<int>> shrink_coord(
vector<vector<T>>& vectors,
unordered_map<T, int>& compressed = unordered_map<T, int>(),
unordered_map<int, T>& reverse_compressed = unordered_map<int, T>()) {
// 座標のsetを作成
set<T> coords;
for (const auto& vec : vectors) {
coords.insert(vec.begin(), vec.end());
}
// 座標圧縮用のmapを作成
int idx = 0;
for (const auto& coord : coords) {
compressed[coord] = idx;
reverse_compressed[idx] = coord;
++idx;
}
// 各vectorを圧縮
vector<vector<int>> result(vectors.size());
for (size_t i = 0; i < vectors.size(); ++i) {
result[i].resize(vectors[i].size());
transform(vectors[i].begin(), vectors[i].end(), result[i].begin(),
[&](const T& val) { return compressed.at(val); });
}
// 圧縮されたvectorを返す
return result;
}
/**
* @brief NQueenクラスはN-Queen問題を解くためのクラスです。
*/
struct NQueen {
public:
explicit NQueen(int n) : N(n) { assert(n > 0);}
/// @brief クイーンを置く @param x x座標 @param y y座標 @return クイーンを置けたかどうか
bool add_queen(int x, int y) {
if (row.count(x) || col.count(y) || diag1.count(x + y) || diag2.count(x - y)) return false;
assert(x < N && y < N);
assert(x >= 0 && y >= 0);
queens.emplace(x, y);
row.insert(x); // x が同じ行にある
col.insert(y); // y が同じ列にある
diag1.insert(x + y); // x + y が同じ斜めの利き筋にある
diag2.insert(x - y); // x - y が同じ斜めの利き筋にある
return true;
}
/// @brief クイーンを取り除く @param x x座標 @param y y座標 @return クイーンを取り除けたかどうか
bool remove_queen(int x, int y) {
if (!row.count(x) || !col.count(y) || !diag1.count(x + y) || !diag2.count(x - y)) return false;
assert(x < N && y < N);
assert(x >= 0 && y >= 0);
queens.erase({x, y});
if (added_queens.count({x, y})) added_queens.erase({x, y});
row.erase(x);
col.erase(y);
diag1.erase(x + y);
diag2.erase(x - y);
return true;
}
/// @brief x, yにクイーンを置けるかどうか @param x x座標 @param y y座標 @return クイーンを置けるかどうか
bool is_valid(int x, int y) { return !row.count(x) && !col.count(y) && !diag1.count(x + y) && !diag2.count(x - y);}
/// @brief クイーンが全てのマスを利き筋に置いているかどうか
bool is_valid() { return (int)row.size() == N;}
/// @brief x行目以降のクイーンを置く @return クイーンを置けたかどうか
bool solve(int x = 0) {
if (x == N) return true;
if (is_valid()) return true;
if (row.count(x)) return solve(x + 1);
for (int y = 0; y < N; y++) {
if (is_valid(x, y)) {
add_queen(x, y);
added_queens.emplace(x, y);
if (solve(x + 1)) return true;
remove_queen(x, y);
}
}
return false;
}
/// @brief チェス盤のサイズを返す
int size() { return N; }
/// @brief クイーンの位置を出力する
friend std::ostream& operator<<(std::ostream& os, const NQueen& nqueen) {
for (auto [x, y] : nqueen.queens) os << x << " " << y << "\n";
return os;
}
/// @brief クイーンの位置をvector<pair<int,int>>に出力する
vector<pair<int, int> > get_queens() { return vector<pair<int, int> >(queens.begin(), queens.end());}
/// @brief 追加したクイーンの位置をvector<pair<int,int>>に出力する
vector<pair<int, int> > get_added_queens() { return vector<pair<int, int> >(added_queens.begin(), added_queens.end());}
/// @brief 盤面を出力する @param c クイーンの文字 @param d 空白の文字 @note デフォルトは 'Q' と '.'
void print(char c = 'Q', char d = '.') {
vector<vector<char> > board(N, vector<char>(N, d));
for (auto [x, y] : queens) {
board[x][y] = c;
}
for (auto& row : board) {
for (auto& c : row) {
cout << c;
}
cout << "\n";
}
}
private:
int N; // チェス盤のサイズ
unordered_set<int> row, col, diag1, diag2; // それぞれの行、列、斜めの利き筋にクイーンがあるかどうか
set<pair<int, int> > queens; // クイーンの位置
set<pair<int, int> > added_queens; // 追加したクイーンの位置
};
/**
* @brief Rolling-Hash(ローリングハッシュ)
* @see https://qiita.com/keymoon/items/11fac5627672a6d6a9f6
* @see https://ei1333.github.io/library/string/rolling-hash.hpp
*/
struct RollingHash {
static const uint64_t mod = (1ull << 61ull) - 1;
using uint128_t = __uint128_t;
const uint64_t base;
vector< uint64_t > power;
/// @brief 加算 @param a 加数 @param b 加数 @return 和
static inline uint64_t add(uint64_t a, uint64_t b) {
if((a += b) >= mod) a -= mod;
return a;
}
/// @brief 乗算 @param a 乗数 @param b 乗数 @return 積
static inline uint64_t mul(uint64_t a, uint64_t b) {
uint128_t c = (uint128_t) a * b;
return add(c >> 61, c & mod);
}
/// @brief 2^61-1 未満の乱数を生成する
static inline uint64_t generate_base() {
mt19937_64 mt(chrono::steady_clock::now().time_since_epoch().count());
uniform_int_distribution< uint64_t > rand(1, RollingHash::mod - 1);
return rand(mt);
}
/// @brief ハッシュテーブルのサイズを拡張する @param sz 拡張するサイズ
inline void expand(size_t sz) {
if(power.size() < sz + 1) {
int pre_sz = (int) power.size();
power.resize(sz + 1);
for(int i = pre_sz - 1; i < (int)sz; i++) {
power[i + 1] = mul(power[i], base);
}
}
}
explicit RollingHash(uint64_t base = generate_base()) : base(base), power{1} {}
/// @brief 文字列sのハッシュ値を計算する @param s 文字列 @return ハッシュ値
vector< uint64_t > build(const string &s) const {
int sz = s.size();
vector< uint64_t > hashed(sz + 1);
for(int i = 0; i < sz; i++) {
hashed[i + 1] = add(mul(hashed[i], base), s[i]);
}
return hashed;
}
/// @brief ベクタsのハッシュ値を計算する @tparam T ベクタの型 @param s ベクタ @return ハッシュ値
template< typename T >
vector< uint64_t > build(const vector< T > &s) const {
int sz = s.size();
vector< uint64_t > hashed(sz + 1);
for(int i = 0; i < sz; i++) {
hashed[i + 1] = add(mul(hashed[i], base), s[i]);
}
return hashed;
}
/// @brief 文字列sの[l, r)のハッシュ値を計算する @param s 文字列 @param l 左端 @param r 右端 @return ハッシュ値
uint64_t query(const vector< uint64_t > &s, int l, int r) {
expand(r - l);
return add(s[r], mod - mul(s[l], power[r - l]));
}
/// @brief ハッシュ値h1とハッシュ値h2を結合する @param h1 ハッシュ値1 @param h2 ハッシュ値2 @param h2len ハッシュ値2の長さ @return 結合されたハッシュ値
uint64_t combine(uint64_t h1, uint64_t h2, size_t h2len) {
expand(h2len);
return add(mul(h1, power[h2len]), h2);
}
/// @brief ハッシュテーブルaの区間[l1,r1)と、ハッシュテーブルbの区間[l2,r2)の最長共通接頭辞の長さを求める @param a ハッシュテーブルa @param l1 左端 @param r1 右端 @param b ハッシュテーブルb @param l2 左端 @param r2 右端 @return 最長共通接頭辞の長さ
int lcp(const vector< uint64_t > &a, int l1, int r1, const vector< uint64_t > &b, int l2, int r2) {
int len = min(r1 - l1, r2 - l2);
int low = 0, high = len + 1;
while(high - low > 1) {
int mid = (low + high) / 2;
if(query(a, l1, l1 + mid) == query(b, l2, l2 + mid)) low = mid;
else high = mid;
}
return low;
}
};
/**
* @brief K-Shortest-PathをYen’s Algorithm により求める関数
* @tparam T グラフの重みの型 @param g グラフ @param s 始点 @param t 終点 @param k 最短経路の数
* @return vector<pair<T, vector<int>>> 最短経路の長さと経路
*/
template <typename T>
vector<pair<T, vector<int>>> k_shortest_path(const Graph<T> &g, int s, int t, int k) {
assert(s != t);
int N = (int)g.size();
int M = 0;
for (int i = 0; i < N; i++) M += (int)g[i].size();
vector<int> latte(M), malta(M);
vector<T> cost(M);
for (int i = 0; i < N; i++) {
for (auto &e : g[i]) {
latte[e.idx] = i;
malta[e.idx] = e.to;
cost[e.idx] = e.cost;
}
}
const auto INF = numeric_limits<T>::max();
vector<int> dame(M, -1);
int timestamp = 0;
// dijkstra
auto shortest_path = [&](vector<T> &dist, vector<int> &from, vector<int> &id, int st) {
using Pi = pair<T, int>;
priority_queue<Pi, vector<Pi>, greater<>> que;
que.emplace(dist[st], st);
while (!que.empty()) {
T cost;
int idx;
tie(cost, idx) = que.top();
que.pop();
if (dist[idx] < cost) continue;
if (idx == t) return;
for (auto &e : g[idx]) {
auto next_cost = cost + e.cost;
if (dist[e.to] <= next_cost) continue;
if (dame[e.idx] == timestamp) continue;
dist[e.to] = next_cost;
from[e.to] = idx;
id[e.to] = e.idx;
que.emplace(dist[e.to], e.to);
}
}
};
auto restore = [](const vector<int> &es, const vector<int> &vs, int from,
int to) {
vector<int> tap;
while (to != from) {
tap.emplace_back(es[to]);
to = vs[to];
}
reverse(begin(tap), end(tap));
return tap;
};
vector<T> dist(g.size(), INF);
vector<int> from(g.size(), -1), id(g.size(), -1);
dist[s] = 0;
shortest_path(dist, from, id, s);
if (dist[t] == INF) return {};
vector<pair<T, vector<int> > > A;
set<pair<T, vector<int> > > B;
A.emplace_back(dist[t], restore(id, from, s, t));
for (int i = 1; i < k; i++) {
dist.assign(g.size(), INF);
from.assign(g.size(), -1);
id.assign(g.size(), -1);
dist[s] = 0;
vector<int> candidate(A.size());
iota(begin(candidate), end(candidate), 0);
auto &last_path = A.back().second;
int cur = s;
for (int j = 0; j < last_path.size(); j++) {
for (auto &k : candidate) {
if (j < A[k].second.size()) dame[A[k].second[j]] = timestamp;
}
vector<T> dist2{dist};
vector<int> from2{from}, id2{id};
shortest_path(dist2, from2, id2, cur);
++timestamp;
if (dist2[t] != INF) {
auto path = restore(id2, from2, s, t);
bool ok = true;
for (auto &p : candidate) {
if (path == A[p].second) {
ok = false;
break;
}
}
if (ok) B.emplace(dist2[t], path);
}
vector<int> accept;
for (auto &k : candidate) {
if (j < A[k].second.size() && A[k].second[j] == last_path[j]) {
accept.emplace_back(k);
}
}
dist[malta[last_path[j]]] =
dist[latte[last_path[j]]] + cost[last_path[j]];
from[malta[last_path[j]]] = latte[last_path[j]];
id[malta[last_path[j]]] = last_path[j];
cur = malta[last_path[j]];
candidate = move(accept);
}
if (B.size()) {
A.emplace_back(*B.begin());
B.erase(B.begin());
}
}
return A;
}
// ---------------------------------------
// ----- segment trees by @rrrrikiOW -----
// ---------------------------------------
// ----- Based on AtCoder Library --------
// -------------- VER.1.1.0 --------------
// ----- Last Update: 2024/03/03 ---------
// ---------------------------------------
/// @brief 2の冪に切り上げる @param n 数 @return 2の冪
int ceil_pow2(int n) {
int x = 0;
while ((1U << x) < (unsigned int)(n)) x++;
return x;
}
#ifndef USE_ACL
/// @brief セグメント木
/// @tparam S セグメント木の型 @tparam op セグメント木の演算 @tparam e セグメント木の単位元
template <class S, S (*op)(S, S), S (*e)()> struct segtree {
public:
segtree() : segtree(0) {}
explicit segtree(int n) : segtree(std::vector<S>(n, e())) {}
explicit segtree(const std::vector<S>& v) : _n(int(v.size())) {
log = ceil_pow2(_n);
size = 1 << log;
d = std::vector<S>(2 * size, e());
for (int i = 0; i < _n; i++) d[size + i] = v[i];
for (int i = size - 1; i >= 1; i--) {
update(i);
}
}
/// @brief 0-indexed で k 番目の要素を x に変更する O(logN)
void set(int p, S x) {
assert(0 <= p && p < _n);
p += size;
d[p] = x;
for (int i = 1; i <= log; i++) update(p >> i);
}
/// @brief 0-indexed で k 番目の要素を取得する O(logN)
S get(int p) const {
assert(0 <= p && p < _n);
return d[p + size];
}
/// @brief op(a[l], ..., a[r - 1]) を、モノイドの性質を満たしていると仮定して計算します。 l=r のときは e() を返します。 O(logN)
S prod(int l, int r) const {
assert(0 <= l && l <= r && r <= _n);
S sml = e(), smr = e();
l += size;
r += size;
while (l < r) {
if (l & 1) sml = op(sml, d[l++]);
if (r & 1) smr = op(d[--r], smr);
l >>= 1;
r >>= 1;
}
return op(sml, smr);
}
/// @brief op(a[0], ..., a[n - 1]) を、モノイドの性質を満たしていると仮定して計算します O(1)
S all_prod() const { return d[1]; }
template <bool (*f)(S)> int max_right(int l) const {
return max_right(l, [](S x) { return f(x); });
}
template <class F> int max_right(int l, F f) const {
assert(0 <= l && l <= _n);
assert(f(e()));
if (l == _n) return _n;
l += size;
S sm = e();
do {
while (l % 2 == 0) l >>= 1;
if (!f(op(sm, d[l]))) {
while (l < size) {
l = (2 * l);
if (f(op(sm, d[l]))) {
sm = op(sm, d[l]);
l++;
}
}
return l - size;
}
sm = op(sm, d[l]);
l++;
} while ((l & -l) != l);
return _n;
}
template <bool (*f)(S)> int min_left(int r) const {
return min_left(r, [](S x) { return f(x); });
}
template <class F> int min_left(int r, F f) const {
assert(0 <= r && r <= _n);
assert(f(e()));
if (r == 0) return 0;
r += size;
S sm = e();
do {
r--;
while (r > 1 && (r % 2)) r >>= 1;
if (!f(op(d[r], sm))) {
while (r < size) {
r = (2 * r + 1);
if (f(op(d[r], sm))) {
sm = op(d[r], sm);
r--;
}
}
return r + 1 - size;
}
sm = op(d[r], sm);
} while ((r & -r) != r);
return 0;
}
private:
int _n, size, log;
std::vector<S> d;
void update(int k) { d[k] = op(d[2 * k], d[2 * k + 1]); }
};
/// @brief 遅延セグメント木
/// @tparam S セグメント木の型 @tparam op セグメント木の演算 @tparam e セグメント木の単位元
/// @tparam F 作用素の型 @tparam mapping 作用素の演算 @tparam composition 作用素の合成 @tparam id 作用素の単位元
template <class S,
S (*op)(S, S),
S (*e)(),
class F,
S (*mapping)(F, S),
F (*composition)(F, F),
F (*id)()>
struct lazy_segtree {
public:
lazy_segtree() : lazy_segtree(0) {}
explicit lazy_segtree(int n) : lazy_segtree(std::vector<S>(n, e())) {}
explicit lazy_segtree(const std::vector<S>& v) : _n(int(v.size())) {
log = ceil_pow2(_n);
size = 1 << log;
d = std::vector<S>(2 * size, e());
lz = std::vector<F>(size, id());
for (int i = 0; i < _n; i++) d[size + i] = v[i];
for (int i = size - 1; i >= 1; i--) {
update(i);
}
}
/// @brief 0-indexed で k 番目の要素を x に変更する O(logN)
void set(int p, S x) {
assert(0 <= p && p < _n);
p += size;
for (int i = log; i >= 1; i--) push(p >> i);
d[p] = x;
for (int i = 1; i <= log; i++) update(p >> i);
}
/// @brief 0-indexed で k 番目の要素を取得する O(logN)
S get(int p) {
assert(0 <= p && p < _n);
p += size;
for (int i = log; i >= 1; i--) push(p >> i);
return d[p];
}
/// @brief op(a[l], ..., a[r - 1]) を、モノイドの性質を満たしていると仮定して計算します。 l=r のときは e() を返します。 O(logN)
S prod(int l, int r) {
assert(0 <= l && l <= r && r <= _n);
if (l == r) return e();
l += size;
r += size;
for (int i = log; i >= 1; i--) {
if (((l >> i) << i) != l) push(l >> i);
if (((r >> i) << i) != r) push((r - 1) >> i);
}
S sml = e(), smr = e();
while (l < r) {
if (l & 1) sml = op(sml, d[l++]);
if (r & 1) smr = op(d[--r], smr);
l >>= 1;
r >>= 1;
}
return op(sml, smr);
}
/// @brief op(a[0], ..., a[n - 1]) を、モノイドの性質を満たしていると仮定して計算します O(1)
S all_prod() { return d[1]; }
/// @brief a[p] = f(a[p])
void apply(int p, F f) {
assert(0 <= p && p < _n);
p += size;
for (int i = log; i >= 1; i--) push(p >> i);
d[p] = mapping(f, d[p]);
for (int i = 1; i <= log; i++) update(p >> i);
}
/// @brief [l, r) の要素に f を作用させます O(logN)
void apply(int l, int r, F f) {
assert(0 <= l && l <= r && r <= _n);
if (l == r) return;
l += size;
r += size;
for (int i = log; i >= 1; i--) {
if (((l >> i) << i) != l) push(l >> i);
if (((r >> i) << i) != r) push((r - 1) >> i);
}
{
int l2 = l, r2 = r;
while (l < r) {
if (l & 1) all_apply(l++, f);
if (r & 1) all_apply(--r, f);
l >>= 1;
r >>= 1;
}
l = l2;
r = r2;
}
for (int i = 1; i <= log; i++) {
if (((l >> i) << i) != l) update(l >> i);
if (((r >> i) << i) != r) update((r - 1) >> i);
}
}
template <bool (*g)(S)> int max_right(int l) {
return max_right(l, [](S x) { return g(x); });
}
template <class G> int max_right(int l, G g) {
assert(0 <= l && l <= _n);
assert(g(e()));
if (l == _n) return _n;
l += size;
for (int i = log; i >= 1; i--) push(l >> i);
S sm = e();
do {
while (l % 2 == 0) l >>= 1;
if (!g(op(sm, d[l]))) {
while (l < size) {
push(l);
l = (2 * l);
if (g(op(sm, d[l]))) {
sm = op(sm, d[l]);
l++;
}
}
return l - size;
}
sm = op(sm, d[l]);
l++;
} while ((l & -l) != l);
return _n;
}
template <bool (*g)(S)> int min_left(int r) {
return min_left(r, [](S x) { return g(x); });
}
template <class G> int min_left(int r, G g) {
assert(0 <= r && r <= _n);
assert(g(e()));
if (r == 0) return 0;
r += size;
for (int i = log; i >= 1; i--) push((r - 1) >> i);
S sm = e();
do {
r--;
while (r > 1 && (r % 2)) r >>= 1;
if (!g(op(d[r], sm))) {
while (r < size) {
push(r);
r = (2 * r + 1);
if (g(op(d[r], sm))) {
sm = op(d[r], sm);
r--;
}
}
return r + 1 - size;
}
sm = op(d[r], sm);
} while ((r & -r) != r);
return 0;
}
private:
int _n, size, log;
std::vector<S> d;
std::vector<F> lz;
void update(int k) { d[k] = op(d[2 * k], d[2 * k + 1]); }
void all_apply(int k, F f) {
d[k] = mapping(f, d[k]);
if (k < size) lz[k] = composition(f, lz[k]);
}
void push(int k) {
all_apply(2 * k, lz[k]);
all_apply(2 * k + 1, lz[k]);
lz[k] = id();
}
};
#endif
/// @brief 双対セグメント木 @tparam T セグメント木の型 @tparam composition セグメント木のマージ関数 @tparam id セグメント木の単位元
/// @fn apply 区間に作用を適用する @fn get 位置pの値を取得する
template <class F, F (*composition)(F, F), F (*id)()> struct dual_segtree {
public:
/// @brief セグメント木を初期化する @param n サイズ
explicit dual_segtree(int n) : dual_segtree(std::vector<F>(n, id())) {}
/// @brief セグメント木を初期化する @param v vector<F>型の配列
explicit dual_segtree(const std::vector<F> &v) : _n(int(v.size())) {
log = ceil_pow2(_n);
size = 1 << log;
lz = std::vector<F>(2 * size, id());
for (int i = 0; i < _n; i++) lz[size + i] = v[i];
}
/// @brief [l, r) の要素に f を作用させます O(logN) @param l 左端 @param r 右端 @param f 作用素
void apply(int l, int r, F f) {
assert(0 <= l && l <= r && r <= _n);
if (l == r) return;
l += size;
r += size;
for (int i = log; i >= 1; i--) { // 遅延評価
if (((l >> i) << i) != l) push(l >> i); // lがiの倍数でない場合は、lを親に移動
if (((r >> i) << i) != r) push((r - 1) >> i); // rがiの倍数でない場合は、rを親に移動
}
while (l < r) {
if (l & 1) all_apply(l++, f); // lが奇数の場合は、lに作用を適用してからlをインクリメント
if (r & 1) all_apply(--r, f); // rが奇数の場合は、rをデクリメントしてからrに作用を適用
l >>= 1; // lを親に移動
r >>= 1; // rを親に移動
}
}
/// @brief 位置pの値を取得する @param p 位置
F get(int p) {
assert(0 <= p && p < _n);
p += size;
for (int i = log; i >= 1; i--) push(p >> i);
return lz[p];
}
private:
int _n, size, log;
std::vector<F> lz;
/// @brief 作用素を遅延評価する @param i 位置 @param f 作用素
void all_apply(int i, F f) {
lz[i] = composition(f, lz[i]);
}
/// @brief 作用素を遅延評価する @param i 位置
void push(int i) {
assert(i < size);
all_apply(2 * i, lz[i]);
all_apply(2 * i + 1, lz[i]);
lz[i] = id();
}
};
// ----- segment trees by @rrrrikiOW -----
/**
* @author @rrrrikiOW
* @brief EulerTourを行い、部分木・パスクエリ・LCAを解く
* @tparam E_T 辺の重みの型 @tparam V_T 頂点の重みの型
*/
template <typename E_T,typename V_T> struct EulerTour {
private:
/// @brief 辺の重みを管理
struct PE{
E_T value; // 辺の重み
int rate; // 0 寄与しない 1 加算する -1 減算
PE(E_T value, int rate) : value(value), rate(rate) {}
PE() : value(0), rate(0) {}
E_T val() { return value*rate; }
};
/// @brief 頂点の重みを管理
struct PV{
V_T value; // 頂点の重み
int rate; // 0 寄与しない 1 加算する -1 減算
PV(V_T value, int rate) : value(value), rate(rate) {}
PV() : value(0), rate(0) {}
V_T val() { return value*rate; }
};
bool build_lca_flag = false;
bool build_seg_edge_flag[2] = {false, false};
bool build_seg_vert_flag[2] = {false, false};
vector<V_T> vert_weight; // 頂点の重み
unordered_map<ll,pair<int,int>> edge_to_index;
vector<PE> edge_weight_tour, edge_weight_tour_minus; // 辺の重み
vector<PV> vert_weight_tour, vert_weight_tour_minus; // 頂点の重み
vector<pair<int, int>> depth; // (depth, node)
using S = pair<int, int>;
static S op_lca(S a, S b) { return min(a, b); } // depthの最小値を求める
static S e_lca() { return {1e9,1e9}; } // e: 単位元
segtree<S, &EulerTour::op_lca, &EulerTour::e_lca> seg_lca; // LCAを求めるセグメント木
static PE op_edge(PE a, PE b) { return PE(a.value*a.rate + b.value*b.rate, 1); }
static PE e_edge() { return PE(0, 0); } // e: 単位元
static PE mapping_edge(E_T f, PE x) { return (x.rate == 0) ? x : PE(x.value + f, x.rate); } // mapping: 作用素 f に対して x に作用させる関数
static E_T composition_edge(E_T f, E_T g) { return f+g; } // composition: 作用素 f, g を合成する関数
static E_T id_edge() { return 0; } // id: 作用素の単位元
lazy_segtree<PE, &EulerTour::op_edge, &EulerTour::e_edge, E_T, &EulerTour::mapping_edge, &EulerTour::composition_edge, &EulerTour::id_edge> seg_edge0, seg_edge1; // 辺の合計を管理するセグメント木
static PV op_vert(PV a, PV b) { return PV(a.value*a.rate + b.value*b.rate, 1); }
static PV e_vert() { return PV(0, 0); } // e: 単位元
static PV mapping_vert(V_T f, PV x) { return (x.rate == 0) ? x : PV(x.value + f, x.rate); } // mapping: 作用素 f に対して x に作用させる関数
static V_T composition_vert(V_T f, V_T g) { return f+g; } // composition: 作用素 f, g を合成する関数
static V_T id_vert() { return 0; } // id: 作用素の単位元
lazy_segtree<PV, &EulerTour::op_vert, &EulerTour::e_vert, V_T, &EulerTour::mapping_vert, &EulerTour::composition_vert, &EulerTour::id_vert> seg_vert0, seg_vert1; // 点の合計を管理するセグメント木
/// @brief lcaを構築する
void build_lca() {
seg_lca = segtree<S, op_lca, e_lca>(depth);
build_lca_flag = true;
}
/// @brief seg_edgeを構築する @param ind 0 部分木 or 1 パスクエリ
void build_seg_edge(int ind) {
build_seg_edge_flag[ind] = true;
if (ind == 0) seg_edge0 = lazy_segtree<PE, &EulerTour::op_edge, &EulerTour::e_edge, E_T, &EulerTour::mapping_edge, &EulerTour::composition_edge, &EulerTour::id_edge>(edge_weight_tour);
else seg_edge1 = lazy_segtree<PE, &EulerTour::op_edge, &EulerTour::e_edge, E_T, &EulerTour::mapping_edge, &EulerTour::composition_edge, &EulerTour::id_edge>(edge_weight_tour_minus);
}
/// @brief seg_vertを構築する @param ind 0 部分木 or 1 パスクエリ
void build_seg_vert(int ind) {
build_seg_vert_flag[ind] = true;
if (ind == 0) seg_vert0 = lazy_segtree<PV, &EulerTour::op_vert, &EulerTour::e_vert, V_T, &EulerTour::mapping_vert, &EulerTour::composition_vert, &EulerTour::id_vert>(vert_weight_tour);
else seg_vert1 = lazy_segtree<PV, &EulerTour::op_vert, &EulerTour::e_vert, V_T, &EulerTour::mapping_vert, &EulerTour::composition_vert, &EulerTour::id_vert>(vert_weight_tour_minus);
}
public:
vector<int> in, out;
// コンストラクタ
EulerTour(Graph<E_T> &g, int root=0, vector<V_T> vert_w=vector<V_T>()) : in(g.size()), out(g.size()){
if (vert_w.size() == 0) vert_weight = vector<V_T>(g.size(), 1);
else vert_weight = vert_w;
int idx = 0;
edge_weight_tour.reserve(2 * g.size());
edge_weight_tour_minus.reserve(2 * g.size());
vert_weight_tour.reserve(2 * g.size());
vert_weight_tour_minus.reserve(2 * g.size());
edge_weight_tour.emplace_back(PE(0, 0));
edge_weight_tour_minus.emplace_back(PE(0, 0));
vert_weight_tour.emplace_back(PV(vert_weight[root], 1));
vert_weight_tour_minus.emplace_back(PV(vert_weight[root], 1));
depth.emplace_back(0, root);
function<void(int, int, int)> dfs = [&](int v, int p, int d) {
in[v] = idx++;
for (auto e : g[v]) {
if (e.to == p) continue;
depth.emplace_back(d + 1, e.to);
edge_weight_tour.emplace_back(PE(e.cost, 1));
edge_weight_tour_minus.emplace_back(PE(e.cost, 1));
vert_weight_tour.emplace_back(PV(vert_weight[e.to], 1));
vert_weight_tour_minus.emplace_back(PV(vert_weight[e.to], 1));
edge_to_index[min(v,(int)e.to) *(ll)1e10 + max(v,(int)e.to)] = {idx, -1};
dfs(e.to, v, d + 1);
edge_to_index[min(v,(int)e.to) *(ll)1e10 + max(v,(int)e.to)].second = idx++;
depth.emplace_back(d, v);
edge_weight_tour.emplace_back(PE(0,0));
edge_weight_tour_minus.emplace_back(PE(e.cost, -1));
vert_weight_tour.emplace_back(PV(0, 0));
vert_weight_tour_minus.emplace_back(PV(vert_weight[e.to], -1));
}
out[v] = idx;
};
dfs(root, -1, 0);
}
/// @brief lcaを求める @param u ノードu @param v ノードv @return lca
int lca(int u, int v) {
if (!build_lca_flag) build_lca();
return seg_lca.prod(min(in[u], in[v]), max(out[u], out[v])).second;
}
/// @brief 辺のパスクエリ @param u ノードu @param v ノードv @return 辺のコスト
E_T query_edge(int u, int v) {
// [0, in[u]+1) + [0, in[v]+1) - 2* [0, in[lca(u,v)+1)]
if (!build_seg_edge_flag[1]) build_seg_edge(1);
if (!build_lca_flag) build_lca();
return (seg_edge1.prod(0, in[u]+1)).val() + (seg_edge1.prod(0, in[v]+1)).val() - 2*(seg_edge1.prod(0, in[lca(u,v)]+1)).val();
}
/// @brief 頂点のパスクエリ @param u ノードu @param v ノードv @return 頂点のコスト
V_T query_vert(int u, int v) {
// [0, in[u]+1) + [0, in[v]+1) - 2* [0, in[lca(u,v)+1)]
if (!build_seg_vert_flag[1]) build_seg_vert(1);
if (!build_lca_flag) build_lca();
return (seg_vert1.prod(0, in[u]+1)).val() + (seg_vert1.prod(0, in[v]+1)).val() - 2*(seg_vert1.prod(0, in[lca(u,v)]+1)).val();
}
/// @brief 部分木の辺コストの合計を求める @param u ノードu @return 部分木の辺のコスト
E_T query_subtree_edge(int u) {
if (!build_seg_edge_flag[0]) build_seg_edge(0);
return (seg_edge0.prod(in[u]+1, out[u])).val();
}
/// @brief 部分木の頂点コストの合計を求める @param u ノードu @return 部分木の頂点のコスト
V_T query_subtree_vert(int u) {
if (!build_seg_vert_flag[0]) build_seg_vert(0);
return (seg_vert0.prod(in[u], out[u])).val();
}
/// @brief 辺のコストを更新する @param u ノードu @param v ノードv @param x 更新する値 @param mode 0 部分木 or 1 パスクエリ or 2 両方
void update_edge(int u, int v, E_T x, int mode=1) {
if (u>v) swap(u,v);
if (mode != 1) {
if (!build_seg_edge_flag[0]) build_seg_edge(0);
seg_edge0.set(edge_to_index[(ll)u *1e10 + v].first, PE(x, 1));
}
if (mode != 0) {
if (!build_seg_edge_flag[1]) build_seg_edge(1);
seg_edge1.set(edge_to_index[(ll)u *1e10 + v].first, PE(x, 1));
seg_edge1.set(edge_to_index[(ll)u *1e10 + v].second, PE(x, -1));
}
}
/// @brief 辺のコストを加算する @param u ノードu @param v ノードv @param x 加算する値 @param mode 0 部分木 or 1 パスクエリ or 2 両方
void add_cost_edge(int u, int v, E_T x, int mode=1) {
if (u>v) swap(u,v);
if (mode != 1) {
if (!build_seg_edge_flag[0]) build_seg_edge(0);
PE target = seg_edge0.get(edge_to_index[(ll)u *1e10 + v].first);
seg_edge0.set(edge_to_index[(ll)u *1e10 + v].first, target.value + x, target.rate);
}
if (mode != 0) {
if (!build_seg_edge_flag[1]) build_seg_edge(1);
PE target = seg_edge1.get(edge_to_index[(ll)u *1e10 + v].first);
PE target2 = seg_edge1.get(edge_to_index[(ll)u *1e10 + v].second);
seg_edge1.set(edge_to_index[(ll)u *1e10 + v].first, target.value + x, target.rate);
seg_edge1.set(edge_to_index[(ll)u *1e10 + v].second, target2.value + x, target2.rate);
}
}
/// @brief 頂点のコストを更新する @param u ノードu @param x 更新する値 @param mode 0 部分木 or 1 パスクエリ or 2 両方
void update_vert(int u, V_T x, int mode=1) {
if (mode != 1) {
if (!build_seg_vert_flag[0]) build_seg_vert(0);
seg_vert0.set(in[u], PV(x, 1));
}
if (mode != 0) {
if (!build_seg_vert_flag[1]) build_seg_vert(1);
seg_vert1.set(in[u], PV(x, 1));
seg_vert1.set(out[u], PV(x, -1));
}
}
/// @brief 頂点のコストを加算する @param u ノードu @param x 加算する値 @param mode 0 部分木 or 1 パスクエリ or 2 両方
void add_cost_vert(int u, V_T x, int mode=1) {
if (mode != 1) {
if (!build_seg_vert_flag[0]) build_seg_vert(0);
PV target = seg_vert0.get(in[u]);
seg_vert0.set(in[u], target.value + x, target.rate);
}
if (mode != 0) {
if (!build_seg_vert_flag[1]) build_seg_vert(1);
PV target = seg_vert1.get(in[u]);
PV target2 = seg_vert1.get(out[u]);
seg_vert1.set(in[u], target.value + x, target.rate);
seg_vert1.set(out[u], target2.value + x, target2.rate);
}
}
/// @brief 部分木の全頂点にコストを加算する @param u ノードu @param x 加算する値 @param mode 0 部分木 or 1 パスクエリ or 2 両方
void add_cost_subtree_vert(int u, V_T x, int mode=1) {
if (mode != 1) {
if (!build_seg_vert_flag[0]) build_seg_vert(0);
seg_vert0.apply(in[u], out[u], x);
}
if (mode != 0) {
if (!build_seg_vert_flag[1]) build_seg_vert(1);
seg_vert1.apply(in[u], out[u], x);
}
}
/// @brief 部分木の全辺にコストを加算する @param u ノードu @param x 加算する値 @param mode 0 部分木 or 1 パスクエリ or 2 両方
void add_cost_subtree_edge(int u, E_T x, int mode=1) {
if (mode != 1) {
if (!build_seg_edge_flag[0]) build_seg_edge(0);
seg_edge0.apply(in[u]+1, out[u], x);
}
if (mode != 0) {
if (!build_seg_edge_flag[1]) build_seg_edge(1);
seg_edge1.apply(in[u]+1, out[u], x);
}
}
/// @brief depthを取得する @param u ノードu @return depth
int get_depth(int u) {
return depth[in[u]].first;
}
/// @brief 頂点のコストを取得する @param u ノードu @param mode 0 部分木 or 1 パスクエリ
V_T get_vert(int u, int mode=1) {
if (mode == 0) {
if (!build_seg_vert_flag[0]) build_seg_vert(0);
return seg_vert0.get(in[u]).val();
} else {
if (!build_seg_vert_flag[1]) build_seg_vert(1);
return seg_vert1.get(in[u]).val();
}
}
/// @brief 辺のコストを取得する @param u ノードu @param v ノードv @param mode 0 部分木 or 1 パスクエリ
E_T get_edge(int u, int v, int mode=1) {
if (u>v) swap(u,v);
if (mode == 0) {
if (!build_seg_edge_flag[0]) build_seg_edge(0);
return seg_edge0.get(edge_to_index[(ll)u *1e10 + v].first).val();
} else {
if (!build_seg_edge_flag[1]) build_seg_edge(1);
return seg_edge1.get(edge_to_index[(ll)u *1e10 + v].first).val();
}
}
};
/// @brief 乱択平衡二分木(基底クラス)
template <typename Node>
struct RBSTBase {
using Ptr = Node *;
template <typename... Args>
inline Ptr my_new(Args... args) {
return new Node(args...);
}
inline void my_del(Ptr t) { delete t; }
/// @brief 空の木を返す。 O(1)
inline Ptr make_tree() const { return nullptr; }
// for avoiding memory leak, activate below
/*
using Ptr = shared_ptr<Node>;
template <typename... Args>
inline Ptr my_new(Args... args) {
return make_shared<Node>(args...);
}
inline void my_del(Ptr t) {}
Ptr make_tree() {return Ptr();}
*/
/// @brief tの大きさを返す。 O(1)
int size(Ptr t) const { return count(t); }
/// @brief lとrをマージして新たに出来た木のポインタを返す。 O(logN)
Ptr merge(Ptr l, Ptr r) {
if (!l || !r) return l ? l : r;
if (int((rng() * (l->cnt + r->cnt)) >> 32) < l->cnt) {
push(l);
l->r = merge(l->r, r);
return update(l);
} else {
push(r);
r->l = merge(l, r->l);
return update(r);
}
}
///@brief tを[0, k)と[k, |t|)の二つの木に分割する。 O(logN)
pair<Ptr, Ptr> split(Ptr t, int k) {
if (!t) return {nullptr, nullptr};
push(t);
if (k <= count(t->l)) {
auto s = split(t->l, k);
t->l = s.second;
return {s.first, update(t)};
} else {
auto s = split(t->r, k - count(t->l) - 1);
t->r = s.first;
return {update(t), s.second};
}
}
Ptr build(int l, int r, const vector<decltype(Node::key)> &v) {
if (l + 1 == r) return my_new(v[l]);
int m = (l + r) >> 1;
Ptr pm = my_new(v[m]);
if (l < m) pm->l = build(l, m, v);
if (m + 1 < r) pm->r = build(m + 1, r, v);
return update(pm);
}
/// @brief 列vを初期値とした新たな木を作る。 O(N)
Ptr build(const vector<decltype(Node::key)> &v) {
return build(0, (int)v.size(), v);
}
/// @brief tのk番目にNode(args...)を挿入する。 O(logN)
template <typename... Args>
void insert(Ptr &t, int k, const Args &... args) {
auto x = split(t, k);
t = merge(merge(x.first, my_new(args...)), x.second);
}
/// @brief tのk番目のノードを削除する。 O(logN)
void erase(Ptr &t, int k) {
auto x = split(t, k);
auto y = split(x.second, 1);
my_del(y.first);
t = merge(x.first, y.second);
}
protected:
static uint64_t rng() {
static uint64_t x_ = 88172645463325252ULL;
return x_ ^= x_ << 7, x_ ^= x_ >> 9, x_ & 0xFFFFFFFFull;
}
inline int count(const Ptr t) const { return t ? t->cnt : 0; }
virtual void push(Ptr) = 0;
virtual Ptr update(Ptr) = 0;
};
/// @brief RBSTのノード キー（値）、部分木の合計、遅延伝搬用の変数、左右の子ノードへのポインタ、部分木のサイズ、反転フラグ
template <typename T, typename E>
struct LazyReversibleRBSTNode {
typename RBSTBase<LazyReversibleRBSTNode>::Ptr l, r;
T key, sum;
E lazy;
int cnt;
bool rev;
LazyReversibleRBSTNode(const T &t = T(), const E &e = E())
: l(), r(), key(t), sum(t), lazy(e), cnt(1), rev(false) {}
};
/// @brief 遅延伝搬反転可能乱択平衡二分木
/// @see https://nyaannyaan.github.io/library/rbst/lazy-reversible-rbst.hpp.html
template <typename T, typename E, T (*f)(T, T), T (*g)(T, E), E (*h)(E, E),
T (*ts)(T)>
struct LazyReversibleRBST : RBSTBase<LazyReversibleRBSTNode<T, E>> {
using Node = LazyReversibleRBSTNode<T, E>;
using base = RBSTBase<LazyReversibleRBSTNode<T, E>>;
using base::merge;
using base::split;
using typename base::Ptr;
LazyReversibleRBST() = default;
/// @brief tを反転する。 O(1)
void toggle(Ptr t) {
if(!t) return;
swap(t->l, t->r);
t->sum = ts(t->sum);
t->rev ^= true;
}
/// @brief fold(t, a, b) : [a, b) の和を求める O(logN)
T fold(Ptr &t, int a, int b) {
auto x = split(t, a);
auto y = split(x.second, b - a);
auto ret = sum(y.first);
t = merge(x.first, merge(y.first, y.second));
return ret;
}
/// @brief reverse(t, a, b) : [a, b) を反転する O(logN)
void reverse(Ptr &t, int a, int b) {
auto x = split(t, a);
auto y = split(x.second, b - a);
toggle(y.first);
t = merge(x.first, merge(y.first, y.second));
}
/// @brief apply(t, a, b, e) : [a, b) に e を作用させる O(logN)
void apply(Ptr &t, int a, int b, const E &e) {
auto x = split(t, a);
auto y = split(x.second, b - a);
propagate(y.first, e);
t = merge(x.first, merge(y.first, y.second));
}
protected:
inline T sum(const Ptr t) const { return t ? t->sum : T(); }
Ptr update(Ptr t) override {
push(t);
t->cnt = 1;
t->sum = t->key;
if (t->l) t->cnt += t->l->cnt, t->sum = f(t->l->sum, t->sum);
if (t->r) t->cnt += t->r->cnt, t->sum = f(t->sum, t->r->sum);
return t;
}
void push(Ptr t) override {
if (t->rev) {
if (t->l) toggle(t->l);
if (t->r) toggle(t->r);
t->rev = false;
}
if (t->lazy != E()) {
if (t->l) propagate(t->l, t->lazy);
if (t->r) propagate(t->r, t->lazy);
t->lazy = E();
}
}
void propagate(Ptr t, const E &x) {
t->lazy = h(t->lazy, x);
t->key = g(t->key, x);
t->sum = g(t->sum, x);
}
};
template <typename T> struct BipartiteGraph : dsu {
vector<int> color; // 0 or 1
BipartiteGraph(int n) : color(n + n, -1), dsu(n + n) {} // n: 頂点数
/// @brief 二部グラフ判定 @param g グラフ @return 二部グラフかどうか
bool bipartite(Graph<T> &g) {
for (int i = 0; i < g.size(); i++) {
for (auto e : g[i]) {
if (same(i, e.to)) return false;
merge(e.from, e.to + g.size());
merge(e.to, e.from + g.size());
}
}
return true;
}
/// @brief 二部グラフ彩色 @return 二部グラフ彩色が可能かどうか
bool bipartite_graph_coloring() {
for (int i = 0; i < color.size() / 2; i++) {
int a = leader(i), b = leader(i + color.size() / 2);
if (a == b) return false;
if (color[a] == -1) color[a] = 0, color[b] = 1;
}
return true;
}
int operator[](int i) { return color[leader(i)]; }
};
// clang-format on
#endif
/*
******* 神龜雖壽 *******
******* 猶有竟時 *******
******* 騰蛇乘霧 *******
******* 終爲土灰 *******
******* 老驥伏櫪 *******
******* 志在千里 *******
******* 烈士暮年 *******
******* 壯心不已 *******
******* 盈縮之期 *******
******* 不但在天 *******
******* 養怡之福 *******
******* 可得永年 *******
******* 幸甚至哉 *******
******* 歌以詠志 *******
*/