1/*2 * Copyright (c) 2026 Bruno Monteiro3 *4 * Permission is hereby granted, free of charge, to any person obtaining a copy5 * of this software and associated documentation files (the "Software"), to deal6 * in the Software without restriction, including without limitation the rights7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell8 * copies of the Software, and to permit persons to whom the Software is9 * furnished to do so, subject to the following conditions:10 *11 * The above copyright notice and this permission notice shall be included in12 * all copies or substantial portions of the Software.13 *14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN20 * THE SOFTWARE.21 */2223#pragma once2425#include <algorithm>26#include <bitset>27#include <cstdint>28#include <functional>29#include <initializer_list>30#include <iomanip>31#include <iostream>32#include <limits>33#include <map>34#include <optional>35#include <queue>36#include <random>37#include <set>38#include <sstream>39#include <stdexcept>40#include <string>41#include <sys/types.h>42#include <type_traits>43#include <unordered_map>44#include <unordered_set>45#include <utility>46#include <vector>4748namespace tgen {4950/**************************51 * *52 * GENERAL OPERATIONS *53 * *54 **************************/5556namespace detail {5758// Type aliases.59using u128 = unsigned __int128;60using i128 = __int128;6162/*63 * Error handling.64 */6566inline void throw_assertion_error(const std::string &condition,67 const std::string &msg, const char *file,68 int line) {69 throw std::runtime_error("tgen: " + msg + " (assertion `" + condition +70 "` failed at " + file + ":" +71 std::to_string(line) + ")");72}73inline void throw_assertion_error(const std::string &condition,74 const char *file, int line) {75 throw std::runtime_error("tgen: assertion `" + condition + "` failed at " +76 std::string(file) + ":" + std::to_string(line));77}78inline std::runtime_error error(const std::string &msg) {79 return std::runtime_error("tgen: " + msg);80}81inline std::runtime_error contradiction_error(const std::string &type,82 const std::string &msg = "") {83 // Tried to generate a contradictory type.84 std::string error_msg =85 type + ": invalid " + type + " (contradictory restrictions)";86 if (!msg.empty())87 error_msg += ": " + msg;88 return error(error_msg);89}90inline std::runtime_error91complex_restrictions_error(const std::string &type,92 const std::string &msg = "") {93 // Tried to generate a type with too many distinct restrictions.94 std::string error_msg =95 type + ": cannot represent " + type + " (complex restrictions)";96 if (!msg.empty())97 error_msg += ": " + msg;98 return error(error_msg);99}100inline void tgen_ensure_against_bug(bool cond, const std::string &msg = "") {101 if (!cond) {102 std::string error_msg;103 if (!msg.empty())104 error_msg = "tgen: " + msg + "\n";105 error_msg += "tgen: THERE IS A BUG IN TGEN; PLEASE CONTACT MAINTAINERS";106 throw std::runtime_error(error_msg);107 }108}109110// Ensures condition is true, with a clear error message on failure.111#define tgen_ensure(cond, ...) \112 if (!(cond)) \113 tgen::detail::throw_assertion_error(#cond, ##__VA_ARGS__, __FILE__, \114 __LINE__)115116// Registering checks.117inline bool registered = false;118inline void ensure_registered() {119 tgen_ensure(registered,120 "tgen was not registered! You should call "121 "tgen::register_gen(argc, argv) before running tgen functions");122}123124// Template magic to detect types at compile time.125126// Detects containers != std::string.127template <typename T, typename = void> struct is_container : std::false_type {};128template <typename T>129struct is_container<T,130 std::void_t<typename std::remove_reference_t<T>::value_type,131 decltype(std::begin(std::declval<T>())),132 decltype(std::end(std::declval<T>()))>>133 : std::true_type {};134// Exclude all basic_string variants135template <typename Char, typename Traits, typename Alloc>136struct is_container<std::basic_string<Char, Traits, Alloc>> : std::false_type {137};138template <typename Char, typename Traits, typename Alloc>139struct is_container<const std::basic_string<Char, Traits, Alloc>>140 : std::false_type {};141template <typename Char, typename Traits, typename Alloc>142struct is_container<std::basic_string<Char, Traits, Alloc> &>143 : std::false_type {};144template <typename Char, typename Traits, typename Alloc>145struct is_container<const std::basic_string<Char, Traits, Alloc> &>146 : std::false_type {};147148// Detects std::pair.149template <typename T> struct is_pair : std::false_type {};150template <typename A, typename B>151struct is_pair<std::pair<A, B>> : std::true_type {};152// Detects std::tuple.153template <typename T> struct is_tuple : std::false_type {};154template <typename... Ts>155struct is_tuple<std::tuple<Ts...>> : std::true_type {};156// Detects scalar (printed atomically).157template <typename T>158struct is_scalar159 : std::bool_constant<!is_container<T>::value and !is_tuple<T>::value and160 !is_pair<T>::value> {};161// Detects complex container.162template <typename T>163struct is_container_multiline164 : std::bool_constant<is_container<T>::value and165 !is_scalar<typename std::remove_cv_t<166 std::remove_reference_t<T>>::value_type>::value> {167};168// Detects complex std::pair.169template <typename T> struct is_pair_multiline : std::false_type {};170template <typename A, typename B>171struct is_pair_multiline<std::pair<A, B>>172 : std::bool_constant<!is_scalar<A>::value or !is_scalar<B>::value> {};173// Detects complex std::tuple.174template <typename Tuple> struct is_tuple_multiline : std::false_type {};175template <typename... Ts>176struct is_tuple_multiline<std::tuple<Ts...>>177 : std::bool_constant<(!is_scalar<Ts>::value or ...)> {};178179// Used to return false at compile time only if evaluated.180template <typename> inline constexpr bool dependent_false_v = false;181182/*183 * Properties of custom types.184 */185186// If type is sequential (list-like).187using is_sequential_tag = void;188189// Detects associative containers.190template <typename T, typename = void>191struct is_associative_container : std::false_type {};192template <typename T>193struct is_associative_container<194 T, std::void_t<typename T::key_type, typename T::key_compare>>195 : std::true_type {};196197// Detects sequential generator values.198template <typename T, typename = void>199struct is_sequential : std::false_type {};200template <typename T>201struct is_sequential<202 T, std::void_t<typename std::decay_t<T>::tgen_is_sequential_tag>>203 : std::true_type {};204205/*206 * Unique rng to use.207 */208209// The single rng to be used by the library.210inline std::mt19937 rng;211212/*213 * Printing.214 */215216// Print view struct for printing either a container or a sequential generator217// element.218template <typename T,219 bool IsCont = detail::is_container<std::decay_t<T>>::value>220struct print_cols_view;221222// Container.223template <typename T> struct print_cols_view<T, true> {224 const T &value;225 decltype(std::begin(std::declval<const T &>())) it;226227 print_cols_view(const T &v) : value(v), it(v.begin()) {}228229 std::size_t size() const { return value.size(); }230 decltype(auto) get(std::size_t) const { return *it; }231 void advance() { ++it; }232};233234// Sequential generator element.235template <typename T> struct print_cols_view<T, false> {236 const T &value;237238 print_cols_view(const T &v) : value(v) {}239240 std::size_t size() const { return value.size(); }241 decltype(auto) get(std::size_t i) const { return value[i]; }242 void advance() {}243};244245/*246 * Distinct generation.247 */248249// Rejection cap is multiplier * |seen|; with one value left, falsely reporting250// exhaustion has probability about e^{-84} < 10^{-36}.251constexpr int distinct_attempt_multiplier = 84;252253// One rejection-sampling step for distinct generation.254// O(T * log k + log^2 k) amortized expected time per call when generating k255// distinct values and next() runs in O(T).256template <typename Seen, typename Fn>257auto try_generate_distinct(Seen &seen, Fn &&next, bool insert = true)258 -> std::optional<std::invoke_result_t<Fn &>> {259 using T = std::invoke_result_t<Fn &>;260 size_t attempts =261 distinct_attempt_multiplier * std::max<size_t>(1, seen.size());262 for (size_t i = 0; i < attempts; ++i) {263 T val = next();264 if (insert) {265 if (seen.insert(val).second)266 return val;267 } else if (seen.count(val) == 0)268 return val;269 }270 return std::nullopt;271}272273} // namespace detail274275/*276 * Compiler configuration (see set_compiler).277 */278279// Kinds of compilers.280enum class compiler_kind { gcc, clang, unknown };281282// Compiler identity and version.283struct compiler_value {284 compiler_kind kind_;285 int major_;286 int minor_;287288 compiler_value(compiler_kind kind = compiler_kind::unknown, int major = 0,289 int minor = 0)290 : kind_(kind), major_(major), minor_(minor) {}291};292293namespace detail {294295// Global C++ version value (0 means unknown).296struct cpp_value {297 int version_;298299 cpp_value(std::optional<int> version = std::nullopt)300 : version_(version ? *version : 0) {301 if (version) {302 tgen_ensure(*version == 17 or *version == 20 or *version == 23,303 "unsupported C++ version (use 17, 20, 23)");304 }305 }306};307308inline cpp_value cpp;309inline compiler_value compiler;310311} // namespace detail312313/*314 * Base classes.315 */316317// Needed for return type of some functions.318template <typename T> struct list;319320// Generates distinct values of a function.321template <typename Func, typename... Args> struct distinct {322 Func func_;323 std::tuple<Args...> args_;324 using T = std::invoke_result_t<Func &, Args &...>;325 std::set<T> seen_;326327 distinct(Func func, Args... args)328 : func_(std::move(func)), args_(std::move(args)...) {}329330 // Generates a distinct value (i.e., one not returned before).331 //332 // Assume gen() produces a uniformly random value in O(T) time.333 // Since duplicates are rejected, the expected number of trials over334 // k successful generations is:335 //336 // sum_{i=1}^k k / i = O(k log k)337 //338 // (coupon collector argument).339 //340 // Each trial additionally performs O(log k) work to check/store341 // previously generated values, yielding a total time of342 // O((T + log k) * k log k).343 //344 // Thus, the amortized expected time per call is345 // O(T * log k + log^2 k).346 //347 // With extremely small probability (< 1e-18), the algorithm may348 // incorrectly report that no more distinct values exist.349 auto gen() {350 auto val = generate_distinct(true);351 if (val)352 return *val;353354 throw detail::error("distinct: no more distinct values");355 }356 template <typename U> auto gen(std::initializer_list<U> il) {357 return gen(std::vector<U>(il));358 }359360 // Generates a list of distinct values.361 auto gen_list(int size) {362 std::vector<T> res;363 for (int i = 0; i < size; ++i)364 res.push_back(gen());365366 return typename list<T>::value(res);367 }368369 // Checks if there are no more distinct values.370 // With extremely small probability (< 1e-18), the algorithm may371 // incorrectly report that there are no more distinct values.372 bool empty() { return generate_distinct(false) == std::nullopt; }373374 // Generates all distinct values.375 auto gen_all() {376 std::vector<T> res;377 while (true) {378 auto val = generate_distinct(true);379 if (val)380 res.push_back(*val);381 else382 break;383 }384 return typename list<T>::value(res);385 }386387 // Nice error for `out << distinct`.388 friend std::ostream &operator<<(std::ostream &out, const distinct &) {389 static_assert(390 detail::dependent_false_v<distinct>,391 "distinct: cannot print a distinct generator. Maybe you forgot to "392 "call `gen()`?");393 return out;394 }395396 private:397 // Generates distinct value and inserts it if `insert` is true.398 // Returns the value if found, otherwise returns std::nullopt.399 auto generate_distinct(bool insert) {400 return detail::try_generate_distinct(401 seen_, [&] { return std::apply(func_, args_); }, insert);402 }403};404template <typename Func, typename... Args>405distinct(Func, Args...) -> distinct<Func, Args...>;406407// Base struct for generators.408template <typename Gen> struct gen_base {409 const Gen &self() const { return *static_cast<const Gen *>(this); }410411 template <typename... Args> auto gen_list(int size, Args &&...args) const {412 std::vector<typename Gen::value> res;413414 for (int i = 0; i < size; ++i)415 res.push_back(static_cast<const Gen *>(this)->gen(416 std::forward<Args>(args)...));417418 return typename list<typename Gen::value>::value(res);419 }420421 // Calls the generator until predicate is true.422 template <typename Pred, typename... Args>423 auto gen_until(Pred predicate, int max_tries, Args &&...args) const {424 for (int i = 0; i < max_tries; ++i) {425 typename Gen::value val = static_cast<const Gen *>(this)->gen(426 std::forward<Args>(args)...);427428 if (predicate(val))429 return val;430 }431432 throw detail::error("could not generate value matching predicate");433 }434 template <typename Pred, typename T, typename... Args>435 auto gen_until(Pred predicate, int max_tries, std::initializer_list<T> il,436 Args &&...args) const {437 return gen_until(predicate, max_tries, std::vector<T>(il),438 std::forward<Args>(args)...);439 }440441 // Distinct for generator.442 template <typename... Args> auto distinct(Args &&...args) const {443 return tgen::distinct(444 [self = self()](auto &&...inner_args) mutable -> decltype(auto) {445 return self.gen(446 std::forward<decltype(inner_args)>(inner_args)...);447 },448 std::forward<Args>(args)...);449 }450 template <typename T, typename... Args>451 auto distinct(std::initializer_list<T> il, Args &&...args) const {452 return distinct(std::vector<T>(il), std::forward<Args>(args)...);453 }454455 // Nice error for `out << generator`.456 friend std::ostream &operator<<(std::ostream &out, const gen_base &) {457 static_assert(detail::dependent_false_v<gen_base>,458 "gen_base: cannot print a generator. Maybe you forgot to "459 "call `gen()`?");460 return out;461 }462};463464// Base class for generator values.465template <typename Val> struct gen_value_base {466 const Val &self() const { return *static_cast<const Val *>(this); }467468 bool operator<(const Val &rhs) const {469 return self().to_std() < rhs.to_std();470 }471};472473namespace detail {474475// Detects generator values.476template <typename T>477struct is_generator_value478 : std::is_base_of<gen_value_base<std::decay_t<T>>, std::decay_t<T>> {};479480} // namespace detail481482/*483 * Easier printing.484 */485486// Struct to print standard types to std::ostream;487struct print {488 std::string s_;489490 template <typename T> print(const T &val, char sep = ' ') {491 std::ostringstream oss;492 write(oss, val, sep);493 s_ = oss.str();494 }495 template <typename T>496 print(const std::initializer_list<T> &il, char sep = ' ') {497 std::ostringstream oss;498 write(oss, std::vector<T>(il), sep);499 s_ = oss.str();500 }501 template <typename T>502 print(const std::initializer_list<std::initializer_list<T>> &il,503 char sep = ' ') {504 std::ostringstream oss;505 std::vector<std::vector<T>> mat;506 for (const auto &i : il)507 mat.push_back(i);508 write(oss, mat, sep);509 s_ = oss.str();510 }511512 template <typename T> void write(std::ostream &os, const T &val, char sep) {513 if constexpr (detail::is_pair<T>::value) {514 if constexpr (detail::is_pair_multiline<T>::value) {515 write(os, val.first, sep);516 os << '\n';517 write(os, val.second, sep);518 } else {519 // Use space for inner separator.520 write(os, val.first, ' ');521 os << sep;522 write(os, val.second, ' ');523 }524 } else if constexpr (detail::is_tuple<T>::value)525 write_tuple(os, val, sep);526 else if constexpr (detail::is_container<T>::value)527 write_container(os, val, sep);528 else if constexpr (std::is_same_v<T, detail::i128> or529 std::is_same_v<T, detail::u128>)530 write_128_number(os, val);531 else532 os << val;533 }534535 // Writes 128 bit number.536 template <typename T> void write_128_number(std::ostream &os, T num) {537 static const long long BASE = 1e18;538539 if (num < 0) {540 os << '-';541 num = -num;542 }543544 if (num >= BASE) {545 write_128_number(os, num / BASE);546 os << std::setw(18) << std::setfill('0')547 << static_cast<long long>(num % BASE);548 } else549 os << static_cast<long long>(num);550 }551 // Writes container, checking separator.552 template <typename C>553 void write_container(std::ostream &os, const C &container, char sep) {554 bool first = true;555556 for (const auto &e : container) {557 if (!first)558 os << (detail::is_container_multiline<C>::value ? '\n' : sep);559 first = false;560 write(os, e, detail::is_container_multiline<C>::value ? sep : ' ');561 }562 }563564 // Writes tuple, checking separator.565 template <typename Tuple, size_t... I>566 void write_tuple_impl(std::ostream &os, const Tuple &tp, char sep,567 std::index_sequence<I...>) {568 bool first = true;569 ((os << (first ? (first = false, "")570 : (detail::is_tuple_multiline<Tuple>::value571 ? "\n"572 : std::string(1, sep))),573 write(os, std::get<I>(tp),574 detail::is_tuple_multiline<Tuple>::value ? sep : ' ')),575 ...);576 }577 template <typename T>578 void write_tuple(std::ostream &os, const T &tp, char sep) {579 write_tuple_impl(os, tp, sep,580 std::make_index_sequence<std::tuple_size<T>::value>{});581 }582583 friend std::ostream &operator<<(std::ostream &out, const print &pr) {584 return out << pr.s_;585 }586};587588// Prints in a new line.589struct println : print {590 template <typename T>591 println(const T &val, char sep = ' ') : print(val, sep) {}592 template <typename T>593 println(const std::initializer_list<T> &il, char sep = ' ')594 : print(il, sep) {}595 template <typename T>596 println(const std::initializer_list<std::initializer_list<T>> &il,597 char sep = ' ')598 : print(il, sep) {}599600 friend std::ostream &operator<<(std::ostream &out, const println &pr) {601 return out << pr.s_ << '\n';602 }603};604605// Prints container / sequential generator value on its own column.606// Example:607// A = {1, 2, 3}, B = {4, 2, 5}608// print_each(A, B) will print:609// "1 4610// 2 2611// 3 5612//",613// that is, it prints the end of the line for all lines.614template <typename... Args> struct print_cols {615 std::string s_;616617 print_cols(const Args &...args) {618 static_assert(619 ((detail::is_container<std::decay_t<Args>>::value or620 detail::is_sequential<std::decay_t<Args>>::value) and621 ...),622 "print_cols: arguments must be containers or sequential generator "623 "values");624 std::ostringstream oss;625 write(oss, args...);626 s_ = oss.str();627 }628629 void write(std::ostream &os, const Args &...args) {630 auto views = std::apply(631 [](const Args &...inner_args) {632 return std::make_tuple(633 detail::print_cols_view<decltype(inner_args)>{634 inner_args}...);635 },636 std::forward_as_tuple(args...));637638 const std::size_t n = std::get<0>(views).size();639640 auto check = [&](const auto &v) {641 tgen_ensure(v.size() == n, "print_cols: sizes should be the same");642 };643 std::apply([&](const auto &...v) { (check(v), ...); }, views);644645 for (std::size_t i = 0; i < n; ++i) {646 bool first = true;647648 std::apply(649 [&](const auto &...v) {650 ((os << (first ? "" : " ") << print(v.get(i)),651 first = false),652 ...);653 },654 views);655656 os << '\n';657658 std::apply([](auto &...v) { (v.advance(), ...); }, views);659 }660 }661662 friend std::ostream &operator<<(std::ostream &out, const print_cols &pr) {663 return out << pr.s_;664 }665};666667/*668 * Global random operations.669 */670671namespace detail {672673// libstdc++ accepts std::uniform_int_distribution with narrow integral types674// (char/signed char/unsigned char/short/bool), but libc++ rejects them with a675// hard static_assert ("IntType must be a supported integer type"). Promote such676// types to a width the standard guarantees, preserving signedness, so the same677// `next<T>` works across both standard libraries (e.g. Apple clang / libc++).678template <typename T>679using uniform_int_t = std::conditional_t<680 (sizeof(T) >= sizeof(short)), T,681 std::conditional_t<std::is_signed_v<T>, int, unsigned int>>;682683} // namespace detail684685// Returns a uniformly random number in [0, right)686// O(1).687template <typename T> T next(T right) {688 detail::ensure_registered();689 if constexpr (std::is_integral_v<T>) {690 tgen_ensure(right >= 1, "value for `next` must be valid");691 return static_cast<T>(692 std::uniform_int_distribution<detail::uniform_int_t<T>>(693 0,694 static_cast<detail::uniform_int_t<T>>(right) - 1)(detail::rng));695 } else if constexpr (std::is_floating_point_v<T>) {696 tgen_ensure(right >= 0, "value for `next` must be valid");697 return std::uniform_real_distribution<T>(0, right)(detail::rng);698 } else699 throw detail::error("invalid type for next (" +700 std::string(typeid(T).name()) + ")");701}702703// Returns a uniformly random number in [left, right].704// For floating-point types, uses uniform_real_distribution ([left, right) in705// C++), equivalent to [left, right] because the right endpoint has probability706// zero.707// O(1).708template <typename T> T next(T left, T right) {709 detail::ensure_registered();710 tgen_ensure(left <= right, "range for `next` must be valid");711 if constexpr (std::is_integral_v<T>)712 return static_cast<T>(713 std::uniform_int_distribution<detail::uniform_int_t<T>>(714 static_cast<detail::uniform_int_t<T>>(left),715 static_cast<detail::uniform_int_t<T>>(right))(detail::rng));716 else if constexpr (std::is_floating_point_v<T>)717 return std::uniform_real_distribution<T>(left, right)(detail::rng);718 else719 throw detail::error("invalid type for next (" +720 std::string(typeid(T).name()) + ")");721}722723// Skewed next.724//725// Returns a random number in [0, right) with a bias controlled by `w`.726// - w = 0:727// Uniform distribution.728// - w > 0:729// Returns the maximum of (w + 1) independent uniform samples.730// Biases the distribution toward larger values.731// The resulting density is proportional to:732// f(x) = x^w733// In particular:734// w = 1 -> linear bias735// w = 2 -> quadratic bias736// w = 3 -> cubic bias737// - w < 0:738// Returns the minimum of (-w + 1) independent uniform samples.739// Symmetric to the w > 0 case.740// The continuous version corresponds to Beta distributions:741// w > 0 -> Beta(w + 1, 1)742// w < 0 -> Beta(1, -w + 1)743// For |w| > 5, the distribution is approximate.744// O(1).745template <typename T> T wnext(T right, int w) {746 // For small |w|, use the naive approach.747 if (abs(w) <= 5) {748 T val = next<T>(right);749 for (int i = 0; i < w; ++i)750 val = std::max(val, next<T>(right));751 for (int i = 0; i < -w; ++i)752 val = std::min(val, next<T>(right));753 return val;754 }755756 // O(1) way.757 double x, r = next<double>(0, 1);758759 if (w >= 0) {760 x = std::pow(r, 1.0 / (w + 1));761 } else {762 x = 1.0 - std::pow(r, 1.0 / (-w + 1));763 }764765 return T(x * right);766}767768// Returns a random number in [left, right] with a bias controlled by `w`.769// O(1).770template <typename T> T wnext(T left, T right, int w) {771 // For small |w|, use the naive approach.772 if (abs(w) <= 5) {773 T val = next<T>(left, right);774 for (int i = 0; i < w; ++i)775 val = std::max(val, next<T>(left, right));776 for (int i = 0; i < -w; ++i)777 val = std::min(val, next<T>(left, right));778 return val;779 }780781 // O(1) way.782 double x, r = next<double>(0, 1);783784 if (w >= 0) {785 x = std::pow(r, 1.0 / (w + 1));786 } else {787 x = 1.0 - std::pow(r, 1.0 / (-w + 1));788 }789790 return left + T(x * (right - left));791}792793namespace detail {794795// Uniformly random 128 bit number in [0, total).796// O(1) expected.797inline u128 next128(u128 total) {798 tgen_ensure(total > 0, "next128: total must be positive");799800 // Largest multiple of total less than 2^128.801 u128 limit = (u128(-1) / total) * total;802803 while (true) {804 // Generate uniform 128-bit random number.805 u128 r = (u128(next<uint64_t>(0, std::numeric_limits<uint64_t>::max()))806 << 64) |807 next<uint64_t>(0, std::numeric_limits<uint64_t>::max());808809 if (r < limit)810 return r % total;811 }812}813814} // namespace detail815816// Weighted sampler.817//818// Generates indices with probability proportional to `distribution`, using819// alias method.820//821// Internally, integral weights are accumulated in unsigned __int128 (exact);822// floating-point weights are accumulated in double.823// <O(n), O(1)>.824template <typename T> struct weighted_sampler {825 static_assert(std::is_arithmetic_v<T>,826 "weighted_sampler requires an arithmetic weight type");827828 // Internal storage type: `u128` for integral inputs (exact arithmetic),829 // `double` for floating-point inputs.830 using storage_t =831 std::conditional_t<std::is_integral_v<T>, detail::u128, double>;832833 int n_;834 std::vector<storage_t> weight_;835 std::vector<int> alias_;836 storage_t total_;837838 // Creates an alias method for generating indices with probability839 // proportional to the distribution.840 // O(n).841 weighted_sampler(const std::vector<T> &distribution)842 : n_(distribution.size()), alias_(n_) {843 tgen_ensure(distribution.size() > 0,844 "weighted_sampler: distribution must be non-empty");845 for (const auto &w : distribution)846 tgen_ensure(w >= 0,847 "weighted_sampler: distribution must be non-negative");848849 total_ = std::accumulate(distribution.begin(), distribution.end(),850 storage_t(0));851852 std::queue<int> big, small;853 for (int i = 0; i < n_; ++i) {854 weight_.push_back(storage_t(n_) * storage_t(distribution[i]));855 if (weight_[i] < total_)856 small.push(i);857 else858 big.push(i);859 }860861 while (!small.empty() and !big.empty()) {862 int s = small.front();863 small.pop();864 int b = big.front();865 big.pop();866867 alias_[s] = b;868869 weight_[b] -= total_ - weight_[s];870 if (weight_[b] < total_)871 small.push(b);872 else873 big.push(b);874 }875876 detail::tgen_ensure_against_bug(877 small.empty(), "weighted_sampler: small must be empty");878879 // The remaining elements should have weight equal to total and be880 // assigned to themselves.881 while (!big.empty()) {882 int b = big.front();883 big.pop();884 if constexpr (std::is_integral_v<T>) {885 detail::tgen_ensure_against_bug(886 weight_[b] == total_,887 "weighted_sampler: weight of big element must be total");888 }889 alias_[b] = b;890 }891 }892 weighted_sampler(const std::initializer_list<T> &distribution)893 : weighted_sampler(std::vector<T>(distribution)) {}894895 // Uniformly random value in [0, total). Overloaded so next() can dispatch896 // at compile time to the right primitive for the chosen `storage_t`.897 static detail::u128 sample_below(detail::u128 total) {898 return detail::next128(total);899 }900 static double sample_below(double total) {901 return tgen::next<double>(0, total);902 }903904 // Generates a random index with probability proportional to the905 // distribution.906 // O(1).907 size_t next() const {908 int i = tgen::next<int>(0, n_ - 1);909 return sample_below(total_) < weight_[i] ? i : alias_[i];910 }911};912template <typename T>913weighted_sampler(const std::vector<T> &) -> weighted_sampler<T>;914template <typename T>915weighted_sampler(const std::initializer_list<T> &) -> weighted_sampler<T>;916917// Returns i with probability proportional to distribution[i].918// O(|distribution|).919template <typename T>920size_t next_by_distribution(const std::vector<T> &distribution) {921 return weighted_sampler(distribution).next();922}923template <typename T>924size_t next_by_distribution(const std::initializer_list<T> &distribution) {925 return next_by_distribution(std::vector<T>(distribution));926}927928// Returns a vector of k indices with probability proportional to929// `distribution`. Uses alias method.930// O(k + |distribution|).931template <typename T>932std::vector<int> many_by_distribution(int k,933 const std::vector<T> &distribution) {934 tgen_ensure(distribution.size() > 0, "distribution must be non-empty");935 tgen_ensure(k >= 0, "number of elements to choose must be non-negative");936937 weighted_sampler am(distribution);938 std::vector<int> res;939 for (int i = 0; i < k; ++i)940 res.push_back(am.next());941 return res;942}943template <typename T>944std::vector<int>945many_by_distribution(int k, const std::initializer_list<T> &distribution) {946 return many_by_distribution(k, std::vector<T>(distribution));947}948949// Shuffles [first, last) inplace uniformly, for RandomAccessIterator.950// O(|container|).951template <typename It> void shuffle(It first, It last) {952 if (first == last)953 return;954955 for (It i = first + 1; i != last; ++i)956 std::iter_swap(i, first + next(0, static_cast<int>(i - first)));957}958959// Shuffles container uniformly.960// O(|container|).961template <typename C> [[nodiscard]] auto shuffled(const C &container) {962 if constexpr (detail::is_associative_container<C>::value) {963 std::vector<typename C::value_type> vec(container.begin(),964 container.end());965 shuffle(vec.begin(), vec.end());966 return vec;967 } else {968 auto new_container = container;969 shuffle(new_container.begin(), new_container.end());970 return new_container;971 }972}973template <typename T>974[[nodiscard]] std::vector<T> shuffled(const std::initializer_list<T> &il) {975 return shuffled(std::vector<T>(il));976}977978// Returns a random element from [first, last) uniformly.979// O(1) for random_access_iterator, O(|last - first|) otherwise.980template <typename It> typename It::value_type pick(It first, It last) {981 int size = std::distance(first, last);982 tgen_ensure(size > 0, "cannot pick from empty range");983 It it = first;984 std::advance(it, next(0, size - 1));985 return *it;986}987988// Returns a random element from container uniformly.989// O(1) for random_access_iterator, O(|container|) otherwise.990template <typename C> typename C::value_type pick(const C &container) {991 return pick(container.begin(), container.end());992}993template <typename T> T pick(const std::initializer_list<T> &il) {994 return pick(std::vector<T>(il));995}996997// Returns container[i] with probability proportional to distribution[i].998// O(1) for random_access_iterator, O(|container|) otherwise.999template <typename C, typename T>1000typename C::value_type pick_by_distribution(const C &container,1001 std::vector<T> distribution) {1002 tgen_ensure(container.size() == distribution.size(),1003 "container and distribution must have the same size");1004 auto it = container.begin();1005 std::advance(it, next_by_distribution(distribution));1006 return *it;1007}1008template <typename C, typename T>1009typename C::value_type1010pick_by_distribution(const C &container,1011 const std::initializer_list<T> &distribution) {1012 return pick_by_distribution(container, std::vector<T>(distribution));1013}1014template <typename T, typename U>1015T pick_by_distribution(const std::initializer_list<T> &il,1016 const std::vector<U> &distribution) {1017 return pick_by_distribution(std::vector<T>(il), distribution);1018}1019template <typename T, typename U>1020T pick_by_distribution(const std::initializer_list<T> &il,1021 const std::initializer_list<U> &distribution) {1022 return pick_by_distribution(std::vector<T>(il),1023 std::vector<U>(distribution));1024}10251026// Chooses k values uniformly from container, as in a subsequence of size k.1027// Returns a copy. O(|container|).1028template <typename C> C choose(const C &container, int k) {1029 tgen_ensure(0 < k and k <= static_cast<int>(container.size()),1030 "number of elements to choose must be valid");1031 std::vector<typename C::value_type> new_vec;1032 C new_container;1033 int need = k, left = container.size();1034 for (auto cur_it = container.begin(); cur_it != container.end(); ++cur_it) {1035 if (next(1, left--) <= need) {1036 new_container.insert(new_container.end(), *cur_it);1037 need--;1038 }1039 }1040 return new_container;1041}1042template <typename T>1043std::vector<T> choose(const std::initializer_list<T> &il, int k) {1044 return choose(std::vector<T>(il), k);1045}10461047// Number distinct generator for integral types.1048// Optimized for performance (unordered_map virtual list; gen_list uses array1049// pool, complement, or sparse sampling).1050template <typename T> struct distinct_range {1051 T left_, right_;1052 T num_available_;1053 std::unordered_map<T, T> virtual_list_;10541055 // When the range fits in memory, sample via array Fisher–Yates.1056 static constexpr size_t array_pool_max = size_t{1} << 23;10571058 // Generator of distinct values in [left, right].1059 distinct_range(T left, T right)1060 : left_(left), right_(right), num_available_(right - left + 1) {}10611062 // Returns the number of distinct values left to generate.1063 T size() const { return num_available_; }10641065 // Generates a random value in [left_, right_] that has not been generated1066 // yet.1067 // O(log n).1068 T gen() {1069 // One iteration of Fisher–Yates.1070 tgen_ensure(size() > 0, "distinct_range: no more values to generate");10711072 T i = next<T>(0, size() - 1);1073 T j = size() - 1;10741075 auto vi_it = virtual_list_.find(i);1076 T vi = vi_it == virtual_list_.end() ? i : vi_it->second;1077 auto vj_it = virtual_list_.find(j);1078 T vj = vj_it == virtual_list_.end() ? j : vj_it->second;1079 virtual_list_[i] = vj;10801081 --num_available_;10821083 return vi + left_;1084 }10851086 // Generates a list of distinct values.1087 // Optimized for performance (array pool, complement, or sparse sampling).1088 // O(size) when the range fits in memory; O(size log range) otherwise.1089 auto gen_list(int count) {1090 tgen_ensure(count >= 0, "distinct_range: size must be nonnegative");1091 tgen_ensure(count <= num_available_,1092 "distinct_range: no more values to generate");10931094 size_t range_size = right_ - left_ + 1;1095 size_t sample_count = count;10961097 std::vector<T> res;1098 if (sample_count > 0) {1099 if (range_size <= array_pool_max)1100 res = sample_from_pool(sample_count, range_size);1101 else if (sample_count * 2 > range_size)1102 res = sample_complement(sample_count, range_size);1103 else1104 res = sample_sparse(sample_count);1105 }11061107 num_available_ -= count;1108 virtual_list_.clear();1109 return typename list<T>::value(res);1110 }11111112 // Generates all distinct values.1113 // O(n) when the range fits in memory; O(n log n) otherwise.1114 auto gen_all() { return gen_list(size()); }11151116 private:1117 // Samples count distinct values via array Fisher–Yates on [left_, right_].1118 // O(range_size) time and memory.1119 std::vector<T> sample_from_pool(size_t count, size_t range_size) {1120 std::vector<T> pool(range_size);1121 std::iota(pool.begin(), pool.end(), left_);1122 for (size_t i = 0; i < count; ++i) {1123 size_t j = next<size_t>(i, range_size - 1);1124 std::swap(pool[i], pool[j]);1125 }1126 pool.resize(count);1127 return pool;1128 }11291130 // Samples count distinct values by excluding range_size - count values.1131 // O(range_size + (range_size - count) log(range_size)).1132 std::vector<T> sample_complement(size_t count, size_t range_size) {1133 size_t exclude_count = range_size - count;1134 std::unordered_set<T> excluded;1135 excluded.reserve(exclude_count * 2);11361137 if (exclude_count <= array_pool_max) {1138 for (T value : sample_from_pool(exclude_count, range_size))1139 excluded.insert(value);1140 } else {1141 for (T value : sample_sparse(exclude_count))1142 excluded.insert(value);1143 }11441145 std::vector<T> res;1146 res.reserve(count);1147 for (T value = left_; value <= right_; ++value) {1148 if (!excluded.count(value))1149 res.push_back(value);1150 }1151 detail::tgen_ensure_against_bug(1152 res.size() == count, "distinct_range: complement sampling failed");1153 return res;1154 }11551156 // Samples count distinct values via sparse-map Fisher–Yates.1157 // O(count log(range_size)).1158 std::vector<T> sample_sparse(size_t count) {1159 std::unordered_map<T, T> local_virtual;1160 local_virtual.reserve(count * 2);1161 T remaining = range_span();1162 std::vector<T> res;1163 res.reserve(count);1164 for (size_t step = 0; step < count; ++step) {1165 T i = next<T>(0, remaining - 1);1166 T j = remaining - 1;11671168 auto vi_it = local_virtual.find(i);1169 T vi = vi_it == local_virtual.end() ? i : vi_it->second;1170 auto vj_it = local_virtual.find(j);1171 T vj = vj_it == local_virtual.end() ? j : vj_it->second;1172 local_virtual[i] = vj;11731174 res.push_back(vi + left_);1175 --remaining;1176 }1177 return res;1178 }11791180 // Returns right_ - left_ + 1.1181 // O(1).1182 T range_span() { return right_ - left_ + 1; }1183};11841185// Distinct generator for containers.1186template <typename T> struct distinct_container {1187 std::vector<T> list_;1188 distinct_range<size_t> idx_;11891190 // Creates a distinct container generator for the given container.1191 template <typename C>1192 distinct_container(const C &container)1193 : list_(container.begin(), container.end()),1194 idx_(0, static_cast<int>(container.size()) - 1) {}1195 distinct_container(const std::initializer_list<T> &il)1196 : distinct_container(std::vector<T>(il)) {}11971198 // Returns the number of distinct elements left to generate.1199 size_t size() const { return idx_.size(); }12001201 // Generates a random element from container uniformly.1202 // O(log n).1203 T gen() { return list_[idx_.gen()]; }12041205 // Generates a list of distinct values.1206 // O(size * log(n)).1207 auto gen_list(int size) {1208 std::vector<T> res;1209 for (int i = 0; i < size; ++i)1210 res.push_back(gen());1211 return typename list<T>::value(res);1212 }12131214 // Generates all distinct values.1215 // O(n log(n))1216 auto gen_all() {1217 std::vector<T> res;1218 while (size() > 0)1219 res.push_back(gen());1220 return typename list<T>::value(res);1221 }1222};1223template <typename C>1224distinct_container(const C &) -> distinct_container<typename C::value_type>;12251226/************1227 * *1228 * OPTS *1229 * *1230 ************/12311232/*1233 * Opts - options given to the generator.1234 *1235 * Incompatible with testlib.1236 *1237 * Opts are a list of either positional or named options.1238 *1239 * Named options are given in one of the following formats:1240 * 1) -keyname=value or --keyname=value (ex. -n=10 , --test-count=20)1241 * 2) -keyname value or --keyname value (ex. -n 10 , --test-count 20)1242 *1243 * Positional options are numbered from 0 sequentially.1244 * For example, for "10 -n=20 str" positional option 1 is the string "str".1245 */12461247/*1248 * C++ version selection.1249 */12501251// Sets C++ version.1252inline void set_cpp_version(int version) {1253 detail::cpp = detail::cpp_value(version);1254}12551256/*1257 * Compiler selection.1258 */12591260// GCC compiler type.1261inline compiler_value gcc(int major = 0, int minor = 0) {1262 return {compiler_kind::gcc, major, minor};1263}12641265// Clang compiler type.1266inline compiler_value clang(int major = 0, int minor = 0) {1267 return {compiler_kind::clang, major, minor};1268}12691270// Sets compiler.1271inline void set_compiler(compiler_value compiler) {1272 detail::compiler.kind_ = compiler.kind_;1273 detail::compiler.major_ = compiler.major_;1274 detail::compiler.minor_ = compiler.minor_;1275}12761277namespace detail {12781279// Processes special opt flags.1280// Returns true if the key is a special opt flag.1281inline bool process_special_opt_flags(std::string &key) {1282 // Checks for gen::CPP=17|20|231283 if (key.find("tgen::CPP:") == 0) {1284 int prefix_len = std::string("tgen::CPP:").size();1285 tgen_ensure(static_cast<int>(key.size()) == prefix_len + 2 and1286 std::isdigit(key[prefix_len]) and1287 std::isdigit(key[prefix_len + 1]),1288 "invalid CPP format");1289 int version = std::stoi(key.substr(prefix_len, 2));1290 set_cpp_version(version);1291 return true;1292 }12931294 // Checks for tgen::(GCC|CLANG) or1295 // tgen::(GCC|CLANG):(version|version.minor).1296 compiler_kind kind;1297 size_t prefix_len = 0;12981299 if (key.find("tgen::GCC") == 0) {1300 kind = compiler_kind::gcc;1301 prefix_len = std::string("tgen::GCC").size();1302 } else if (key.find("tgen::CLANG") == 0) {1303 kind = compiler_kind::clang;1304 prefix_len = std::string("tgen::CLANG").size();1305 } else {1306 return false;1307 }13081309 if (key.size() == prefix_len) {1310 set_compiler(compiler_value(kind, 0, 0));1311 return true;1312 }13131314 tgen_ensure(key[prefix_len] == ':', "invalid compiler format");1315 ++prefix_len; // for ':'.13161317 std::string inside = key.substr(prefix_len, key.size() - prefix_len);1318 int major = 0, minor = 0;13191320 size_t dot = inside.find('.');1321 if (dot == std::string::npos) {1322 tgen_ensure(!inside.empty() and1323 std::all_of(inside.begin(), inside.end(), ::isdigit),1324 "invalid compiler version");1325 major = std::stoi(inside);1326 } else {1327 std::string maj = inside.substr(0, dot);1328 std::string min = inside.substr(dot + 1);13291330 tgen_ensure(!maj.empty() and1331 std::all_of(maj.begin(), maj.end(), ::isdigit) and1332 maj.size() <= 3,1333 "invalid compiler major version");1334 tgen_ensure(!min.empty() and1335 std::all_of(min.begin(), min.end(), ::isdigit) and1336 min.size() <= 3,1337 "invalid compiler minor version");13381339 major = std::stoi(maj);1340 minor = std::stoi(min);1341 }13421343 set_compiler(compiler_value(kind, major, minor));13441345 return true;1346}13471348inline std::vector<std::string>1349 pos_opts; // Dictionary containing the positional parsed opts.1350inline std::map<std::string, std::string>1351 named_opts; // Global dictionary the named parsed opts.13521353template <typename T> T get_opt(const std::string &value) {1354 try {1355 if constexpr (std::is_same_v<T, bool>) {1356 if (value == "true" or value == "1")1357 return true;1358 if (value == "false" or value == "0")1359 return false;1360 } else if constexpr (std::is_integral_v<T>) {1361 if constexpr (std::is_unsigned_v<T>)1362 return static_cast<T>(std::stoull(value));1363 else1364 return static_cast<T>(std::stoll(value));1365 } else if constexpr (std::is_floating_point_v<T>)1366 return static_cast<T>(std::stold(value));1367 else1368 return value; // Default: std::string.1369 } catch (...) {1370 }13711372 throw error("invalid value `" + value + "` for type " + typeid(T).name());1373}13741375inline void parse_opts(int argc, char **argv) {1376 // Parses the opts into `pos_opts` vector and `named_opts`1377 // map. Starting from 1 to ignore the name of the executable.1378 for (int i = 1; i < argc; ++i) {1379 std::string key(argv[i]);13801381 if (process_special_opt_flags(key))1382 continue;13831384 if (key[0] == '-') {1385 tgen_ensure(key.size() > 1,1386 "invalid opt (" + std::string(argv[i]) + ")");1387 if ('0' <= key[1] and key[1] <= '9') {1388 // This case is a positional negative number argument.1389 pos_opts.push_back(key);1390 continue;1391 }13921393 // Pops first char '-'.1394 key = key.substr(1);1395 } else {1396 // This case is a positional argument that does not start with '-'.1397 pos_opts.push_back(key);1398 continue;1399 }14001401 // Pops a possible second char '-'.1402 if (key[0] == '-') {1403 tgen_ensure(key.size() > 1,1404 "invalid opt (" + std::string(argv[i]) + ")");14051406 // Pops first char '-'.1407 key = key.substr(1);1408 }14091410 // Assumes that, if it starts with '-' and second char is not a digit,1411 // then it is a <key, value> pair.1412 // 1 or 2 chars '-' have already been popped.14131414 std::size_t eq = key.find('=');1415 if (eq != std::string::npos) {1416 // This is the '--key=value' case.1417 std::string value = key.substr(eq + 1);1418 key = key.substr(0, eq);1419 tgen_ensure(!key.empty() and !value.empty(),1420 "expected non-empty key/value in opt (" +1421 std::string(argv[i]) + ")");1422 tgen_ensure(named_opts.count(key) == 0,1423 "cannot have repeated keys");1424 named_opts[key] = value;1425 } else {1426 // This is the '--key value' case.1427 tgen_ensure(named_opts.count(key) == 0,1428 "cannot have repeated keys");1429 tgen_ensure(argv[i + 1], "value cannot be empty");1430 named_opts[key] = std::string(argv[i + 1]);1431 ++i;1432 }1433 }1434}14351436inline void set_seed(int argc, char **argv) {1437 std::vector<uint32_t> seed;14381439 // Starting from 1 to ignore the name of the executable.1440 for (int i = 1; i < argc; ++i) {1441 // We append the number of chars, and then the list of chars.1442 int size_pos = seed.size();1443 seed.push_back(0);1444 for (char *s = argv[i]; *s != '\0'; ++s) {1445 ++seed[size_pos];1446 seed.push_back(*s);1447 }1448 }1449 std::seed_seq seq(seed.begin(), seed.end());1450 rng.seed(seq);1451}14521453} // namespace detail14541455// Returns true if there is an opt at a given index.1456inline bool has_opt(std::size_t index) {1457 detail::ensure_registered();1458 return index < detail::pos_opts.size();1459}14601461// Returns true if there is an opt with a given key.1462inline bool has_opt(const std::string &key) {1463 detail::ensure_registered();1464 return detail::named_opts.count(key) != 0;1465}1466template <typename K>1467std::enable_if_t<std::is_same_v<K, char>, bool> has_opt(K key) {1468 return has_opt(std::string(1, key));1469}14701471// Returns the parsed opt by a given index. If no opts with the given index are1472// found, returns the given default_value.1473template <typename T>1474T opt(size_t index, std::optional<T> default_value = std::nullopt) {1475 detail::ensure_registered();1476 if (!has_opt(index)) {1477 if (default_value)1478 return *default_value;1479 throw detail::error("cannot find opt at index " +1480 std::to_string(index));1481 }1482 return detail::get_opt<T>(detail::pos_opts[index]);1483}14841485// Returns the parsed opt by a given key. If no opts with the given key are1486// found, returns the given default_value.1487template <typename T>1488T opt(const std::string &key, std::optional<T> default_value = std::nullopt) {1489 detail::ensure_registered();1490 if (!has_opt(key)) {1491 if (default_value)1492 return *default_value;1493 throw detail::error("cannot find opt with key " + key);1494 }1495 return detail::get_opt<T>(detail::named_opts[key]);1496}1497template <typename T, typename K>1498std::enable_if_t<std::is_same_v<K, char>, T>1499opt(K key, std::optional<T> default_value = std::nullopt) {1500 return opt<T>(std::string(1, key), default_value);1501}15021503// Registers generator by initializing rng and parsing opts.1504inline void register_gen(int argc, char **argv) {1505 detail::set_seed(argc, argv);15061507 detail::pos_opts.clear();1508 detail::named_opts.clear();1509 detail::parse_opts(argc, argv);15101511 detail::registered = true;1512}15131514// Registers generator by initializing rng with a given seed.1515inline void register_gen(std::optional<long long> seed = std::nullopt) {1516 if (seed)1517 detail::rng.seed(*seed);1518 else1519 detail::rng.seed();15201521 detail::pos_opts.clear();1522 detail::named_opts.clear();15231524 detail::registered = true;1525}15261527/************1528 * *1529 * LIST *1530 * *1531 ************/15321533/*1534 * List generator.1535 *1536 * List of integral types.1537 */15381539template <typename T> struct list : gen_base<list<T>> {1540 int size_; // Size of list.1541 T value_l_, value_r_; // Range of defined values.1542 std::set<T> values_; // Set of values. If empty, use range; if not,1543 // represents the possible values, and the range1544 // represents the index in this set.1545 std::map<T, int>1546 value_idx_in_set_; // Index of every value in the set above.1547 mutable std::vector<std::pair<T, T>>1548 val_range_; // Range of values of each index.1549 mutable std::vector<std::vector<int>> neigh_; // Adjacency list of equality.1550 std::vector<std::set<int>>1551 diff_restrictions_; // All different restrictions.1552 bool index_constraints_{1553 false}; // True after fix/equal narrows per-index generation.1554 mutable bool uses_full_range_{1555 false}; // If true, every index uses [value_l_, value_r_] lazily.15561557 // Creates generator for lists of size 'size', with random T in [value_left,1558 // value_right].1559 list(int size, T value_left, T value_right)1560 : size_(size), value_l_(value_left), value_r_(value_right),1561 uses_full_range_(true) {1562 tgen_ensure(size_ > 0, "list: size must be positive");1563 tgen_ensure(value_l_ <= value_r_, "list: value range must be valid");1564 }15651566 // Creates list with value set.1567 list(int size, std::set<T> values)1568 : size_(size), values_(values), index_constraints_(true) {1569 tgen_ensure(size_ > 0, "list: size must be positive");1570 tgen_ensure(!values.empty(), "list: value set must be non-empty");1571 value_l_ = 0, value_r_ = values.size() - 1;1572 val_range_.assign(size_, {value_l_, value_r_});1573 int idx = 0;1574 for (T val : values_)1575 value_idx_in_set_[val] = idx++;1576 }15771578 // Restricts lists for list[idx] = val.1579 list &fix(int idx, T val) {1580 tgen_ensure(0 <= idx and idx < size_, "list: index must be valid");1581 ensure_val_range_materialized();1582 if (values_.size() == 0) {1583 auto &[left, right] = val_range_[idx];1584 if (left == right and value_l_ != value_r_) {1585 tgen_ensure(left == val,1586 "list: must not set to two different values");1587 } else {1588 tgen_ensure(left <= val and val <= right,1589 "list: value must be in the defined range");1590 }1591 left = right = val;1592 } else {1593 tgen_ensure(values_.count(val),1594 "list: value must be in the set of values");1595 auto &[left, right] = val_range_[idx];1596 int new_val = value_idx_in_set_[val];1597 tgen_ensure(left <= new_val and new_val <= right,1598 "list: must not set to two different values");1599 left = right = new_val;1600 }1601 index_constraints_ = true;1602 return *this;1603 }16041605 // Restricts lists for list[idx_1] = list[idx_2].1606 list &equal(int idx_1, int idx_2) {1607 tgen_ensure(0 <= std::min(idx_1, idx_2) and1608 std::max(idx_1, idx_2) < size_,1609 "list: indices must be valid");1610 if (idx_1 == idx_2)1611 return *this;16121613 ensure_val_range_materialized();1614 ensure_neigh_allocated();1615 index_constraints_ = true;1616 neigh_[idx_1].push_back(idx_2);1617 neigh_[idx_2].push_back(idx_1);1618 return *this;1619 }16201621 // Restricts lists for list[S] to be equal, for given subset S of indices.1622 list &equal(std::set<int> indices) {1623 if (!indices.empty()) {1624 std::set<int>::iterator beg = indices.begin();1625 for (auto it = std::next(beg); it != indices.end(); ++it)1626 equal(*beg, *it);1627 }1628 return *this;1629 }16301631 // Restricts lists for list[left..right] to have all equal values.1632 list &equal_range(int left, int right) {1633 tgen_ensure(0 <= left and left <= right and right < size_,1634 "list: range indices must be valid");1635 for (int i = left; i < right; ++i)1636 equal(i, i + 1);1637 return *this;1638 }16391640 // Restricts lists for all equal elements.1641 list &all_equal() { return equal_range(0, size_ - 1); }16421643 // Restricts lists for list[S] to be different (distinct), for given subset1644 // S of indices. You cannot add two of these restrictions on sets that1645 // intersect.1646 list &different(std::set<int> indices) {1647 if (!indices.empty())1648 diff_restrictions_.push_back(indices);1649 return *this;1650 }16511652 // Restricts lists for list[idx_1] != list[idx_2].1653 list &different(int idx_1, int idx_2) {1654 std::set<int> indices = {idx_1, idx_2};1655 return different(indices);1656 }16571658 // Restricts lists for list[left..right] to have all different values.1659 list &different_range(int left, int right) {1660 tgen_ensure(0 <= left and left <= right and right < size_,1661 "list: range indices must be valid");1662 std::vector<int> indices(right - left + 1);1663 std::iota(indices.begin(), indices.end(), left);1664 return different(std::set<int>(indices.begin(), indices.end()));1665 }16661667 // Restricts lists for all different elements.1668 list &all_different() {1669 std::vector<int> indices(size_);1670 std::iota(indices.begin(), indices.end(), 0);1671 return different(std::set<int>(indices.begin(), indices.end()));1672 }16731674 // List value.1675 // Operations on a value are not random.1676 struct value : gen_value_base<value> {1677 using tgen_is_sequential_tag = detail::is_sequential_tag;16781679 using value_type = T; // Value type, for templates.1680 using std_type = std::vector<T>; // std type for value.16811682 std::vector<T> vec_; // list.1683 char sep_; // Separator for printing.16841685 value(const std::vector<T> &vec) : vec_(vec), sep_(' ') {}1686 value(const std::initializer_list<T> &il) : value(std::vector<T>(il)) {}16871688 // Fetches size.1689 int size() const { return vec_.size(); }16901691 // Fetches position idx.1692 T &operator[](int idx) {1693 tgen_ensure(0 <= idx and idx < size(),1694 "list: value: index out of bounds");1695 return vec_[idx];1696 }1697 const T &operator[](int idx) const {1698 tgen_ensure(0 <= idx and idx < size(),1699 "list: value: index out of bounds");1700 return vec_[idx];1701 }17021703 // Sorts values in non-decreasing order.1704 // O(n log n).1705 value &sort() {1706 std::sort(vec_.begin(), vec_.end());1707 return *this;1708 }17091710 // Reverses list.1711 // O(n).1712 value &reverse() {1713 std::reverse(vec_.begin(), vec_.end());1714 return *this;1715 }17161717 // Sets the separator for the list, for printing.1718 // O(1).1719 value &separator(char sep) {1720 sep_ = sep;1721 return *this;1722 }17231724 // Concatenates two values.1725 // Linear.1726 value operator+(const value &rhs) const {1727 std::vector<T> new_vec = vec_;1728 for (int i = 0; i < rhs.size(); ++i)1729 new_vec.push_back(rhs[i]);1730 return value(new_vec);1731 }17321733 // Shuffles list uniformly.1734 // O(n).1735 value &shuffle() {1736 for (int i = 0; i < size(); ++i)1737 std::swap(vec_[i], vec_[next(0, size() - 1)]);1738 return *this;1739 }17401741 // Returns a random element uniformly.1742 // O(1).1743 T pick() const { return vec_[next<int>(0, size() - 1)]; }17441745 // Returns vec_[i] with probability proportional to distribution[i].1746 // O(1).1747 template <typename Dist>1748 T pick_by_distribution(const std::vector<Dist> &distribution) const {1749 tgen_ensure(static_cast<size_t>(size()) == distribution.size(),1750 "value and distribution must have the same size");1751 return vec_[next_by_distribution(distribution)];1752 }1753 template <typename Dist>1754 T pick_by_distribution(1755 const std::initializer_list<Dist> &distribution) const {1756 return pick_by_distribution(std::vector<Dist>(distribution));1757 }17581759 // Chooses k values uniformly, as in a subsequence of size k.1760 // O(n).1761 value choose(int k) const {1762 tgen_ensure(0 < k and k <= size(),1763 "number of elements to choose must be valid");1764 std::vector<T> new_vec;1765 int need = k;1766 for (int i = 0; need > 0; ++i) {1767 int left = size() - i;1768 if (next(1, left) <= need) {1769 new_vec.push_back(vec_[i]);1770 need--;1771 }1772 }1773 return value(new_vec);1774 }17751776 // Prints to std::ostream, separated by sep_.1777 friend std::ostream &operator<<(std::ostream &out, const value &val) {1778 for (int i = 0; i < val.size(); ++i) {1779 if (i > 0)1780 out << val.sep_;1781 out << val[i];1782 }1783 return out;1784 }17851786 // Gets a std::vector representing the value.1787 auto to_std() const {1788 if constexpr (!detail::is_generator_value<T>::value) {1789 return vec_;1790 } else {1791 std::vector<typename T::std_type> vec;1792 for (const auto &i : vec_)1793 vec.push_back(i.to_std());1794 return vec;1795 }1796 }1797 };17981799 // Generates list value.1800 // Optimized for performance (unconstrained and all-different fast paths).1801 // O(n log n).1802 value gen() const {1803 if (diff_restrictions_.empty()) {1804 if (auto unconstrained = try_gen_unconstrained())1805 return *unconstrained;1806 }1807 if (auto all_different = try_gen_all_different())1808 return *all_different;18091810 ensure_neigh_allocated();1811 std::vector<T> vec(size_);1812 std::vector<bool> defined_idx(1813 size_, false); // For every index, if it has been set in `vec`.18141815 std::vector<int> comp_id(size_, -1); // Component id of each index.1816 std::vector<std::vector<int>> comp(size_); // Component of each comp-id.1817 int comp_count = 0; // Number of different components.18181819 // Defines value of entire component.1820 auto define_comp = [&](int cur_comp, T val) {1821 for (int idx : comp[cur_comp]) {1822 tgen_ensure(!defined_idx[idx]);1823 vec[idx] = val;1824 defined_idx[idx] = true;1825 }1826 };18271828 // Groups = components.1829 {1830 std::vector<bool> vis(size_, false); // Visited for each index.1831 for (int idx = 0; idx < size_; ++idx)1832 if (!vis[idx]) {1833 T new_value;1834 bool value_defined = false;18351836 // BFS to visit the connected component, grouping equal1837 // values.1838 std::queue<int> q({idx});1839 vis[idx] = true;1840 std::vector<int> component;1841 while (!q.empty()) {1842 int cur_idx = q.front();1843 q.pop();18441845 component.push_back(cur_idx);18461847 // Checks value.1848 auto [l, r] = val_range_at(cur_idx);1849 if (l == r) {1850 if (!value_defined) {1851 // We found the value.1852 value_defined = true;1853 new_value = l;1854 } else if (new_value != l) {1855 // We found a contradiction1856 throw detail::contradiction_error(1857 "list",1858 "tried to set value to `" +1859 std::to_string(new_value) +1860 "`, but it was already set as `" +1861 std::to_string(l) + "`");1862 }1863 }18641865 for (int nxt_idx : neigh_[cur_idx]) {1866 if (!vis[nxt_idx]) {1867 vis[nxt_idx] = true;1868 q.push(nxt_idx);1869 }1870 }1871 }18721873 // Group entire component, checking if value is defined.1874 for (int cur_idx : component) {1875 comp_id[cur_idx] = comp_count;1876 comp[comp_id[cur_idx]].push_back(cur_idx);1877 }18781879 // Defines value if needed.1880 if (value_defined)1881 define_comp(comp_count, new_value);18821883 ++comp_count;1884 }1885 }18861887 // Initial parsing of different restrictions.1888 std::vector<std::set<int>> diff_containing_comp_idx(comp_count);1889 {1890 int dist_id = 0;1891 for (const std::set<int> &diff : diff_restrictions_) {1892 // Checks if there are too many different values.1893 if (static_cast<uint64_t>(diff.size() - 1) +1894 static_cast<uint64_t>(value_l_) >1895 static_cast<uint64_t>(value_r_))1896 throw detail::contradiction_error(1897 "list", "tried to generate " +1898 std::to_string(diff.size()) +1899 " different values, but the maximum is " +1900 std::to_string(value_r_ - value_l_ + 1));19011902 // Checks if two values in same component are marked as1903 // different.1904 std::set<int> comp_ids;1905 for (int idx : diff) {1906 if (comp_ids.count(comp_id[idx]))1907 throw detail::contradiction_error(1908 "list", "tried to set two indices as equal and "1909 "different");1910 comp_ids.insert(comp_id[idx]);19111912 diff_containing_comp_idx[comp_id[idx]].insert(dist_id);1913 }1914 ++dist_id;1915 }1916 }19171918 // If some value is in >= 3 sets, then there is a cycle.1919 for (auto &diff_containing : diff_containing_comp_idx)1920 if (diff_containing.size() >= 3)1921 throw detail::complex_restrictions_error(1922 "list",1923 "one index cannot be in >= 3 'different' restrictions");19241925 std::vector<bool> vis_diff(diff_restrictions_.size(), false);1926 std::vector<bool> initially_defined_comp_idx(comp_count, false);19271928 // Fills the value in a tree defined by "different" restrictions.1929 auto define_tree = [&](int diff_id) {1930 // The set `diff_restrictions_[diff_id]` can have some1931 // values that are defined.19321933 // Generates set of already defined values.1934 std::set<T> defined_values;1935 for (int idx : diff_restrictions_[diff_id])1936 if (defined_idx[idx]) {1937 // Checks if two values in `diff_restrictions_[dist_id]`1938 // have been set to the same value1939 if (defined_values.count(vec[idx]))1940 throw detail::contradiction_error(1941 "list",1942 "tried to set two indices as equal and different");19431944 defined_values.insert(vec[idx]);1945 }19461947 // Generates values in this root "different" restriction.1948 {1949 int new_value_count = diff_restrictions_[diff_id].size() -1950 static_cast<int>(defined_values.size());1951 std::vector<T> generated_values =1952 generate_distinct_values(new_value_count, defined_values);1953 auto val_it = generated_values.begin();1954 for (int idx : diff_restrictions_[diff_id])1955 if (defined_idx[idx]) {1956 // The root can cover these components, but there should1957 // not be any other defined in this tree.1958 initially_defined_comp_idx[comp_id[idx]] = false;1959 } else {1960 define_comp(comp_id[idx], *val_it);1961 ++val_it;1962 }1963 }19641965 // BFS on the tree of "different" restrictions.1966 std::queue<std::pair<int, int>> q; // {id, parent id}1967 q.emplace(diff_id, -1);1968 vis_diff[diff_id] = true;1969 while (!q.empty()) {1970 auto [cur_diff, parent] = q.front();1971 q.pop();19721973 std::set<int> neigh_diff;1974 for (int idx : diff_restrictions_[cur_diff])1975 for (int nxt_diff :1976 diff_containing_comp_idx[comp_id[idx]]) {1977 if (nxt_diff == cur_diff or nxt_diff == parent)1978 continue;19791980 // Cycle found.1981 if (vis_diff[nxt_diff])1982 throw detail::complex_restrictions_error(1983 "list",1984 "cycle found in 'different' restrictions");19851986 neigh_diff.insert(nxt_diff);1987 }19881989 for (int nxt_diff : neigh_diff) {1990 vis_diff[nxt_diff] = true;1991 q.emplace(nxt_diff, cur_diff);19921993 // Generates this "different" restriction.1994 std::set<T> nxt_defined_values;1995 for (int idx2 : diff_restrictions_[nxt_diff])1996 if (defined_idx[idx2]) {1997 // There cannot be any more defined. This case is1998 // when there are values not covered by a single1999 // "different" restriction in the tree.2000 if (initially_defined_comp_idx[comp_id[idx2]])2001 throw detail::complex_restrictions_error(2002 "list");20032004 nxt_defined_values.insert(vec[idx2]);2005 }2006 int new_value_count =2007 diff_restrictions_[nxt_diff].size() -2008 static_cast<int>(nxt_defined_values.size());2009 std::vector<T> generated_values = generate_distinct_values(2010 new_value_count, nxt_defined_values);2011 auto val_it = generated_values.begin();2012 for (int idx2 : diff_restrictions_[nxt_diff])2013 if (!defined_idx[idx2]) {2014 define_comp(comp_id[idx2], *val_it);2015 ++val_it;2016 }2017 }2018 }2019 };20202021 // Loops through "different" restrictions, sorts "different"2022 // restrictions by number of defined components (non-increasing). This2023 // guarantees that if there is a valid root (that covers all 'defined'),2024 // we find it.2025 {2026 std::vector<std::pair<int, int>> defined_cnt_and_diff_idx;2027 int dist_id = 0;2028 for (const std::set<int> &diff : diff_restrictions_) {2029 int defined_cnt = 0;2030 for (int idx : diff)2031 if (defined_idx[idx]) {2032 ++defined_cnt;2033 initially_defined_comp_idx[comp_id[idx]] = true;2034 }2035 defined_cnt_and_diff_idx.emplace_back(defined_cnt, dist_id);2036 ++dist_id;2037 }20382039 std::sort(defined_cnt_and_diff_idx.rbegin(),2040 defined_cnt_and_diff_idx.rend());2041 for (auto [defined_cnt, diff_idx] : defined_cnt_and_diff_idx)2042 if (!vis_diff[diff_idx])2043 define_tree(diff_idx);2044 }20452046 // Loops through "different" restrictions do define the rest.2047 for (std::size_t dist_id = 0; dist_id < diff_restrictions_.size();2048 ++dist_id)2049 if (!vis_diff[dist_id])2050 define_tree(dist_id);20512052 // Define final values. These values all should be random in [l, r], and2053 // the "different" restrictions have already been processed. However,2054 // there can be still equality restrictions, so we define entire2055 // components.2056 for (int idx = 0; idx < size_; ++idx)2057 if (!defined_idx[idx])2058 define_comp(comp_id[idx], next<T>(value_l_, value_r_));20592060 if (!values_.empty()) {2061 // Needs to fetch the values from the value set.2062 std::vector<T> value_vec(values_.begin(), values_.end());2063 for (T &val : vec)2064 val = value_vec[val];2065 }20662067 return value(vec);2068 }20692070 private:2071 // Materializes neigh_ after the first equality restriction.2072 void ensure_neigh_allocated() const {2073 if (neigh_.size() == static_cast<size_t>(size_))2074 return;2075 neigh_.assign(size_, {});2076 }20772078 // Materializes val_range_ after the first per-index restriction.2079 void ensure_val_range_materialized() const {2080 if (!uses_full_range_)2081 return;2082 val_range_.assign(size_, {value_l_, value_r_});2083 uses_full_range_ = false;2084 }20852086 // Returns the allowed value range at index idx.2087 std::pair<T, T> val_range_at(int idx) const {2088 if (uses_full_range_)2089 return {value_l_, value_r_};2090 return val_range_[idx];2091 }20922093 // Generates a uniformly random list of k distinct values in `[value_l,2094 // value_r]`, such that no value is in `forbidden_values`.2095 std::vector<T>2096 generate_distinct_values(int k, const std::set<T> &forbidden_values) const {2097 for (auto forbidden : forbidden_values)2098 tgen_ensure(value_l_ <= forbidden and forbidden <= value_r_);2099 const T num_available =2100 (value_r_ - value_l_ + 1) - forbidden_values.size();2101 if (num_available < k)2102 throw detail::complex_restrictions_error(2103 "list", "not enough distinct values");2104 if (forbidden_values.empty())2105 return distinct_range<T>(value_l_, value_r_).gen_list(k).to_std();21062107 std::map<T, T> virtual_list;2108 std::vector<T> gen_list;2109 for (int i = 0; i < k; ++i) {2110 T j = next<T>(i, num_available - 1);2111 T vj = virtual_list.count(j) ? virtual_list[j] : j;2112 T vi = virtual_list.count(i) ? virtual_list[i] : i;21132114 virtual_list[j] = vi, virtual_list[i] = vj;21152116 gen_list.push_back(virtual_list[i]);2117 }21182119 for (T &val : gen_list)2120 val += value_l_;21212122 std::vector<std::pair<T, int>> values_sorted;2123 for (std::size_t i = 0; i < gen_list.size(); ++i)2124 values_sorted.emplace_back(gen_list[i], i);2125 std::sort(values_sorted.begin(), values_sorted.end());2126 auto cur_it = forbidden_values.begin();2127 int smaller_forbidden_count = 0;2128 for (auto [val, idx] : values_sorted) {2129 while (cur_it != forbidden_values.end() and2130 *cur_it <= val + smaller_forbidden_count)2131 ++cur_it, ++smaller_forbidden_count;2132 gen_list[idx] += smaller_forbidden_count;2133 }21342135 return gen_list;2136 }21372138 // If this generator has no constraints beyond [value_l_, value_r_],2139 // returns independent uniform samples; otherwise returns std::nullopt.2140 // O(n).2141 std::optional<value> try_gen_unconstrained() const {2142 if (!values_.empty() or index_constraints_)2143 return std::nullopt;21442145 std::vector<T> vec(size_);2146 for (int i = 0; i < size_; ++i)2147 vec[i] = next<T>(value_l_, value_r_);2148 return value(vec);2149 }21502151 // If this generator is exactly all-distinct in [value_l_, value_r_],2152 // returns a uniformly random list; otherwise returns std::nullopt.2153 // Optimized for performance (distinct_range fast path).2154 // O(n log n).2155 std::optional<value> try_gen_all_different() const {2156 if (!values_.empty() or diff_restrictions_.size() != 1)2157 return std::nullopt;21582159 const std::set<int> &diff = diff_restrictions_[0];2160 if (static_cast<int>(diff.size()) != size_ or *diff.begin() != 0 or2161 *diff.rbegin() != size_ - 1)2162 return std::nullopt;21632164 if (!neigh_.empty()) {2165 for (const auto &adj : neigh_) {2166 if (!adj.empty())2167 return std::nullopt;2168 }2169 }21702171 if (index_constraints_)2172 return std::nullopt;21732174 if (static_cast<long long>(size_) >2175 static_cast<long long>(value_r_) - value_l_ + 1)2176 throw detail::contradiction_error(2177 "list", "tried to generate " + std::to_string(size_) +2178 " different values, but the maximum is " +2179 std::to_string(value_r_ - value_l_ + 1));21802181 return distinct_range<T>(value_l_, value_r_).gen_list(size_);2182 }2183};21842185/*******************2186 * *2187 * PERMUTATION *2188 * *2189 *******************/21902191/*2192 * Permutation generation.2193 *2194 * Permutation are defined always as numbers in [0, n), that is, 0-based.2195 */21962197struct permutation : gen_base<permutation> {2198 int size_; // Size of permutation.2199 std::vector<std::pair<int, int>> defs_; // {idx, value}.2200 std::optional<std::vector<int>> cycle_sizes_; // Cycle sizes.22012202 // Creates generator for permutation of size 'size'.2203 permutation(int size) : size_(size) {2204 tgen_ensure(size_ > 0, "permutation: size must be positive");2205 }22062207 // Restricts permutations for permutation[idx] = val.2208 permutation &fix(int idx, int val) {2209 tgen_ensure(0 <= idx and idx < size_,2210 "permutation: index must be valid");2211 defs_.emplace_back(idx, val);2212 return *this;2213 }22142215 // Restricts permutations for permutation to have cycle sizes.2216 permutation &cycles(const std::vector<int> &cycle_sizes) {2217 tgen_ensure(2218 size_ == std::accumulate(cycle_sizes.begin(), cycle_sizes.end(), 0),2219 "permutation: cycle sizes must add up to size of permutation");2220 cycle_sizes_ = cycle_sizes;2221 return *this;2222 }2223 permutation &cycles(const std::initializer_list<int> &cycle_sizes) {2224 return cycles(std::vector<int>(cycle_sizes));2225 }22262227 // Permutation value.2228 // Operations on a value are not random.2229 struct value : gen_value_base<value> {2230 using tgen_is_sequential_tag = detail::is_sequential_tag;22312232 using std_type = std::vector<int>; // std type for value.2233 std::vector<int> vec_; // Permutation.2234 char sep_; // Separator for printing.2235 bool add_1_; // If should add 1, for printing.22362237 value(const std::vector<int> &vec)2238 : vec_(vec), sep_(' '), add_1_(false) {2239 tgen_ensure(!vec_.empty(), "permutation: value: cannot be empty");2240 std::vector<bool> vis(vec_.size(), false);2241 for (int i = 0; i < size(); ++i) {2242 tgen_ensure(0 <= vec_[i] and2243 vec_[i] < static_cast<int>(vec_.size()),2244 "permutation: value: values must be from `0` to "2245 "`size-1`");2246 tgen_ensure(!vis[vec_[i]],2247 "permutation: value: cannot have repeated values");2248 vis[vec_[i]] = true;2249 }2250 }2251 value(const std::initializer_list<int> &il)2252 : value(std::vector<int>(il)) {}22532254 // Fetches size.2255 int size() const { return vec_.size(); }22562257 // Fetches position idx.2258 const int &operator[](int idx) const {2259 tgen_ensure(0 <= idx and idx < size(),2260 "permutation: value: index out of bounds");2261 return vec_[idx];2262 }22632264 // Returns parity of the permutation (+1 if even, -1 if odd).2265 // O(n).2266 int parity() const {2267 std::vector<bool> vis(size(), false);2268 int cycles = 0;22692270 for (int i = 0; i < size(); ++i)2271 if (!vis[i]) {2272 ++cycles;2273 for (int j = i; !vis[j]; j = vec_[j])2274 vis[j] = true;2275 }2276 // Even iff (n - cycles) is even.2277 return ((size() - cycles) % 2 == 0) ? +1 : -1;2278 }22792280 // Sorts values in increasing order.2281 // O(n).2282 value &sort() {2283 for (int i = 0; i < size(); ++i)2284 vec_[i] = i;2285 return *this;2286 }22872288 // Reverses permutation.2289 // O(n).2290 value &reverse() {2291 std::reverse(vec_.begin(), vec_.end());2292 return *this;2293 }22942295 // Inverse of the permutation.2296 // O(n).2297 value &inverse() {2298 std::vector<int> inv(size());2299 for (int i = 0; i < size(); ++i)2300 inv[vec_[i]] = i;2301 swap(vec_, inv);2302 return *this;2303 }23042305 // Sets the separator, for printing.2306 // O(1).2307 value &separator(char sep) {2308 sep_ = sep;2309 return *this;2310 }23112312 // Sets that should print values 1-based.2313 // O(1).2314 value &add_1() {2315 add_1_ = true;2316 return *this;2317 }23182319 // Shuffles permutation uniformly.2320 // O(n).2321 value &shuffle() {2322 for (int i = 0; i < size(); ++i)2323 std::swap(vec_[i], vec_[next(0, size() - 1)]);2324 return *this;2325 }23262327 // Returns a random element uniformly.2328 // O(1).2329 int pick() const { return vec_[next<int>(0, size() - 1)]; }23302331 // Returns vec_[i] with probability proportional to distribution[i].2332 // O(1).2333 template <typename Dist>2334 int pick_by_distribution(const std::vector<Dist> &distribution) const {2335 tgen_ensure(static_cast<size_t>(size()) == distribution.size(),2336 "value and distribution must have the same size");2337 return vec_[next_by_distribution(distribution)];2338 }2339 template <typename Dist>2340 int pick_by_distribution(2341 const std::initializer_list<Dist> &distribution) const {2342 return pick_by_distribution(std::vector<Dist>(distribution));2343 }23442345 // Prints to std::ostream, separated by sep_.2346 friend std::ostream &operator<<(std::ostream &out, const value &val) {2347 for (int i = 0; i < val.size(); ++i) {2348 if (i > 0)2349 out << val.sep_;2350 out << val[i] + val.add_1_;2351 }2352 return out;2353 }23542355 // Gets a std::vector representing the value.2356 std::vector<int> to_std() const { return std_type(vec_); }2357 };23582359 // Generates permutation value.2360 // O(n).2361 value gen() const {2362 if (!cycle_sizes_) {2363 // Cycle sizes not specified.2364 std::vector<int> idx_to_val(size_, -1), val_to_idx(size_, -1);2365 for (auto [idx, val] : defs_) {2366 tgen_ensure(2367 0 <= val and val < size_,2368 "permutation: value in permutation must be in [0, " +2369 std::to_string(size_) + ")");23702371 if (idx_to_val[idx] != -1) {2372 tgen_ensure(idx_to_val[idx] == val,2373 "permutation: cannot set an index to two "2374 "different values");2375 } else2376 idx_to_val[idx] = val;23772378 if (val_to_idx[val] != -1) {2379 tgen_ensure(val_to_idx[val] == idx,2380 "permutation: cannot set two indices to the "2381 "same value");2382 } else2383 val_to_idx[val] = idx;2384 }23852386 std::vector<int> perm(size_);2387 std::iota(perm.begin(), perm.end(), 0);2388 shuffle(perm.begin(), perm.end());2389 int cur_idx = 0;2390 for (int &i : idx_to_val)2391 if (i == -1) {2392 // While this value is used, skip.2393 while (val_to_idx[perm[cur_idx]] != -1)2394 ++cur_idx;2395 i = perm[cur_idx++];2396 }2397 return idx_to_val;2398 }23992400 // Creates cycles.2401 std::vector<int> order(size_);2402 std::iota(order.begin(), order.end(), 0);2403 shuffle(order.begin(), order.end());2404 int idx = 0;2405 std::vector<std::vector<int>> cycles;2406 for (int cycle_size : *cycle_sizes_) {2407 cycles.emplace_back();2408 for (int i = 0; i < cycle_size; ++i)2409 cycles.back().push_back(order[idx++]);2410 }24112412 // Retrieves permutation from cycles.2413 std::vector<int> perm(size_, -1);2414 for (const std::vector<int> &cycle : cycles) {2415 int cur_size = cycle.size();2416 for (int i = 0; i < cur_size; ++i)2417 perm[cycle[i]] = cycle[(i + 1) % cur_size];2418 }24192420 return value(perm);2421 }2422};24232424/************2425 * *2426 * MATH *2427 * *2428 ************/24292430namespace math {24312432namespace detail {24332434using namespace tgen::detail;24352436inline int popcount(uint64_t x) { return __builtin_popcountll(x); }24372438inline int ctzll(uint64_t x) {2439 // Mystery code found on the internet.2440 // Uses de Bruijn sequence.2441 static const unsigned char index64[64] = {2442 0, 1, 2, 53, 3, 7, 54, 27, 4, 38, 41, 8, 34, 55, 48, 28,2443 62, 5, 39, 46, 44, 42, 22, 9, 24, 35, 59, 56, 49, 18, 29, 11,2444 63, 52, 6, 26, 37, 40, 33, 47, 61, 45, 43, 21, 23, 58, 17, 10,2445 51, 25, 36, 32, 60, 20, 57, 16, 50, 31, 19, 15, 30, 14, 13, 12};2446 return index64[((x & -x) * 0x022FDD63CC95386D) >> 58];2447}24482449inline uint64_t mul_mod(uint64_t a, uint64_t b, uint64_t m) {2450 return static_cast<u128>(a) * b % m;2451}24522453// O(log n).2454// 0 <= x < m.2455inline uint64_t expo_mod(uint64_t x, uint64_t y, uint64_t m) {2456 if (!y)2457 return 1;2458 uint64_t ans = expo_mod(mul_mod(x, x, m), y / 2, m);2459 return y % 2 ? mul_mod(x, ans, m) : ans;2460}24612462} // namespace detail24632464// O(log^2 n).2465inline bool is_prime(uint64_t n) {2466 if (n < 2)2467 return false;2468 if (n == 2 or n == 3)2469 return true;2470 if (n % 2 == 0)2471 return false;24722473 uint64_t r = detail::ctzll(n - 1), d = n >> r;2474 // These bases are guaranteed to work for n <= 2^64.2475 for (int a : {2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {2476 uint64_t x = detail::expo_mod(a, d, n);2477 if (x == 1 or x == n - 1 or a % n == 0)2478 continue;24792480 for (uint64_t j = 0; j < r - 1; ++j) {2481 x = detail::mul_mod(x, x, n);2482 if (x == n - 1)2483 break;2484 }2485 if (x != n - 1)2486 return false;2487 }2488 return true;2489}24902491namespace detail {24922493inline uint64_t pollard_rho(uint64_t n) {2494 if (n == 1 or is_prime(n))2495 return n;2496 auto f = [n](uint64_t x) { return mul_mod(x, x, n) + 1; };24972498 uint64_t x = 0, y = 0, t = 30, prd = 2, x0 = 1, q;2499 while (t % 40 != 0 or std::gcd(prd, n) == 1) {2500 if (x == y)2501 x = ++x0, y = f(x);2502 q = mul_mod(prd, x > y ? x - y : y - x, n);2503 if (q != 0)2504 prd = q;2505 x = f(x), y = f(f(y)), ++t;2506 }2507 return std::gcd(prd, n);2508}25092510inline std::vector<uint64_t> factor(uint64_t n) {2511 if (n == 1)2512 return {};2513 if (is_prime(n))2514 return {n};2515 uint64_t d = pollard_rho(n);2516 std::vector<uint64_t> l = factor(d), r = factor(n / d);2517 l.insert(l.end(), r.begin(), r.end());2518 return l;2519}25202521// Error handling.2522template <typename T>2523std::runtime_error there_is_no_in_range_error(const std::string &type, T l,2524 T r) {2525 return error("math: there is no " + type + " in range [" +2526 std::to_string(l) + ", " + std::to_string(r) + "]");2527}2528template <typename T>2529std::runtime_error there_is_no_from_error(const std::string &type, T r) {2530 return error("math: there is no " + type + " from " + std::to_string(r));2531}2532template <typename T>2533std::runtime_error there_is_no_upto_error(const std::string &type, T r) {2534 return error("math: there is no " + type + " up to " + std::to_string(r));2535}25362537// O(log mod).2538// 0 < a < mod.2539// gcd(a, mod) = 1.2540inline i128 modular_inverse_128(i128 a, i128 mod) {2541 tgen_ensure(0 < a and a < mod,2542 "math: modular inverse requires 0 < value < mod");25432544 i128 t = 0, new_t = 1;2545 i128 r = mod, new_r = a;25462547 while (new_r != 0) {2548 i128 q = r / new_r;25492550 auto tmp_t = t - q * new_t;2551 t = new_t;2552 new_t = tmp_t;25532554 auto tmp_r = r - q * new_r;2555 r = new_r;2556 new_r = tmp_r;2557 }25582559 tgen_ensure(r == 1, "math: remainder and mod must be coprime");25602561 if (t < 0)2562 t += mod;2563 return t;2564}25652566// checks if a * b <= limit, for positive numbers.2567inline bool mul_leq(uint64_t a, uint64_t b, uint64_t limit) {2568 if (a == 0 or b == 0)2569 return true;2570 return a <= limit / b;2571}25722573// base^exp, or null if base^exp > limit.2574inline std::optional<uint64_t> expo(uint64_t base, uint64_t exp,2575 uint64_t limit) {2576 uint64_t result = 1;25772578 while (exp) {2579 if (exp & 1) {2580 if (!mul_leq(result, base, limit))2581 return std::nullopt;2582 result *= base;2583 }25842585 exp >>= 1;2586 // Necessary for correctness.2587 if (!exp)2588 break;25892590 if (!mul_leq(base, base, limit))2591 return std::nullopt;2592 base *= base;2593 }2594 return result;2595}25962597// O(log n log k).2598// 0 < k.2599inline uint64_t kth_root_floor(uint64_t n, uint64_t k) {2600 tgen_ensure_against_bug(k > 0, "math: value must be valid");2601 if (k == 1 or n <= 1)2602 return n;26032604 uint64_t lo = 1, hi = 1ULL << ((64 + k - 1) / k);26052606 while (lo < hi) {2607 uint64_t mid = lo + (hi - lo + 1) / 2;26082609 if (expo(mid, k, n)) {2610 lo = mid;2611 } else {2612 hi = mid - 1;2613 }2614 }2615 return lo;2616}26172618// gcd(a, b).2619// O(log a).2620inline i128 gcd128(i128 a, i128 b) {2621 if (a < 0)2622 a = -a;2623 if (b < 0)2624 b = -b;2625 while (b != 0) {2626 i128 t = a % b;2627 a = b;2628 b = t;2629 }2630 return a;2631}26322633// min(2^64, a*b).2634// O(log a).2635// a, b >= 0.2636inline i128 mul_saturate(i128 a, i128 b) {2637 tgen_ensure(a >= 0 and b >= 0);2638 static const i128 LIMIT = static_cast<i128>(1) << 64;2639 if (a == 0 or b == 0)2640 return 0;2641 if (a > LIMIT / b)2642 return LIMIT;2643 return a * b;2644}26452646struct crt {2647 using T = i128;2648 T a, m;26492650 crt() : a(0), m(1) {}2651 crt(T a_, T m_) : a(a_), m(m_) {}2652 crt operator*(crt C) {2653 if (m == 0 or C.m == 0)2654 return {-1, 0};26552656 T g = gcd128(m, C.m);2657 if ((C.a - a) % g != 0)2658 return {-1, 0};26592660 T m1 = m / g;2661 T m2 = C.m / g;26622663 if (m2 == 1)2664 return {a, m};26652666 T inv = modular_inverse_128(m1 % m2, m2);26672668 T k = ((C.a - a) / g) % m2;2669 if (k < 0)2670 k += m2;26712672 k = static_cast<u128>(k) * inv % m2;26732674 T lcm = mul_saturate(m, m2);26752676 T res = (a + static_cast<T>((static_cast<u128>(k) * m) % lcm)) % lcm;2677 if (res < 0)2678 res += lcm;26792680 return {res, lcm};2681 }2682};26832684// Math hacks to operate on log space.26852686inline constexpr long double LOG_ZERO = -INFINITY;2687inline constexpr long double LOG_ONE = 0.0;26882689inline long double log_space(long double x) {2690 return x == 0.0 ? LOG_ZERO : std::log(x);2691}26922693// Math hack to add two values in log space.2694inline long double add_log_space(long double a, long double b) {2695 if (a < b)2696 std::swap(a, b);2697 if (b == LOG_ZERO)2698 return a;2699 return a + log1p(exp(b - a));2700}27012702// Math hack to subtract two values in log space.2703// a >= b.2704inline long double sub_log_space(long double a, long double b) {2705 if (b >= a)2706 return LOG_ZERO;2707 if (b == LOG_ZERO)2708 return a;2709 return a + log1p(-exp(b - a));2710}27112712} // namespace detail27132714// Sorted.2715// O(n^(1/4) log n) expected.2716// 0 < n.2717inline std::vector<uint64_t> factor(uint64_t n) {2718 tgen_ensure(n > 0, "math: number to factor must be positive");2719 auto factors = detail::factor(n);2720 std::sort(factors.begin(), factors.end());2721 return factors;2722}27232724// Sorted.2725// O(n^(1/4) log n) expected.2726// 0 < n.2727inline std::vector<std::pair<uint64_t, int>> factor_by_prime(uint64_t n) {2728 tgen_ensure(n > 0, "math: number to factor must be positive");2729 std::vector<std::pair<uint64_t, int>> primes;2730 for (uint64_t p : factor(n)) {2731 if (!primes.empty() and primes.back().first == p)2732 ++primes.back().second;2733 else2734 primes.emplace_back(p, 1);2735 }2736 return primes;2737}27382739// O(log mod).2740// 0 < a < mod.2741// gcd(a, mod) = 1.2742inline uint64_t modular_inverse(uint64_t a, uint64_t mod) {2743 return detail::modular_inverse_128(a, mod);2744}27452746// O(n^(1/4) log n) expected.2747// 0 < n.2748inline uint64_t totient(uint64_t n) {2749 tgen_ensure(n > 0, "math: totient(0) is undefined");2750 uint64_t phi = n;27512752 for (auto [p, e] : factor_by_prime(n))2753 phi -= phi / p;27542755 return phi;2756}27572758// Returns `(p_i, g_i)`: `p_i` is the prime, `g_i` is the gap.2759inline const std::pair<std::vector<uint64_t>, std::vector<uint64_t>> &2760prime_gaps() {2761 // From https://en.wikipedia.org/wiki/Prime_gap.2762 static const std::pair<std::vector<uint64_t>, std::vector<uint64_t>> value{2763 /* clang-format off */ {2764 2, 3, 7, 23, 89, 113, 523, 887, 1129, 1327, 9551, 15683, 19609,2765 31397, 155921, 360653, 370261, 492113, 1349533, 1357201, 2010733,2766 4652353, 17051707, 20831323, 47326693, 122164747, 189695659,2767 191912783, 387096133, 436273009, 1294268491, 1453168141,2768 2300942549, 3842610773, 4302407359, 10726904659, 20678048297,2769 22367084959, 25056082087, 42652618343, 127976334671, 182226896239,2770 241160624143, 297501075799, 303371455241, 304599508537,2771 416608695821, 461690510011, 614487453523, 738832927927,2772 1346294310749, 1408695493609, 1968188556461, 2614941710599,2773 7177162611713, 13829048559701, 19581334192423, 42842283925351,2774 90874329411493, 171231342420521, 218209405436543, 1189459969825483,2775 1686994940955803, 1693182318746371, 43841547845541059,2776 55350776431903243, 80873624627234849, 203986478517455989,2777 218034721194214273, 305405826521087869, 352521223451364323,2778 401429925999153707, 418032645936712127, 804212830686677669,2779 1425172824437699411, 5733241593241196731, 67879889996577777972780 }, /* clang-format on */2781 {1, 2, 4, 6, 8, 14, 18, 20, 22, 34, 36,2782 44, 52, 72, 86, 96, 112, 114, 118, 132, 148, 154,2783 180, 210, 220, 222, 234, 248, 250, 282, 288, 292, 320,2784 336, 354, 382, 384, 394, 456, 464, 468, 474, 486, 490,2785 500, 514, 516, 532, 534, 540, 582, 588, 602, 652, 674,2786 716, 766, 778, 804, 806, 906, 916, 924, 1132, 1184, 1198,2787 1220, 1224, 1248, 1272, 1328, 1356, 1370, 1442, 1476, 1488, 1510}};27882789 return value;2790}27912792// Returns pair (first_composite_in_gap, last_composite_in_gap).2793// O(log(right)) approximately.2794inline std::pair<uint64_t, uint64_t> prime_gap_upto(uint64_t right) {2795 if (right < 4)2796 throw detail::there_is_no_upto_error("prime gap", right);27972798 const auto &[P, G] = prime_gaps();2799 for (int i = P.size() - 1;; --i) {2800 if (P[i] >= right)2801 continue;28022803 uint64_t real_right = std::min(right, P[i] + G[i] - 1);2804 uint64_t prev = i > 0 ? G[i - 1] : 0;2805 uint64_t curr = real_right - P[i];28062807 if (curr >= prev)2808 return {P[i] + 1, real_right};2809 }2810}28112812// From https://oeis.org/A002182/b002182.txt.2813inline const std::vector<uint64_t> &highly_composites() {2814 /* clang-format off */2815 static const std::vector<uint64_t> highly_composites = {2816 1, 2, 4, 6, 12, 24, 36, 48, 60, 120, 180, 240, 360, 720, 840, 1260, 1680,2817 2520, 5040, 7560, 10080, 15120, 20160, 25200, 27720, 45360, 50400, 55440,2818 83160, 110880, 166320, 221760, 277200, 332640, 498960, 554400, 665280,2819 720720, 1081080, 1441440, 2162160, 2882880, 3603600, 4324320, 6486480,2820 7207200, 8648640, 10810800, 14414400, 17297280, 21621600, 32432400,2821 36756720, 43243200, 61261200, 73513440, 110270160, 122522400, 147026880,2822 183783600, 245044800, 294053760, 367567200, 551350800, 698377680, 735134400,2823 1102701600, 1396755360, 2095133040, 2205403200, 2327925600, 2793510720,2824 3491888400, 4655851200, 5587021440, 6983776800, 10475665200, 13967553600,2825 20951330400, 27935107200, 41902660800, 48886437600, 64250746560,2826 73329656400, 80313433200, 97772875200, 128501493120, 146659312800,2827 160626866400, 240940299600, 293318625600, 321253732800, 481880599200,2828 642507465600, 963761198400, 1124388064800, 1606268664000, 1686582097200,2829 1927522396800, 2248776129600, 3212537328000, 3373164194400, 4497552259200,2830 6746328388800, 8995104518400, 9316358251200, 13492656777600, 18632716502400,2831 26985313555200, 27949074753600, 32607253879200, 46581791256000,2832 48910880818800, 55898149507200, 65214507758400, 93163582512000,2833 97821761637600, 130429015516800, 195643523275200, 260858031033600,2834 288807105787200, 391287046550400, 577614211574400, 782574093100800,2835 866421317361600, 1010824870255200, 1444035528936000, 1516237305382800,2836 1732842634723200, 2021649740510400, 2888071057872000, 3032474610765600,2837 4043299481020800, 6064949221531200, 8086598962041600, 10108248702552000,2838 12129898443062400, 18194847664593600, 20216497405104000, 24259796886124800,2839 30324746107656000, 36389695329187200, 48519593772249600, 60649492215312000,2840 72779390658374400, 74801040398884800, 106858629141264000,2841 112201560598327200, 149602080797769600, 224403121196654400,2842 299204161595539200, 374005201994424000, 448806242393308800,2843 673209363589963200, 748010403988848000, 897612484786617600,2844 1122015605983272000, 1346418727179926400, 1795224969573235200,2845 2244031211966544000, 2692837454359852800, 3066842656354276800,2846 4381203794791824000, 4488062423933088000, 6133685312708553600,2847 8976124847866176000, 9200527969062830400, 12267370625417107200ULL,2848 15334213281771384000ULL, 18401055938125660800ULL}; /* clang-format on */2849 return highly_composites;2850}28512852// O(log(right)) approximately.2853inline uint64_t highly_composite_upto(uint64_t right) {2854 for (int i = highly_composites().size() - 1; i >= 0; --i)2855 if (highly_composites()[i] <= right)2856 return highly_composites()[i];28572858 throw detail::there_is_no_upto_error("highly composite number", right);2859}28602861// O(log^3 (right)) expected.2862// Generates a random prime in [left, right].2863inline uint64_t gen_prime(uint64_t left, uint64_t right) {2864 if (right < left or right < 2)2865 throw detail::there_is_no_in_range_error("prime", left, right);2866 left = std::max<uint64_t>(left, 2);2867 auto [l_gap, r_gap] = prime_gap_upto(right);2868 if (right - left + 1 <= r_gap - l_gap + 1) {2869 // There might be no primes in the range.2870 std::vector<uint64_t> vals(right - left + 1);2871 iota(vals.begin(), vals.end(), left);2872 shuffle(vals.begin(), vals.end());2873 for (uint64_t i : vals)2874 if (is_prime(i))2875 return i;2876 throw detail::there_is_no_in_range_error("prime", left, right);2877 }28782879 uint64_t n;2880 do {2881 n = next(left, right);2882 } while (!is_prime(n));2883 return n;2884}28852886// O(log^3 (left)) expected.2887// left <= 2^64 - 59.2888inline uint64_t prime_from(uint64_t left) {2889 tgen_ensure(left <= std::numeric_limits<uint64_t>::max() - 58,2890 "math: invalid bound");2891 for (uint64_t i = std::max<uint64_t>(2, left);; ++i)2892 if (is_prime(i))2893 return i;2894}28952896// O(log^3 (right)) expected.2897inline uint64_t prime_upto(uint64_t right) {2898 if (right >= 2)2899 for (uint64_t i = right; i >= 2; --i)2900 if (is_prime(i))2901 return i;2902 throw detail::there_is_no_upto_error("prime", right);2903}29042905// O(n^(1/4) log n) expected.2906// 0 < n.2907inline int num_divisors(uint64_t n) {2908 int divisors = 1;2909 for (auto [p, e] : factor_by_prime(n))2910 divisors *= (e + 1);2911 return divisors;2912}29132914// Random number in [left, right] with `divisor_count` divisors.2915// O(log(right) log(divisor_count)).2916// divisor_count must be prime.2917inline uint64_t gen_divisor_count(uint64_t left, uint64_t right,2918 int divisor_count) {2919 tgen_ensure(divisor_count > 0 and is_prime(divisor_count),2920 "math: divisor count must be prime");2921 int root = divisor_count - 1;2922 uint64_t lo = detail::kth_root_floor(left, root);2923 if (*detail::expo(lo, root, left) < left)2924 ++lo;2925 uint64_t p = gen_prime(lo, detail::kth_root_floor(right, root));2926 return *detail::expo(p, root, right);2927}29282929// O(|mods| + log (right)).2930// |rems| = |mods|.2931// rems_i < mods_i.2932inline uint64_t gen_congruent(uint64_t left, uint64_t right,2933 std::vector<uint64_t> rems,2934 std::vector<uint64_t> mods) {2935 if (left > right)2936 throw detail::there_is_no_in_range_error("congruent number", left,2937 right);2938 tgen_ensure(rems.size() == mods.size(),2939 "math: number of remainders and mods must be the same");2940 tgen_ensure(rems.size() > 0, "math: must have at least one congruence");29412942 detail::crt crt;2943 for (int i = 0; i < static_cast<int>(rems.size()); ++i) {2944 tgen_ensure(rems[i] < mods[i],2945 "math: remainder must be smaller than the mod");2946 crt = crt * detail::crt(rems[i], mods[i]);29472948 if (crt.a == -1)2949 throw detail::there_is_no_in_range_error("congruent number", left,2950 right);2951 if (crt.m > right) {2952 if (!(left <= crt.a and crt.a <= right))2953 throw detail::there_is_no_in_range_error("congruent number",2954 left, right);29552956 for (int j = 0; j < static_cast<int>(rems.size()); ++j)2957 if (crt.a % mods[j] != rems[j])2958 throw detail::there_is_no_in_range_error("congruent number",2959 left, right);2960 return crt.a;2961 }2962 }29632964 uint64_t k_min = crt.a >= left ? 0 : ((left - crt.a) + crt.m - 1) / crt.m;2965 uint64_t k_max = (right - crt.a) / crt.m;29662967 if (k_min > k_max)2968 throw detail::there_is_no_in_range_error("congruent number", left,2969 right);29702971 return crt.a + next(k_min, k_max) * crt.m;2972}29732974// O(log (right)).2975// rem < mod.2976inline uint64_t gen_congruent(uint64_t left, uint64_t right, uint64_t rem,2977 uint64_t mod) {2978 return gen_congruent(left, right, std::vector<uint64_t>({rem}),2979 std::vector<uint64_t>({mod}));2980}29812982// First congruent number >= left.2983// O(|mods| + log (left)).2984// |rems| = |mods|.2985// rems_i < mods_i.2986inline uint64_t congruent_from(uint64_t left, std::vector<uint64_t> rems,2987 std::vector<uint64_t> mods) {2988 tgen_ensure(rems.size() == mods.size(),2989 "math: number of remainders and mods must be the same");2990 tgen_ensure(rems.size() > 0, "math: must have at least one congruence");29912992 detail::crt crt;2993 for (int i = 0; i < static_cast<int>(rems.size()); ++i) {2994 tgen_ensure(rems[i] < mods[i],2995 "math: remainder must be smaller than the mod");2996 crt = crt * detail::crt(rems[i], mods[i]);29972998 if (crt.a == -1)2999 throw detail::there_is_no_from_error("congruent number", left);3000 if (crt.m > std::numeric_limits<uint64_t>::max()) {3001 if (crt.a < left)3002 throw detail::error(3003 "math: congruent number does not exist or is too large");30043005 for (int j = 0; j < static_cast<int>(rems.size()); ++j)3006 if (crt.a % mods[j] != rems[j])3007 throw detail::error("math: congruent number does "3008 "not exist or is too large");3009 return crt.a;3010 }3011 }30123013 uint64_t k = 0;3014 if (crt.a < left)3015 k = ((left - crt.a) + crt.m - 1) / crt.m;3016 detail::i128 result = crt.a + k * crt.m;30173018 if (result > std::numeric_limits<uint64_t>::max())3019 throw detail::error("math: congruent number is too large");3020 return result;3021}30223023// O(log (left))3024// rem < mod.3025inline uint64_t congruent_from(uint64_t left, uint64_t rem, uint64_t mod) {3026 return congruent_from(left, std::vector<uint64_t>{rem},3027 std::vector<uint64_t>{mod});3028}30293030// Last congruent number <= right.3031// O(|mods| + log (right)).3032// |rems| = |mods|.3033// rems_i < mods_i.3034inline uint64_t congruent_upto(uint64_t right, std::vector<uint64_t> rems,3035 std::vector<uint64_t> mods) {3036 tgen_ensure(rems.size() == mods.size(),3037 "math: number of remainders and mods must be the same");3038 tgen_ensure(rems.size() > 0, "math: must have at least one congruence");30393040 detail::crt crt;3041 for (int i = 0; i < static_cast<int>(rems.size()); ++i) {3042 tgen_ensure(rems[i] < mods[i],3043 "math: remainder must be smaller than the mod");30443045 crt = crt * detail::crt(rems[i], mods[i]);30463047 if (crt.a == -1)3048 throw detail::there_is_no_upto_error("congruent number", right);3049 if (crt.m > right) {3050 if (!(crt.a <= right))3051 throw detail::there_is_no_upto_error("congruent number", right);30523053 for (int j = 0; j < static_cast<int>(rems.size()); ++j)3054 if (crt.a % mods[j] != rems[j])3055 throw detail::there_is_no_upto_error("congruent number",3056 right);3057 return crt.a;3058 }3059 }30603061 if (crt.a > right)3062 throw detail::there_is_no_upto_error("congruent number", right);30633064 uint64_t k = (right - crt.a) / crt.m;3065 detail::i128 result = crt.a + k * crt.m;30663067 if (result < 0)3068 throw detail::there_is_no_upto_error("congruent number", right);3069 return result;3070}30713072// O(log r)3073// rem < mod.3074inline uint64_t congruent_upto(uint64_t right, uint64_t rem, uint64_t mod) {3075 return congruent_upto(right, std::vector<uint64_t>{rem},3076 std::vector<uint64_t>{mod});3077}30783079// Mod used for FFT/NTT.3080inline constexpr int FFT_MOD = 998244353;30813082// Fibonacci sequence up to 2^64.3083inline const std::vector<uint64_t> &fibonacci() {3084 static const std::vector<uint64_t> fib = [] {3085 std::vector<uint64_t> v = {0, 1};3086 while (v.back() <=3087 std::numeric_limits<uint64_t>::max() - v[v.size() - 2])3088 v.push_back(v.back() + v[v.size() - 2]);3089 return v;3090 }();3091 return fib;3092}30933094// Partition is ordered (composition), that is, (1, 1, 2) != (1, 2, 1).3095// O(n).3096// 0 < n.3097// 0 < part_left.3098inline std::vector<int>3099gen_partition(int n, int part_left = 1,3100 std::optional<int> part_right = std::nullopt) {3101 if (!part_right.has_value())3102 part_right = n;3103 part_right = std::min(*part_right, n);3104 tgen_ensure(n > 0 and part_left > 0,3105 "math: invalid parameters to gen_partition");3106 tgen_ensure(part_left <= n and *part_right > 0, "math: no such partition");31073108 // dp[i] = log(number of ways to add to i).3109 std::vector<long double> dp(n + 1, detail::LOG_ZERO);3110 dp[0] = detail::LOG_ONE;3111 long double window = detail::LOG_ZERO;3112 for (int i = 1; i <= n; ++i) {3113 if (i >= part_left)3114 window = detail::add_log_space(window, dp[i - part_left]);3115 if (i >= *part_right + 1)3116 window = detail::sub_log_space(window, dp[i - *part_right - 1]);3117 dp[i] = window;3118 }3119 tgen_ensure(dp[n] >= 0, "math: no such partition");31203121 // Crazy math tricks ahead.3122 auto dp_pref = dp;3123 for (int i = 1; i <= n; ++i)3124 dp_pref[i] = detail::add_log_space(dp_pref[i - 1], dp[i]);31253126 std::vector<int> part;3127 int sum = n;3128 while (sum > 0) {3129 // Will generate a number such that what remains is in [l, r].3130 int l = std::max(0, sum - *part_right), r = sum - part_left;3131 detail::tgen_ensure_against_bug(r >= 0, "math: r < 0 in gen_partition");31323133 int nxt_sum = std::min(sum, r);3134 long double random = next<long double>(0, 1);31353136 // We generate a value X (log space), and then choose nxt_sum such3137 // that dp_pref[nxt_sum-1] < X <= dp_pref[nxt_sum].31383139 // Math hack:3140 // Let A = pref[l-1], B = pref[r], U = rand().3141 // X = log[exp(A) + U * (exp(B) - exp(A))]3142 // = log{exp(B) * [exp(A) / exp(B) + U * (1 - exp(A) / exp(B))]}3143 // = B + log[exp(A - B) + U - U * exp(A - B))]3144 // = B + log[U + (1 - U) * exp(A - B)].3145 long double val_l = l ? dp_pref[l - 1] : detail::LOG_ZERO,3146 val_r = dp_pref[r];3147 while (nxt_sum > l and3148 dp_pref[nxt_sum - 1] >=3149 val_r + detail::log_space(random +3150 (1 - random) * exp(val_l - val_r)))3151 --nxt_sum;31523153 part.push_back(sum - nxt_sum);3154 sum = nxt_sum;3155 }31563157 return part;3158}31593160// Partition is ordered (composition), that is, (1, 1, 2) != (1, 2, 1).3161// O(n) time/memory if part_right is not set, O(n * k) time/memory otherwise.3162// 0 < k <= n.3163// 0 <= part_left.3164inline std::vector<int>3165gen_partition_fixed_size(int n, int k, int part_left = 0,3166 std::optional<int> part_right = std::nullopt) {3167 if (!part_right.has_value())3168 part_right = n;3169 part_right = std::min(*part_right, n);3170 tgen_ensure(0 < k and k <= n and part_left >= 0,3171 "math: invalid parameters to gen_partition_fixed_size");3172 tgen_ensure(static_cast<long long>(k) * part_left <= n and3173 n <= static_cast<long long>(k) * (*part_right),3174 "math: no such partition");31753176 // What we need to distribute to the parts.3177 int s = n - k * part_left;31783179 std::vector<int> part(k);3180 if (*part_right == n) {3181 // Stars and bars - O(n).3182 std::vector<int> cuts = {-1};31833184 int total = s + k - 1, bars = k - 1;3185 for (int i = 0; i < total and bars > 0; ++i)3186 if (next<long double>(0, 1) <3187 static_cast<long double>(bars) / (total - i)) {3188 cuts.push_back(i);3189 --bars;3190 }3191 cuts.push_back(total);31923193 // Recovers parts.3194 for (int i = 0; i < k; ++i)3195 part[i] = cuts[i + 1] - cuts[i] - 1;3196 } else {3197 // DP with log trick - O(nk).3198 int u = *part_right - part_left;31993200 // dp[i][j] = log(#ways to fill i parts with sum j)3201 std::vector<std::vector<long double>> dp(3202 k + 1, std::vector<long double>(s + 1, detail::LOG_ZERO));3203 dp[0][0] = detail::LOG_ONE;32043205 for (int i = 1; i <= k; ++i) {3206 std::vector<long double> pref = dp[i - 1];3207 for (int j = 1; j <= s; ++j)3208 pref[j] = detail::add_log_space(pref[j - 1], dp[i - 1][j]);32093210 for (int j = 0; j <= s; ++j) {3211 dp[i][j] = pref[j];3212 if (j >= u + 1)3213 dp[i][j] = detail::sub_log_space(dp[i][j], pref[j - u - 1]);3214 }3215 }32163217 // Recovers parts backwards.3218 int left_to_distribute = s;3219 for (int i = k; i >= 1; --i) {3220 long double log_total = detail::LOG_ZERO;3221 for (int j = 0; j <= u and j <= left_to_distribute; ++j)3222 log_total = detail::add_log_space(3223 log_total, dp[i - 1][left_to_distribute - j]);3224 detail::tgen_ensure_against_bug(3225 log_total != detail::LOG_ZERO,3226 "math: total == 0 in gen_partition_fixed_size");32273228 // Now we choose a number with probability proportional to3229 // dp[i-1][.].32303231 // log(rand() * total) = log(rand()) + log(total).3232 long double random =3233 detail::log_space(next<long double>(0, 1)) + log_total;32343235 long double cur_prob = detail::LOG_ZERO;3236 int chosen = 0;3237 for (int j = 0; j <= u and j <= left_to_distribute; ++j) {3238 cur_prob = detail::add_log_space(3239 cur_prob, dp[i - 1][left_to_distribute - j]);3240 if (random < cur_prob) {3241 chosen = j;3242 break;3243 }3244 }32453246 part[k - i] = chosen;3247 left_to_distribute -= chosen;3248 }3249 }32503251 for (int &i : part)3252 i += part_left;3253 return part;3254}32553256// Partition is ordered (composition), that is, (1, 1, 2) != (1, 2, 1).3257// Inspired by jngen rndm.partition: random delimiters, sort, gap recovery;3258// omits jngen's part reordering, shuffles, and two-pass redistribution.3259// 0 < k <= n.3260// 0 <= part_left.3261// Not uniformly random; optimized for speed.3262// O(k log k).3263inline std::vector<uint64_t> gen_partition_fixed_size_fast(3264 uint64_t n, int k, uint64_t part_left = 0,3265 std::optional<uint64_t> part_right = std::nullopt) {3266 if (!part_right.has_value())3267 part_right = n;3268 part_right = std::min(*part_right, n);32693270 detail::u128 n128 = n;3271 detail::u128 k128 = k;3272 detail::u128 part_left128 = part_left;3273 detail::u128 part_right128 = *part_right;32743275 tgen_ensure(k > 0 and k128 <= n128,3276 "math: invalid parameters to gen_partition_fixed_size_fast");3277 tgen_ensure(part_right128 >= part_left128 and3278 k128 * part_left128 <= n128 and3279 k128 * part_right128 >= n128,3280 "math: no such partition");32813282 uint64_t slack_total = n128 - k128 * part_left128;3283 uint64_t slack_max = part_right128 - part_left128;32843285 std::vector<uint64_t> part(k);3286 if (k == 1) {3287 part[0] = slack_total;3288 } else {3289 std::vector<uint64_t> cuts(k - 1);3290 for (uint64_t &d : cuts)3291 d = next<uint64_t>(0, slack_total);3292 std::sort(cuts.begin(), cuts.end());32933294 uint64_t prev = 0;3295 for (int i = 0; i + 1 < k; ++i) {3296 part[i] = cuts[i] - prev;3297 prev = cuts[i];3298 }3299 part[k - 1] = slack_total - prev;3300 }33013302 auto add_part_left = [part_left](uint64_t x) -> uint64_t {3303 detail::u128 val = x + part_left;3304 detail::tgen_ensure_against_bug(3305 val <= std::numeric_limits<uint64_t>::max(),3306 "math: part + part_left exceeds uint64_t in "3307 "gen_partition_fixed_size_fast");3308 return val;3309 };33103311 if (slack_max >= slack_total) {3312 for (uint64_t &x : part)3313 x = add_part_left(x);3314 return part;3315 }33163317 detail::u128 remaining = 0;3318 for (uint64_t &x : part) {3319 if (x > slack_max) {3320 remaining += x - slack_max;3321 x = slack_max;3322 }3323 x = add_part_left(x);3324 }33253326 if (remaining > 0) {3327 for (uint64_t &x : part) {3328 if (x < *part_right && remaining > 0) {3329 detail::u128 room = *part_right - x;3330 detail::u128 add = std::min(remaining, room);3331 detail::u128 val = x + add;3332 detail::tgen_ensure_against_bug(3333 val <= *part_right,3334 "math: part exceeds part_right after redistribution in "3335 "gen_partition_fixed_size_fast");3336 x = val;3337 remaining -= add;3338 }3339 }3340 detail::tgen_ensure_against_bug(3341 remaining == 0, "math: remaining mass after redistribution in "3342 "gen_partition_fixed_size_fast");3343 }33443345 return part;3346}33473348// Random partition of elements into k ordered groups (input order preserved).3349// If max_size is unset, part sizes are uniform via gen_partition_fixed_size.3350// If max_size is set, uses gen_partition_fixed_size_fast (not uniform).3351// O(n) if max_size is unset; O(n + k log k) if max_size is set.3352template <typename T>3353std::vector<std::vector<T>>3354partition_elements(std::vector<T> elements, int k, int min_size = 0,3355 std::optional<uint64_t> max_size = std::nullopt) {3356 size_t n = elements.size();3357 tgen_ensure(k > 0, "math: partition_elements: k must be positive");3358 tgen_ensure(min_size >= 0,3359 "math: partition_elements: min_size must be non-negative");33603361 std::vector<uint64_t> sizes;3362 if (max_size.has_value()) {3363 sizes = gen_partition_fixed_size_fast(n, k, min_size, max_size);3364 } else {3365 for (int sz : gen_partition_fixed_size(n, k, min_size))3366 sizes.push_back(sz);3367 }33683369 std::vector<std::vector<T>> groups;3370 groups.reserve(k);3371 size_t pos = 0;3372 for (uint64_t sz : sizes) {3373 groups.emplace_back(elements.begin() + pos,3374 elements.begin() + pos + sz);3375 pos += sz;3376 }3377 return groups;3378}33793380}; // namespace math33813382/**************3383 * *3384 * STRING *3385 * *3386 **************/33873388namespace detail {33893390/*3391 * Regex.3392 *3393 * Compatible with testlib's regex.3394 *3395 * Operations:3396 * - A single character yields itself ("a", "3").3397 * - A list of characters inside square braces yields any a random element3398 * from the list ("[abc123]").3399 * - A range of characters is equivalent to listing them ("[a-z1-9A-Z]").3400 * - A pattern followed by {n} yields the pattern repeated n times ("a{3}").3401 * - A pattern followed by {l,r} yields the pattern repeated between l and r3402 * times, uniformly at random ("a{3,5}").3403 * - A list of patterns separated by | yields a random pattern from the3404 * list, uniformly at random ("abc|def|ghi").3405 * - Parentheses can be used for grouping ("a((a|b){3})").3406 *3407 * Examples:3408 * 1. str("[1-9][0-9]{1,2}") generates two- or three-digit numbers.3409 * 2. str("a[b-d]{2}|e") generates "e" or a random string of length 3, with3410 * the first character being 'a' and the second and3411 * third characters being 'b', 'c', or 'd'.3412 * 3. str("[1-9][0-9]{%d}", n-1) generates n-digit numbers.3413 *3414 * Operations defined by {n} and {l,r} are applied from left to right, and3415 * the pattern that comes before has its delimiters defined either by () or3416 * [] at its end or is taken from the beginning of the pattern (in3417 * "a[bc]{2}", "{2}" is applied to "[bc]", and in "[01]abc{3}", the "{3}" is3418 * applied to "[01]abc").3419 */34203421// If it has children, it is either a SEQ or an OR group, defined by the3422// pattern_ field.3423struct regex_node {3424 // Considered to be repetition of left_bound != -1, pattern if3425 // children_.empty(), otherwise "SEQ" or "OR", defined by the pattern_3426 // field.3427 std::string3428 pattern_; // Either pattern, or "SEQ" or "OR" (if !children_.empty()).3429 std::vector<regex_node> children_; // Children, when SEQ or OR.3430 int left_bound_, right_bound_; // Left and right bounds of the repetition,3431 // or -1 if not a repetition.3432 double3433 log_space_num_ways_; // Log space number of ways to match the pattern.3434 std::optional<distinct_container<char>>3435 distinct_; // Distinct generator for the pattern, for [chars].34363437 // c or [chars].3438 regex_node(const std::string &pattern)3439 : pattern_(pattern), left_bound_(-1), right_bound_(-1) {3440 if (pattern.size() == 1) {3441 log_space_num_ways_ = math::detail::LOG_ONE;3442 return;3443 }3444 tgen_ensure_against_bug(pattern[0] == '[' and pattern.back() == ']',3445 "str: invalid regex: expected character class");3446 int size = pattern.size() - 2;3447 log_space_num_ways_ = math::detail::log_space(size);3448 distinct_ = distinct_container<char>(pattern.substr(1, size));3449 }3450 // SEQ or OR.3451 regex_node(const std::string &pattern, std::vector<regex_node> &children)3452 : pattern_(pattern), left_bound_(-1), right_bound_(-1) {3453 if (pattern == "SEQ") {3454 // Multiply the number of ways.3455 log_space_num_ways_ = math::detail::LOG_ONE;3456 for (const auto &child : children)3457 log_space_num_ways_ += child.log_space_num_ways_;3458 } else if (pattern == "OR") {3459 // Add the number of ways.3460 log_space_num_ways_ = math::detail::LOG_ZERO;3461 for (const auto &child : children)3462 log_space_num_ways_ = math::detail::add_log_space(3463 log_space_num_ways_, child.log_space_num_ways_);3464 } else3465 tgen_ensure_against_bug("str: invalid regex: expected SEQ or OR");34663467 children_ = std::move(children);3468 children.clear();3469 }3470 // REP.3471 regex_node(int left_bound, int right_bound, regex_node &child)3472 : pattern_("REP"), left_bound_(left_bound), right_bound_(right_bound) {3473 log_space_num_ways_ = math::detail::LOG_ZERO;3474 for (int i = left_bound; i <= right_bound; ++i)3475 log_space_num_ways_ = math::detail::add_log_space(3476 log_space_num_ways_, i * child.log_space_num_ways_);34773478 children_.push_back(std::move(child));3479 }3480};34813482// State of the regex parser.3483struct regex_state {3484 std::vector<regex_node> cur; // Current sequence of nodes.3485 std::vector<regex_node> branches; // Branches of the current OR group.3486};34873488// Creates a SEQ node from the current state.3489inline regex_node make_regex_seq(regex_state &st) {3490 return regex_node("SEQ", st.cur);3491}34923493// Finishes current state.3494inline regex_node finish_regex_state(regex_state &st) {3495 // SEQ.3496 if (st.branches.empty())3497 return make_regex_seq(st);34983499 // OR.3500 st.branches.push_back(make_regex_seq(st));3501 return regex_node("OR", st.branches);3502}35033504// Parses a regex pattern into a tree, computing the number of ways to match the3505// pattern.3506inline regex_node parse_regex(std::string regex) {3507 std::string new_regex;3508 for (char c : regex)3509 if (c != ' ')3510 new_regex += c;3511 swap(regex, new_regex);3512 regex_state cur;3513 std::vector<regex_state> stack;35143515 for (size_t i = 0; i < regex.size(); ++i) {3516 char c = regex[i];35173518 if (c == '(') {3519 // Pushes the current state to the stack.3520 stack.push_back(std::move(cur));3521 cur = regex_state();3522 } else if (c == ')') {3523 // Finishes the current state, and adds it to the parent.3524 regex_node node = finish_regex_state(cur);35253526 tgen_ensure(!stack.empty(), "str: invalid regex: unmatched `)`");3527 cur = std::move(stack.back());3528 stack.pop_back();35293530 cur.cur.push_back(std::move(node));3531 } else if (c == '|') {3532 // Starts a new OR group.3533 regex_node node = make_regex_seq(cur);3534 cur.branches.push_back(std::move(node));3535 } else if (c == '[') {3536 // Parses a character class.3537 std::string chars;35383539 for (++i; i < regex.size() and regex[i] != ']'; ++i) {3540 if (i + 2 < regex.size() and regex[i + 1] == '-') {3541 char a = regex[i], b = regex[i + 2];3542 if (a > b)3543 std::swap(a, b);3544 for (char x = a; x <= b; ++x)3545 chars += x;3546 i += 2;3547 } else3548 chars += regex[i];3549 }35503551 tgen_ensure(i < regex.size() and regex[i] == ']',3552 "str: invalid regex: unmatched `[`");3553 cur.cur.emplace_back("[" + chars + "]");3554 } else if (c == '{') {3555 // Parses a repetition.3556 ++i;3557 int l = -1, r = -1;35583559 while (i < regex.size() and3560 isdigit(static_cast<unsigned char>(regex[i]))) {3561 if (l == -1)3562 l = 0;3563 tgen_ensure(l <= static_cast<int>(1e8),3564 "str: invalid regex: number too large inside `{}`");3565 l = 10 * l + (regex[i] - '0');3566 ++i;3567 }35683569 if (i < regex.size() and regex[i] == ',') {3570 ++i;3571 while (i < regex.size() and3572 isdigit(static_cast<unsigned char>(regex[i]))) {3573 if (r == -1)3574 r = 0;3575 tgen_ensure(3576 r <= static_cast<int>(1e8),3577 "str: invalid regex: number too large inside `{}`");3578 r = 10 * r + (regex[i] - '0');3579 ++i;3580 }3581 } else3582 r = l;35833584 tgen_ensure(i < regex.size() and regex[i] == '}',3585 "str: invalid regex: unmatched `{`");3586 tgen_ensure(l != -1 and r != -1,3587 "str: invalid regex: missing number inside `{}`");3588 tgen_ensure(l <= r,3589 "str: invalid regex: invalid range inside `{}`");35903591 // Creates a REP node from the previous node.3592 tgen_ensure(!cur.cur.empty(),3593 "str: invalid regex: expected expression before `{}`");35943595 regex_node rep(l, r, cur.cur.back());3596 cur.cur.pop_back();3597 cur.cur.push_back(std::move(rep));3598 } else {3599 // Creates a char node.3600 cur.cur.emplace_back(std::string(1, c));3601 }3602 }36033604 tgen_ensure(stack.empty(), "str: invalid regex: unmatched `(`");3605 return finish_regex_state(cur);3606}36073608// Generates a uniformly random string that matches the given regex.3609inline void gen_regex(const regex_node &node, std::string &str) {3610 // For [chars], generate a random character from the list.3611 if (node.pattern_[0] == '[') {3612 str += node.pattern_[1 + next<int>(0, node.pattern_.size() - 3)];3613 return;3614 }36153616 // For REP, generate a random number of times to repeat the pattern.3617 if (node.left_bound_ != -1) {3618 // Generates a random value W from 0 to num_ways.3619 // log(W) = log(random(0, 1) * num_ways)3620 // = log(random(0, 1)) + log(num_ways).3621 double log_rand = math::detail::log_space(next<double>(0, 1)) +3622 node.log_space_num_ways_;3623 double cur_prob = math::detail::LOG_ZERO;3624 double child_num_ways = node.children_[0].log_space_num_ways_;36253626 for (int i = node.left_bound_; i <= node.right_bound_; ++i) {3627 cur_prob =3628 math::detail::add_log_space(cur_prob, i * child_num_ways);3629 if (log_rand <= cur_prob) {3630 for (int j = 0; j < i; ++j)3631 gen_regex(node.children_[0], str);3632 return;3633 }3634 }36353636 tgen_ensure_against_bug(false,3637 "str: log_rand > cur_prob in REP gen_regex");3638 }36393640 // For SEQ, generate all children.3641 if (!node.children_.empty() and node.pattern_ == "SEQ") {3642 for (const regex_node &child : node.children_)3643 gen_regex(child, str);3644 return;3645 }36463647 // For OR, generate a random child.3648 if (!node.children_.empty() and node.pattern_ == "OR") {3649 // Generates a random value W from 0 to num_ways.3650 // log(W) = log(random(0, 1) * num_ways)3651 // = log(random(0, 1)) + log(num_ways).3652 double log_rand = math::detail::log_space(next<double>(0, 1)) +3653 node.log_space_num_ways_;3654 double cur_prob = math::detail::LOG_ZERO;36553656 for (const regex_node &child : node.children_) {3657 cur_prob = math::detail::add_log_space(cur_prob,3658 child.log_space_num_ways_);3659 if (log_rand <= cur_prob) {3660 gen_regex(child, str);3661 return;3662 }3663 }36643665 tgen_ensure_against_bug(false,3666 "str: log_rand > cur_prob in OR gen_regex");3667 }36683669 // For char, generate the character.3670 detail::tgen_ensure_against_bug(3671 node.pattern_.size() == 1,3672 "str: invalid regex: expected single character, but got `" +3673 node.pattern_ + "`");3674 str += node.pattern_[0];3675}36763677// Formats a regex string with given arguments.3678template <typename... Args>3679std::string regex_format(const std::string &s, Args &&...args) {3680 if constexpr (sizeof...(Args) == 0) {3681 return s;3682 } else {3683 int size = std::snprintf(nullptr, 0, s.c_str(), args...) + 1;3684 std::string buf(size, '\0');3685 std::snprintf(buf.data(), size, s.c_str(), args...);3686 buf.pop_back(); // remove '\0'3687 return buf;3688 }3689}36903691} // namespace detail36923693/*3694 * String generator.3695 */36963697struct str : gen_base<str> {3698 std::optional<list<char>> list_; // List of characters.3699 std::optional<detail::regex_node>3700 root_; // Root node of the regex tree for the whole string.37013702 // Creates generator for strings of size 'size', with random characters in3703 // [value_left, value_right].3704 str(int size, char value_left = 'a', char value_right = 'z') {3705 tgen_ensure(size > 0, "str: size must be positive");3706 list_ = list<char>(size, value_left, value_right);3707 }37083709 // Creates generator for strings of size 'size', with random characters in3710 // 'chars'.3711 str(int size, std::set<char> chars) {3712 tgen_ensure(size > 0, "str: size must be positive");3713 list_ = list<char>(size, chars);3714 }37153716 // Creates generator for strings that match the given regex.3717 template <typename... Args> str(const std::string ®ex, Args &&...args) {3718 tgen_ensure(regex.size() > 0, "str: regex must be non-empty");37193720 root_ = detail::parse_regex(3721 detail::regex_format(regex, std::forward<Args>(args)...));3722 }37233724 // Restricts strings for str[idx] = value.3725 str &fix(int idx, char character) {3726 tgen_ensure(!root_, "str: cannot add restriction for regex");3727 list_->fix(idx, character);3728 return *this;3729 }37303731 // Restricts strings for list[S] to be equal, for given subset S of indices.3732 str &equal(std::set<int> indices) {3733 tgen_ensure(!root_, "str: cannot add restriction for regex");3734 list_->equal(indices);3735 return *this;3736 }37373738 // Restricts strings for str[idx_1] = str[idx_2].3739 str &equal(int idx_1, int idx_2) {3740 tgen_ensure(!root_, "str: cannot add restriction for regex");3741 list_->equal(idx_1, idx_2);3742 return *this;3743 }37443745 // Restricts strings for str[left..right] to have all equal values.3746 str &equal_range(int left, int right) {3747 tgen_ensure(!root_, "str: cannot add restriction for regex");3748 list_->equal_range(left, right);3749 return *this;3750 }37513752 // Restricts strings for all equal chars.3753 str &all_equal() {3754 tgen_ensure(!root_, "str: cannot add restriction for regex");3755 list_->all_equal();3756 return *this;3757 }37583759 // Restricts strings for str[left..right] to be a palindrome.3760 str &palindrome(int left, int right) {3761 tgen_ensure(!root_, "str: cannot add restriction for regex");3762 tgen_ensure(0 <= left and left <= right and right < list_->size_,3763 "str: range indices must be valid");3764 for (int i = left; i < right - (i - left); ++i)3765 equal(i, right - (i - left));3766 return *this;3767 }37683769 // Restricts strings for the entire string to be a palindrome.3770 str &palindrome() {3771 tgen_ensure(!root_, "str: cannot add restriction for regex");3772 return palindrome(0, list_->size_ - 1);3773 }37743775 // Restricts strings for str[S] to be different (distinct), for given subset3776 // S of indices.3777 str &different(std::set<int> indices) {3778 tgen_ensure(!root_, "str: cannot add restriction for regex");3779 list_->different(indices);3780 return *this;3781 }37823783 // Restricts strings for str[idx_1] != str[idx_2].3784 str &different(int idx_1, int idx_2) {3785 tgen_ensure(!root_, "str: cannot add restriction for regex");3786 list_->different(idx_1, idx_2);3787 return *this;3788 }37893790 // Restricts lists for list[left..right] to have all different chars.3791 str &different_range(int left, int right) {3792 tgen_ensure(!root_, "str: cannot add restriction for regex");3793 list_->different_range(left, right);3794 return *this;3795 }37963797 // Restricts strings for all chars to be different.3798 str &all_different() {3799 tgen_ensure(!root_, "str: cannot add restriction for regex");3800 list_->all_different();3801 return *this;3802 }38033804 // str value.3805 struct value : gen_value_base<value> {3806 using tgen_is_sequential_tag = detail::is_sequential_tag;38073808 using value_type = char;3809 using std_type = std::string;3810 std::string str_;38113812 value(const std::string &str) : str_(str) {3813 tgen_ensure(!str_.empty(), "str: value: cannot be empty");3814 }38153816 // Fetches size.3817 int size() const { return str_.size(); }38183819 // Fetches position idx.3820 char &operator[](int idx) {3821 tgen_ensure(0 <= idx and idx < size(),3822 "str: value: index out of bounds");3823 return str_[idx];3824 }3825 const char &operator[](int idx) const {3826 tgen_ensure(0 <= idx and idx < size(),3827 "str: value: index out of bounds");3828 return str_[idx];3829 }38303831 // Sorts characters in non-decreasing order.3832 // O(n log n).3833 value &sort() {3834 std::sort(str_.begin(), str_.end());3835 return *this;3836 }38373838 // Reverses string.3839 // O(n).3840 value &reverse() {3841 std::reverse(str_.begin(), str_.end());3842 return *this;3843 }38443845 // Lowercases all characters.3846 // O(n).3847 value &lowercase() {3848 for (char &c : str_)3849 c = std::tolower(c);3850 return *this;3851 }38523853 // Uppercases all characters.3854 // O(n).3855 value &uppercase() {3856 for (char &c : str_)3857 c = std::toupper(c);3858 return *this;3859 }38603861 // Concatenates two values.3862 // Linear.3863 value operator+(const value &rhs) const {3864 return value(str_ + rhs.str_);3865 }38663867 // Shuffles string uniformly.3868 // O(n).3869 value &shuffle() {3870 for (int i = 0; i < size(); ++i)3871 std::swap(str_[i], str_[next(0, size() - 1)]);3872 return *this;3873 }38743875 // Returns a random character uniformly.3876 // O(1).3877 char pick() const { return str_[next<int>(0, size() - 1)]; }38783879 // Returns str_[i] with probability proportional to distribution[i].3880 // O(1).3881 template <typename Dist>3882 char pick_by_distribution(const std::vector<Dist> &distribution) const {3883 tgen_ensure(static_cast<size_t>(size()) == distribution.size(),3884 "value and distribution must have the same size");3885 return str_[next_by_distribution(distribution)];3886 }3887 template <typename Dist>3888 char pick_by_distribution(3889 const std::initializer_list<Dist> &distribution) const {3890 return pick_by_distribution(std::vector<Dist>(distribution));3891 }38923893 // Chooses k characters uniformly, as in a subsequence of size k.3894 // O(n).3895 value choose(int k) const {3896 tgen_ensure(0 < k and k <= size(),3897 "number of elements to choose must be valid");3898 std::string new_str;3899 int need = k;3900 for (int i = 0; need > 0; ++i) {3901 int left = size() - i;3902 if (next(1, left) <= need) {3903 new_str.push_back(str_[i]);3904 need--;3905 }3906 }3907 return value(new_str);3908 }39093910 // Prints to std::ostream.3911 friend std::ostream &operator<<(std::ostream &out, const value &val) {3912 return out << val.str_;3913 }39143915 // Gets a std::string representing the value.3916 std::string to_std() const { return std_type(str_); }3917 };39183919 // Generates str value.3920 // If created from restrictions: O(n log n).3921 // If created from regex: expected linear.3922 value gen() const {3923 if (root_) {3924 // Regex.3925 std::string ret_str;3926 gen_regex(*root_, ret_str);3927 return value(ret_str);3928 } else {3929 // List.3930 std::vector<char> vec = list_->gen().to_std();3931 return value(std::string(vec.begin(), vec.end()));3932 }3933 }3934};39353936/************3937 * *3938 * PAIR *3939 * *3940 ************/39413942namespace detail {39433944// Generates pair first == second.3945// O(1).3946template <typename T> std::pair<T, T> gen_eq(T L1, T R1, T L2, T R2) {3947 T L = std::max(L1, L2);3948 T R = std::min(R1, R2);39493950 tgen_ensure(L <= R, "pair: no valid values to generate");3951 T x = next<T>(L, R);3952 return {x, x};3953}39543955// Returns {R1-L1+1, R2-L2+1}.3956template <typename T>3957std::pair<u128, u128> get_n_and_m(T L1, T R1, T L2, T R2) {3958 u128 n = static_cast<i128>(R1) - L1 + 1;3959 u128 m = static_cast<i128>(R2) - L2 + 1;3960 return {n, m};3961}39623963// Returns first + first+1 + ... + last,3964// num_terms terms. Avoids overflow.3965static u128 pos_arith_sum(u128 first, u128 last, u128 num_terms) {3966 u128 x = first + last, y = num_terms;39673968 // x * y / 2, avoiding overflow.3969 if (x % 2 == 0)3970 x /= 2;3971 else3972 y /= 2;39733974 return x * y;3975}39763977// Generates pair first != second.3978// O(1) expected.3979template <typename T> std::pair<T, T> gen_neq(T L1, T R1, T L2, T R2) {3980 auto [n, m] = get_n_and_m(L1, R1, L2, R2);39813982 T L_intersect = std::max(L1, L2);3983 T R_intersect = std::min(R1, R2);3984 u128 inter = static_cast<i128>(R_intersect) - L_intersect + 1;39853986 u128 total = n * m - inter;3987 tgen_ensure(total > 0, "pair: no valid values to generate");39883989 // Runs O(1) expected times in the worst case.3990 T a, b;3991 do {3992 a = next<T>(L1, R1);3993 b = next<T>(L2, R2);3994 } while (a == b);39953996 return {a, b};3997}39983999// For lt, splits 'second' into two regions:4000// 1) second <= R1 -> number of 'first' is (second - L1)4001// 2) second > R1 -> number of 'first' is (R1 - L1 + 1)4002// Returns {count_region1, count_region2}.4003// O(1).4004template <typename T>4005std::pair<u128, u128> count_lt_regions(T L1, T R1, T L2, T R2) {4006 auto [n, m] = get_n_and_m(L1, R1, L2, R2);40074008 // 'second' must be >= L1 + 1.4009 i128 L_second = std::max<i128>(L2, static_cast<i128>(L1) + 1);4010 i128 R_second = R2;40114012 // Split point for 'second'.4013 i128 split = std::min<i128>(R_second, R1);40144015 // Region 1: b in [L_second, split].4016 u128 len1 = std::max<i128>(0, split - L_second + 1);40174018 u128 count_region1 = 0;4019 if (len1 > 0) {4020 // For b in [L_second, split], there are (b - L1) ways.4021 i128 first = L_second - L1;4022 i128 last = split - L1;40234024 // Arithmetic series first + (first + 1) + ... + last, len1 terms.4025 count_region1 = pos_arith_sum(first, last, len1);4026 }40274028 // Region 2: b > R1.4029 // For b in [R1+1, R_second], there are 'n' ways.4030 i128 L_second_region2 = std::max(L_second, static_cast<i128>(R1) + 1);40314032 u128 len2 = std::max<i128>(0, R_second - L_second_region2 + 1);4033 u128 count_region2 = len2 * n;40344035 return {count_region1, count_region2};4036}40374038// Generates pair first < second.4039// O(log(R1 - L1 + 1) + log(R2 - L2 + 1)).4040template <typename T> std::pair<T, T> gen_lt(T L1, T R1, T L2, T R2) {4041 auto [n, m] = get_n_and_m(L1, R1, L2, R2);40424043 // 'second' needs to be at least L1 + 1 to have a valid value for4044 // 'first'.4045 i128 L_second = std::max<i128>(L2, static_cast<i128>(L1) + 1);4046 i128 R_second = R2;40474048 // Splits 'second' into two regions:4049 // 1) b <= R1 -> number of 'first' is (b - L1);4050 // 2) b > R1 -> number of 'first' is (R1 - L1 + 1).4051 i128 split = std::min<i128>(R_second, R1);40524053 auto [count_region1, count_region2] = count_lt_regions(L1, R1, L2, R2);4054 u128 total = count_region1 + count_region2;4055 tgen_ensure(total > 0, "pair: no valid values to generate");40564057 u128 k = detail::next128(total);4058 if (k < count_region1) {4059 // Region 1: invert arithmetic series.40604061 // For b in [L_second, split].4062 u128 len1 = std::max<i128>(0, split - L_second + 1);40634064 // We consider b in [L_second, L_second + d].4065 // Each b contributes (b - L1) = base + (b - L_second).4066 // So we sum: base + (base+1) + ... + (base+d)4067 // d in [0, len1).40684069 i128 base = L_second - L1;4070 i128 lo = 0, hi = static_cast<i128>(len1) - 1;40714072 while (lo < hi) {4073 i128 mid = lo + (hi - lo) / 2;40744075 if (pos_arith_sum(base, base + mid, mid + 1) <= k)4076 lo = mid + 1;4077 else4078 hi = mid;4079 }4080 i128 d = lo;40814082 // Subtracts prefix sum with d-1 terms from k.4083 if (d > 0)4084 k -= pos_arith_sum(base, base + d - 1, d);40854086 return {L1 + static_cast<T>(k), L_second + d};4087 } else {4088 // Region 2: uniform block of size n.4089 k -= count_region1;40904091 // For b in [R1+1, R_second], there are 'n' ways.4092 i128 L_second_region2 = std::max(L_second, static_cast<i128>(R1) + 1);40934094 return {L1 + static_cast<T>(k % n),4095 L_second_region2 + static_cast<T>(k / n)};4096 }4097}40984099// Generates pair first > second.4100// O(log(R1 - L1 + 1) + log(R2 - L2 + 1)).4101template <typename T> std::pair<T, T> gen_gt(T L1, T R1, T L2, T R2) {4102 auto [first, second] = gen_lt(L2, R2, L1, R1);4103 return {second, first};4104}41054106// Generates pair first <= second.4107// O(log(R1 - L1 + 1) + log(R2 - L2 + 1)).4108template <typename T> std::pair<T, T> gen_leq(T L1, T R1, T L2, T R2) {4109 // Counts how many pairs are there with first = second.4110 i128 L_intersect = std::max(L1, L2);4111 i128 R_intersect = std::min(R1, R2);4112 u128 eq_count = std::max<i128>(0, R_intersect - L_intersect + 1);41134114 // Counts how many pairs are there with first < second.4115 auto [lt_region1, lt_region2] = count_lt_regions(L1, R1, L2, R2);4116 u128 lt_count = lt_region1 + lt_region2;41174118 u128 total = eq_count + lt_count;4119 tgen_ensure(total > 0, "pair: no valid values to generate");41204121 if (detail::next128(total) < eq_count)4122 return gen_eq(L1, R1, L2, R2);4123 return gen_lt(L1, R1, L2, R2);4124}41254126// Generates pair first >= second.4127// O(log(R1 - L1 + 1) + log(R2 - L2 + 1)).4128template <typename T> std::pair<T, T> gen_geq(T L1, T R1, T L2, T R2) {4129 auto [first, second] = gen_leq(L2, R2, L1, R1);4130 return {second, first};4131}41324133}; // namespace detail41344135/*4136 * Pair generator.4137 *4138 * Pairs of integral types.4139 */41404141template <typename T> struct pair : gen_base<pair<T>> {4142 std::pair<T, T> first_, second_; // Range of first and second values.4143 // Type of restriction.4144 enum class restriction_type { eq, neq, lt, gt, leq, geq, unspecified };4145 restriction_type type_ = restriction_type::unspecified;41464147 // Creates a pair with random values in [first_l, first_r] and [second_l,4148 // second_r].4149 pair(T first_left, T first_right, T second_left, T second_right)4150 : first_(first_left, first_right), second_(second_left, second_right) {4151 tgen_ensure(first_left <= first_right,4152 "pair: first range must be valid");4153 tgen_ensure(second_left <= second_right,4154 "pair: second range must be valid");4155 }41564157 // Creates a pair with random values in [both_l, both_r].4158 pair(T both_left, T both_right)4159 : pair(both_left, both_right, both_left, both_right) {}41604161 // Restricts pair for first = second.4162 pair &eq() {4163 type_ = restriction_type::eq;4164 return *this;4165 }41664167 // Restricts pair for first != second.4168 pair &neq() {4169 type_ = restriction_type::neq;4170 return *this;4171 }41724173 // Restricts pair for first < second.4174 pair <() {4175 type_ = restriction_type::lt;4176 return *this;4177 }41784179 // Restricts pair for first > second.4180 pair >() {4181 type_ = restriction_type::gt;4182 return *this;4183 }41844185 // Restricts pair for first <= second.4186 pair &leq() {4187 type_ = restriction_type::leq;4188 return *this;4189 }41904191 // Restricts pair for first >= second.4192 pair &geq() {4193 type_ = restriction_type::geq;4194 return *this;4195 }41964197 // Pair value.4198 struct value : gen_value_base<value> {4199 using value_type = T;4200 using std_type = std::pair<T, T>;42014202 std::pair<T, T> pair_;4203 char sep_;42044205 value(const std::pair<T, T> &pair) : pair_(pair), sep_(' ') {}4206 value(const T &first, const T &second)4207 : pair_(first, second), sep_(' ') {}42084209 T first() const { return pair_.first; }4210 T second() const { return pair_.second; }42114212 // Sets the separator for the pair, for printing.4213 value &separator(char sep) {4214 sep_ = sep;4215 return *this;4216 }42174218 // Prints to std::ostream, separated by sep_.4219 friend std::ostream &operator<<(std::ostream &out, const value &val) {4220 return out << val.pair_.first << val.sep_ << val.pair_.second;4221 }42224223 // Gets a std::pair representing the value.4224 auto to_std() const {4225 if constexpr (!detail::is_generator_value<T>::value) {4226 return pair_;4227 } else {4228 std::pair<typename T::std_type, typename T::std_type> pair(4229 pair_.first.to_std(), pair_.second.to_std());4230 return pair;4231 }4232 }4233 };42344235 // Generates a random pair.4236 // O(log(R1 - L1 + 1) + log(R2 - L2 + 1)).4237 value gen() const {4238 T L1 = first_.first, R1 = first_.second;4239 T L2 = second_.first, R2 = second_.second;42404241 switch (type_) {4242 case restriction_type::unspecified:4243 return {next<T>(L1, R1), next<T>(L2, R2)};4244 case restriction_type::eq:4245 return detail::gen_eq<T>(L1, R1, L2, R2);4246 case restriction_type::neq:4247 return detail::gen_neq<T>(L1, R1, L2, R2);4248 case restriction_type::lt:4249 return detail::gen_lt<T>(L1, R1, L2, R2);4250 case restriction_type::gt:4251 return detail::gen_gt<T>(L1, R1, L2, R2);4252 case restriction_type::leq:4253 return detail::gen_leq<T>(L1, R1, L2, R2);4254 case restriction_type::geq:4255 return detail::gen_geq<T>(L1, R1, L2, R2);4256 }4257 throw detail::error("pair: unknown restriction type");4258 }4259};42604261/************4262 * *4263 * TREE *4264 * *4265 ************/42664267namespace detail {42684269// Generates edges from Prufer sequence.4270// O(n).4271inline std::vector<std::pair<int, int>> edges_from_prufer(std::vector<int> p) {4272 int n = p.size() + 2;42734274 // Degrees.4275 std::vector<int> d(n, 1);4276 for (int i : p)4277 d[i]++;42784279 // Adds last vertex.4280 p.push_back(n - 1);42814282 // Finds first vertex with degree 1.4283 int idx, u;4284 idx = u = find(d.begin(), d.end(), 1) - d.begin();42854286 // Generates edges.4287 std::vector<std::pair<int, int>> edges;4288 for (int v : p) {4289 edges.emplace_back(u, v);4290 if (--d[v] == 1 and v < idx)4291 u = v;4292 else4293 idx = u = find(d.begin() + idx + 1, d.end(), 1) - d.begin();4294 }4295 return edges;4296}42974298// Disjoint set union (union-find) for connectivity queries.4299struct dsu {4300 std::vector<int> parent_;4301 std::vector<unsigned char> rank_;43024303 // Creates a dsu with `n` elements, indexed from 0 to n-1.4304 // Initially every element is in its own set.4305 // O(n).4306 dsu(int n) : parent_(n), rank_(n, 0) {4307 for (int i = 0; i < n; ++i)4308 parent_[i] = i;4309 }43104311 // Adds new elements to the dsu, each in their own new set.4312 // O(k) amortized.4313 void add_elements(int k) {4314 for (int i = 0; i < k; ++i) {4315 int new_id = parent_.size();4316 parent_.push_back(new_id);4317 rank_.push_back(0);4318 }4319 }43204321 // Finds representative of set containing i.4322 // O(alpha(n)) amortized, O(log n) worst case.4323 int find(int i) {4324 return parent_[i] == i ? i : parent_[i] = find(parent_[i]);4325 }43264327 // Merges components of `a` and `b`. Returns if the sets were united, and4328 // false if a and b were in the same set.4329 // O(alpha(n)) amortized, O(log n) worst case.4330 bool unite(int a, int b) {4331 a = find(a);4332 b = find(b);4333 if (a == b)4334 return false;4335 if (rank_[a] > rank_[b])4336 std::swap(a, b);4337 parent_[a] = b;4338 if (rank_[a] == rank_[b])4339 ++rank_[b];4340 return true;4341 }4342};43434344} // namespace detail43454346// Forward declaration of wgraph.4347template <typename VWeight, typename EWeight> struct wgraph;43484349/*4350 * Tree generator.4351 *4352 * Unrooted trees with `n` vertices, indexed from 0 to n-1.4353 * These are unrooted undirected labeled trees, that is, isomorphism is not4354 * taken into account. VWeight is the type of vertex weights, and EWeight is4355 * the type of edge weights. Generator does not generate weights. The weights4356 * are to be set in the wtree::value.4357 */43584359template <typename VWeight, typename EWeight>4360struct wtree : gen_base<wtree<VWeight, EWeight>> {4361 int n_; // Number of vertices.4362 std::set<std::pair<int, int>> edges_; // Edges that were set.43634364 // Creates tree generator with `n` vertices.4365 // O(1).4366 wtree(int n) : n_(n) {4367 tgen_ensure(n > 0, "wtree: number of vertices must be positive");4368 }43694370 // Adds edge between u and v (this edge must be generated).4371 // O(log n).4372 wtree &add_edge(int u, int v) {4373 tgen_ensure(0 <= std::min(u, v) and std::max(u, v) < n_,4374 "wtree: vertices must be indexed in [0, n)");4375 tgen_ensure(u != v, "wtree: cannot add self loop to tree");43764377 if (u > v)4378 std::swap(u, v);4379 edges_.emplace(u, v);4380 return *this;4381 }43824383 // Tree value.4384 //4385 // Edges are stored in both directions in adjacency list, but only u < v in4386 // edge list.4387 struct value : gen_value_base<value> {4388 using std_type = std::pair<int, std::vector<std::set<int>>>;43894390 int n_; // Number of vertices.4391 std::vector<std::set<int>> adj_; // Adjacency list.4392 std::vector<std::pair<int, int>> edges_; // Edge list.4393 bool add_1_; // If should add 1 for printing vertex ids.4394 bool print_n_; // If should print n.4395 std::optional<int> print_parents_; // If should print in parent style4396 // (stores the root).4397 std::optional<std::vector<VWeight>> vertex_weights_; // Vertex weights.4398 std::optional<std::vector<EWeight>>4399 edge_weights_; // Edge weights (in same order as edges_).4400 detail::dsu dsu_; // Connectivity of current edges (for cycle checks).44014402 // Creates value from adjacency list.4403 // O(n).4404 value(const std::vector<std::set<int>> &adj)4405 : n_(static_cast<int>(adj.size())), adj_(adj), add_1_(false),4406 print_n_(false), dsu_(n_) {4407 for (int u = 0; u < n_; ++u)4408 for (auto v : adj[u]) {4409 tgen_ensure(4410 0 <= v and v < n_,4411 "wtree: value: vertices must be indexed in [0, n)");4412 // Symmetric adjacency: count each undirected edge once.4413 if (u < v) {4414 edges_.emplace_back(u, v);4415 tgen_ensure(4416 dsu_.unite(u, v),4417 "wtree: value: initial graph must form a tree");4418 }4419 }4420 }44214422 // Creates value from `n` and edge list.4423 // O(n).4424 value(int n, const std::vector<std::pair<int, int>> &edges)4425 : n_(n), adj_(n), add_1_(false), print_n_(false), dsu_(n) {4426 edges_.reserve(edges.size());4427 for (auto [u, v] : edges) {4428 tgen_ensure(0 <= std::min(u, v) and std::max(u, v) < n,4429 "wtree: value: vertices must be indexed in [0, n)");4430 tgen_ensure(dsu_.unite(u, v),4431 "wtree: value: initial graph must form a tree");4432 if (u > v)4433 std::swap(u, v);4434 edges_.emplace_back(u, v);4435 adj_[u].insert(v);4436 adj_[v].insert(u);4437 }4438 }4439 value(int n, const std::set<std::pair<int, int>> &edges)4440 : value(n, std::vector<std::pair<int, int>>(edges.begin(),4441 edges.end())) {}4442 value(int n, const std::initializer_list<std::pair<int, int>> &edges)4443 : value(n, std::vector<std::pair<int, int>>(edges)) {}44444445 // Creates tree from graph via Kruskal-like random spanning tree.4446 // Implemented after wgraph definition.4447 // O(n + m alpha(n)).4448 value(const typename wgraph<VWeight, EWeight>::value &g);44494450 // Weight type conversion.4451 // O(n).4452 template <typename NewVWeight, typename NewEWeight>4453 typename wtree<NewVWeight, NewEWeight>::value4454 convert_weight_types() const {4455 tgen_ensure(!vertex_weights_.has_value() and4456 !edge_weights_.has_value(),4457 "wtree: value: cannot convert weight type after "4458 "assigning weights");44594460 typename wtree<NewVWeight, NewEWeight>::value new_tree(adj_);4461 new_tree.add_1_ = add_1_;4462 new_tree.print_n_ = print_n_;4463 new_tree.print_parents_ = print_parents_;4464 return new_tree;4465 }44664467 // Fetches number of vertices.4468 int n() const { return n_; }44694470 // Fetches a const ref. to adjacency list.4471 const std::vector<std::set<int>> &adj() const { return adj_; }44724473 // Fetches a const ref. to edge list.4474 const std::vector<std::pair<int, int>> &edges() const { return edges_; }44754476 // Fetches a const ref. to vertex weights.4477 const std::optional<std::vector<VWeight>> &vertex_weights() const {4478 return vertex_weights_;4479 }44804481 // Fetches a const ref. to edge weights.4482 const std::optional<std::vector<EWeight>> &edge_weights() const {4483 return edge_weights_;4484 }44854486 // Sets vertex weights.4487 // O(n).4488 template <typename NewVWeight = VWeight>4489 typename wtree<NewVWeight, EWeight>::value set_vertex_weights(4490 const std::vector<NewVWeight> &vertex_weights) const {4491 tgen_ensure(static_cast<int>(vertex_weights.size()) == n(),4492 "wtree: value: must give `n` vertex weights");44934494 auto new_tree = convert_weight_types<NewVWeight, EWeight>();4495 new_tree.vertex_weights_ = vertex_weights;4496 return new_tree;4497 }44984499 // Sets edge weights.4500 // O(n).4501 template <typename NewEWeight = EWeight>4502 typename wtree<VWeight, NewEWeight>::value4503 set_edge_weights(const std::vector<NewEWeight> &edge_weights) const {4504 tgen_ensure(4505 edge_weights.size() == edges().size(),4506 "wtree: value: must give `edges().size()` edge weights");45074508 auto new_tree = convert_weight_types<VWeight, NewEWeight>();4509 new_tree.edge_weights_ = edge_weights;4510 return new_tree;4511 }45124513 // Enables edge-weighted mode before adding weighted edges4514 // incrementally. The tree must have no edges yet. O(1).4515 value &edge_weighted() {4516 tgen_ensure(edges().size() == 0,4517 "wtree: value: edge_weighted requires a tree with no "4518 "edges");4519 tgen_ensure(!edge_weights_.has_value(),4520 "wtree: value: tree is already edge-weighted");45214522 edge_weights_ = std::vector<EWeight>();4523 return *this;4524 }45254526 // Adds 1 to vertex ids, for printing.4527 // O(1).4528 value &add_1() {4529 add_1_ = true;4530 return *this;4531 }45324533 // Prints `n` on a new line before printing the tree.4534 // O(1).4535 value &print_n() {4536 print_n_ = true;4537 return *this;4538 }45394540 // Prints the tree in parent style.4541 // If root = -1, the root is considered to be 0, and its parent is not4542 // printed. Otherwise, prints the parent of the root as -1. If root = n,4543 // randomizes the root. O(1).4544 value &print_parents(int root = -1) {4545 tgen_ensure(root == -1 or (0 <= root and root < n()) or root == n(),4546 "wtree: value: root must be -1, `n`, or in [0, n)");4547 print_parents_ = root;4548 return *this;4549 }45504551 // Shuffles the tree's vertex labels (except those in `indices`,4552 // which keep their current label) and edge order. The change is4553 // applied eagerly to the underlying adjacency list, edge list,4554 // vertex weights and edge weights.4555 // O(n).4556 value &shuffle_except(std::set<int> indices) {4557 // Builds the relabeling: for each vertex `i`, `new_label[i]` is4558 // its new id. Vertices in `indices` keep their label; the others4559 // are permuted among themselves.4560 std::vector<int> new_label(n());4561 std::vector<int> shuffled;4562 for (int i = 0; i < n(); ++i) {4563 if (indices.count(i))4564 new_label[i] = i;4565 else4566 shuffled.push_back(i);4567 }4568 std::vector<int> targets = shuffled;4569 tgen::shuffle(targets.begin(), targets.end());4570 for (size_t k = 0; k < shuffled.size(); ++k)4571 new_label[shuffled[k]] = targets[k];45724573 // Rewrites adjacency list with new labels.4574 std::vector<std::set<int>> new_adj(n());4575 for (int u = 0; u < n(); ++u)4576 for (int v : adj_[u])4577 new_adj[new_label[u]].insert(new_label[v]);4578 adj_ = std::move(new_adj);45794580 // Rewrites edges with new labels (canonical undirected order).4581 for (auto &[u, v] : edges_) {4582 u = new_label[u];4583 v = new_label[v];4584 if (u > v)4585 std::swap(u, v);4586 }45874588 // Permutes vertex weights to match the new labels.4589 if (vertex_weights_.has_value()) {4590 std::vector<VWeight> new_vw(n());4591 for (int i = 0; i < n(); ++i)4592 new_vw[new_label[i]] = (*vertex_weights_)[i];4593 vertex_weights_ = std::move(new_vw);4594 }45954596 // Rebuilds the dsu so future `add_edge` calls see the new labels.4597 dsu_ = detail::dsu(n());4598 for (auto [u, v] : edges_)4599 dsu_.unite(u, v);46004601 // Shuffles edge order, keeping edge weights aligned.46024603 std::vector<int> perm(edges_.size());4604 std::iota(perm.begin(), perm.end(), 0);4605 tgen::shuffle(perm.begin(), perm.end());46064607 std::vector<std::pair<int, int>> new_edges;4608 std::optional<std::vector<EWeight>> new_ew;4609 if (edge_weights_.has_value())4610 new_ew = std::vector<EWeight>();4611 for (int i : perm) {4612 new_edges.push_back(edges_[i]);4613 if (new_ew.has_value())4614 new_ew->push_back((*edge_weights_)[i]);4615 }4616 edges_ = new_edges;4617 if (new_ew.has_value())4618 edge_weights_ = new_ew;46194620 return *this;4621 }46224623 // Shuffles the tree's vertices and edge order.4624 // O(n).4625 value &shuffle() { return shuffle_except({}); }46264627 // Adds edge (u, v).4628 // O(log n) amortized.4629 value &add_edge(int u, int v, std::optional<EWeight> w = std::nullopt) {4630 tgen_ensure(0 <= std::min(u, v) and std::max(u, v) < n(),4631 "wtree: value: vertex ids must be valid");46324633 if (u > v)4634 std::swap(u, v);46354636 if (adj_[u].count(v))4637 return *this;46384639 adj_[u].insert(v);4640 adj_[v].insert(u);4641 edges_.emplace_back(u, v);4642 tgen_ensure(dsu_.unite(u, v),4643 "wtree: value: added edge must not create a cycle");46444645 if (w.has_value()) {4646 tgen_ensure(edge_weights().has_value(),4647 "wtree: value: cannot add weighted edge to "4648 "edge-unweighted tree");46494650 edge_weights_->push_back(*w);4651 } else4652 tgen_ensure(!edge_weights().has_value(),4653 "wtree: value: cannot add unweighted edge to "4654 "edge-weighted tree");46554656 return *this;4657 }46584659 // Links tree with another `rhs`, adding the edge between u (in left4660 // tree) and v (in right tree). Ids for added vertices are updated4661 // accordingly.4662 // O(rhs.n + rhs.m * log n) amortized.4663 value &link(const value &rhs, int new_u, int new_v,4664 std::optional<EWeight> new_w = std::nullopt) {4665 tgen_ensure(0 <= new_u and new_u < n() and 0 <= new_v and4666 new_v < rhs.n(),4667 "wtree: value: vertex ids must be valid");46684669 // Edges from right-hand side.4670 int shift = n();4671 add_vertices(rhs.n(), rhs.vertex_weights());4672 for (int i = 0; i < static_cast<int>(rhs.edges().size()); ++i) {4673 auto [u, v] = rhs.edges()[i];4674 add_edge(shift + u, shift + v,4675 rhs.edge_weights().has_value()4676 ? std::optional<EWeight>((*rhs.edge_weights())[i])4677 : std::nullopt);4678 }46794680 // New edge.4681 add_edge(new_u, shift + new_v, new_w);46824683 return *this;4684 }46854686 // Glues the tree with another `rhs` such that index_pairs[i].first is4687 // considered to be the same as index_pairs[i].second. Ids for added4688 // vertices are updated accordingly.4689 // O(rhs.n + rhs.m * log n) amortized.4690 value &glue(const value &rhs,4691 std::set<std::pair<int, int>> index_pairs) {4692 // Checks validity of indices.4693 std::set<int> idx_left, idx_right;4694 std::vector<int> right_id_to_left(rhs.n(), -1);4695 for (auto [l, r] : index_pairs) {4696 tgen_ensure(4697 0 <= l and l < n() and 0 <= r and r < rhs.n(),4698 "wtree: value: vertex indices to glue must be valid");4699 tgen_ensure(idx_left.count(l) == 0 and idx_right.count(r) == 0,4700 "wtree: value: must not have repeated indices "4701 "on the same side to glue");47024703 idx_left.insert(l);4704 idx_right.insert(r);4705 right_id_to_left[r] = l;4706 }47074708 // Computes new ids of right vertices.4709 std::vector<int> new_right_id(rhs.n(), -1);4710 int intersection_lt = 0;4711 std::optional<std::vector<VWeight>> rhs_vertex_weights;4712 for (int i = 0; i < rhs.n(); ++i) {4713 if (right_id_to_left[i] != -1) {4714 // Is in intersection.4715 ++intersection_lt;4716 new_right_id[i] = right_id_to_left[i];4717 } else {4718 // New id.4719 new_right_id[i] = n() + i - intersection_lt;4720 if (rhs.vertex_weights().has_value()) {4721 if (!rhs_vertex_weights.has_value())4722 rhs_vertex_weights = std::vector<VWeight>();4723 rhs_vertex_weights->push_back(4724 (*rhs.vertex_weights())[i]);4725 }4726 }4727 }47284729 // Adds new vertices and edges.4730 add_vertices(rhs.n() - intersection_lt, rhs_vertex_weights);4731 for (int i = 0; i < static_cast<int>(rhs.edges().size()); ++i) {4732 auto [u, v] = rhs.edges()[i];4733 add_edge(new_right_id[u], new_right_id[v],4734 rhs.edge_weights().has_value()4735 ? std::optional<EWeight>((*rhs.edge_weights())[i])4736 : std::nullopt);4737 }47384739 return *this;4740 }4741 value &glue(const value &rhs,4742 std::initializer_list<std::pair<int, int>> il) {4743 return glue(rhs, std::set<std::pair<int, int>>(il));4744 }47454746 // Glues the tree with another `rhs` at `indices`. That is, idx in4747 // `indices` are considered to be the same vertex. Ids for added4748 // vertices are updated accordingly.4749 // O(rhs.n).4750 value &glue(const value &rhs, std::set<int> indices) {4751 std::set<std::pair<int, int>> index_pairs;4752 for (auto i : indices)4753 index_pairs.emplace(i, i);4754 return glue(rhs, index_pairs);4755 }4756 value &glue(const value &rhs, const std::initializer_list<int> &il) {4757 return glue(rhs, std::set<int>(il));4758 }47594760 // Prints to std::ostream.4761 // O(n).4762 friend std::ostream &operator<<(std::ostream &out, const value &val) {4763 if (val.print_n_)4764 out << val.n() << '\n';47654766 // Prints vertex weights.4767 if (val.vertex_weights()) {4768 for (int i = 0; i < val.n(); ++i) {4769 if (i > 0)4770 out << " ";4771 out << (*val.vertex_weights())[i];4772 }4773 out << '\n';4774 }47754776 tgen_ensure(static_cast<int>(val.edges().size()) == val.n() - 1,4777 "wtree: value: invalid tree to print (number of edges "4778 "must be `n` - 1)");47794780 // Prints in parent style.4781 if (val.print_parents_.has_value()) {4782 tgen_ensure(!val.edge_weights().has_value(),4783 "wtree: value: cannot print parent style if edges "4784 "are weighted");47854786 int root = *val.print_parents_;4787 bool skip_parent_0 = root == -1;4788 if (root == -1)4789 root = 0;4790 if (root == val.n())4791 root = next(0, val.n() - 1);47924793 std::vector<int> parent(val.n(), -1);47944795 std::queue<int> q;4796 std::vector<int> vis(val.n(), false);4797 q.push(root);4798 vis[root] = true;47994800 while (q.size()) {4801 int u = q.front();4802 q.pop();4803 for (int v : val.adj()[u])4804 if (!vis[v]) {4805 vis[v] = true;4806 q.push(v);4807 parent[v] = u;4808 }4809 }48104811 if (skip_parent_0) {4812 for (int i = 1; i < val.n(); ++i) {4813 tgen_ensure(4814 parent[i] < i,4815 "wtree: value: parent of i must be less than i for "4816 "printing in parent style if root is -1");48174818 if (i > 1)4819 out << " ";4820 out << parent[i] + val.add_1_;4821 }4822 } else {4823 for (int i = 0; i < val.n(); ++i) {4824 if (i > 0)4825 out << " ";4826 out << (parent[i] == -1 ? -1 : parent[i]) + val.add_1_;4827 }4828 }48294830 out << '\n';4831 return out;4832 }48334834 // Prints edges.4835 for (int i = 0; i < static_cast<int>(val.edges().size()); ++i) {4836 auto [u, v] = val.edges()[i];4837 out << (u + val.add_1_) << " " << (v + val.add_1_);48384839 // Edge weight.4840 if (val.edge_weights().has_value())4841 out << " " << (*val.edge_weights())[i];48424843 out << '\n';4844 }48454846 return out;4847 }48484849 // Gets a std::pair<n, adj> representing the value.4850 std::pair<int, std::vector<std::set<int>>> to_std() const {4851 return std_type(n_, adj_);4852 }48534854 private:4855 // Adds `k` vertices to the tree (labeled n, n+1, ...n+k-1). Updates4856 // `n` accordingly. This makes the tree invalid (not a tree anymore).4857 // O(k) amortized.4858 value &add_vertices(int k, std::optional<std::vector<VWeight>>4859 new_vertex_weights = std::nullopt) {4860 n_ += k;4861 adj_.resize(n());4862 if (new_vertex_weights.has_value()) {4863 tgen_ensure(vertex_weights().has_value(),4864 "wtree: value: cannot add weighted vertices to "4865 "vertex-unweighted tree");4866 tgen_ensure(4867 static_cast<int>(new_vertex_weights->size()) == k,4868 "wtree: value: number of vertex weights must be equal "4869 "to number of added vertices");48704871 vertex_weights_->insert(vertex_weights_->end(),4872 new_vertex_weights->begin(),4873 new_vertex_weights->end());4874 } else4875 tgen_ensure(!vertex_weights().has_value(),4876 "wtree: value: cannot add unweighted vertices to "4877 "vertex-weighted tree");48784879 dsu_.add_elements(k);48804881 return *this;4882 }4883 };48844885 // Generates tree value.4886 // O(n).4887 value gen() const {4888 // Constructs adjacency list.4889 std::vector<std::vector<int>> adj(n_);4890 for (auto [u, v] : edges_) {4891 adj[u].push_back(v);4892 adj[v].push_back(u);4893 }48944895 std::vector<int> comp_size;4896 std::vector<std::vector<int>> component_ids;4897 std::vector<bool> vis(n_, false);4898 std::queue<int> q;48994900 for (int i = 0; i < n_; ++i) {4901 if (vis[i])4902 continue;49034904 vis[i] = true;4905 q.push(i);4906 comp_size.push_back(0);4907 component_ids.emplace_back();4908 while (q.size()) {4909 int u = q.front();4910 q.pop();4911 ++comp_size.back();4912 component_ids.back().push_back(u);4913 for (int v : adj[u]) {4914 if (!vis[v]) {4915 vis[v] = true;4916 q.push(v);4917 }4918 }4919 }4920 }49214922 // Creates edges connecting the connected components by treating them as4923 // vertices.4924 std::vector<std::pair<int, int>> new_edges(edges_.begin(),4925 edges_.end());4926 if (comp_size.size() > 1) {4927 std::vector<int> prufer_values =4928 many_by_distribution(comp_size.size() - 2, comp_size);4929 for (auto [u, v] : detail::edges_from_prufer(prufer_values))4930 new_edges.emplace_back(pick(component_ids[u]),4931 pick(component_ids[v]));4932 }49334934 return value(n_, new_edges);4935 }49364937 // Generates a (not uniformly) random skewed tree.4938 // Vertex 0 is the root. For each i in 1 .. n-1, parent(i) is4939 // wnext(i, elongation), i.e. a value in [0, i) with skew controlled by4940 // elongation (see wnext).4941 // If elongation is small enough, generates a star (center 0).4942 // If elongation is large enough, generates a path (endpoints 0 and n-1).4943 // O(n).4944 static value gen_skewed(int n, int elongation) {4945 std::vector<std::pair<int, int>> edges;4946 for (int i = 1; i < n; ++i)4947 edges.emplace_back(i, wnext<int>(i, elongation));4948 return value(n, edges);4949 }49504951 // Kruskal-like random tree: random vertex pairs until connected.4952 // Not uniformly random.4953 // O(n log(n) alpha(n)) expected.4954 static value gen_kruskal(int n) {4955 tgen_ensure(n > 0, "wtree: gen_kruskal: n must be positive");4956 if (n == 1)4957 return value(1, {});49584959 detail::dsu components(n);4960 std::vector<std::pair<int, int>> edges;4961 edges.reserve(n - 1);4962 while (edges.size() < size_t(n - 1)) {4963 int u = next(0, n - 1);4964 int v = next(0, n - 1);4965 if (u == v)4966 continue;4967 if (components.unite(u, v))4968 edges.emplace_back(u, v);4969 }4970 return value(n, edges);4971 }4972};49734974/*4975 * Other types of weighted-ness.4976 */49774978// Vertex weighted tree.4979template <typename VWeight> using vtree = wtree<VWeight, int>;49804981// Edge weighted tree.4982template <typename EWeight> using etree = wtree<int, EWeight>;49834984// Unweighted tree.4985using tree = wtree<int, int>;49864987/*************4988 * *4989 * GRAPH *4990 * *4991 *************/49924993namespace detail {49944995// Canonical undirected edge key for duplicate detection; stores (min(u, v),4996// max(u, v)). O(1).4997inline uint64_t undirected_edge_key(int u, int v) {4998 if (u > v)4999 std::swap(u, v);5000 return (static_cast<uint64_t>(u) << 32) |5001 static_cast<uint64_t>(static_cast<uint32_t>(v));5002}50035004// Directed edge key for duplicate detection; stores (u, v).5005// O(1).5006inline uint64_t directed_edge_key(int u, int v) {5007 return (static_cast<uint64_t>(u) << 32) |5008 static_cast<uint64_t>(static_cast<uint32_t>(v));5009}50105011// Maximum number of edges in a simple graph on n vertices.5012// O(1).5013inline long long max_graph_edges(int n, bool directed, bool self_loops) {5014 if (n <= 0)5015 return 0;5016 if (directed)5017 return self_loops ? static_cast<long long>(n) * n5018 : static_cast<long long>(n) * (n - 1);5019 return self_loops ? static_cast<long long>(n) * (n + 1) / 25020 : static_cast<long long>(n) * (n - 1) / 2;5021}50225023// Uniform random edge for rejection sampling.5024// O(1) expected.5025inline std::pair<int, int> get_random_graph_edge(int n, bool directed,5026 bool self_loops) {5027 if (directed) {5028 if (self_loops)5029 return {next<int>(0, n - 1), next<int>(0, n - 1)};5030 int u = next<int>(0, n - 1);5031 int v = next<int>(0, n - 1);5032 while (u == v)5033 v = next<int>(0, n - 1);5034 return {u, v};5035 }5036 if (self_loops) {5037 int u = next<int>(0, n - 1);5038 int v = next<int>(0, n - 1);5039 if (u > v)5040 std::swap(u, v);5041 return {u, v};5042 }5043 int u = next<int>(0, n - 1);5044 int v = next<int>(0, n - 1);5045 while (u == v)5046 v = next<int>(0, n - 1);5047 if (u > v)5048 std::swap(u, v);5049 return {u, v};5050}50515052// Decodes a linear edge index to (u, v) for an undirected simple graph,5053// with u < v.5054// O(log n).5055inline std::pair<int, int> decode_undirected_simple_edge(int n, long long idx) {5056 auto base = [&](int u) -> long long {5057 return static_cast<long long>(u) * (n - 1) -5058 static_cast<long long>(u) * (u - 1) / 2;5059 };5060 int lo = 0, hi = n - 2;5061 while (lo < hi) {5062 int mid = (lo + hi + 1) / 2;5063 if (base(mid) <= idx)5064 lo = mid;5065 else5066 hi = mid - 1;5067 }5068 return {lo, lo + 1 + int(idx - base(lo))};5069}50705071// Decodes a linear edge index to (u, v) for an undirected graph with loops,5072// with u <= v.5073// O(log n).5074inline std::pair<int, int> decode_undirected_loops_edge(int n, long long idx) {5075 auto base = [&](int u) -> long long {5076 return static_cast<long long>(u) * n -5077 static_cast<long long>(u) * (u - 1) / 2;5078 };5079 int lo = 0, hi = n - 1;5080 while (lo < hi) {5081 int mid = (lo + hi + 1) / 2;5082 if (base(mid) <= idx)5083 lo = mid;5084 else5085 hi = mid - 1;5086 }5087 return {lo, lo + int(idx - base(lo))};5088}50895090// Decodes a linear edge index to (u, v) for a directed simple graph (no loops).5091// O(1).5092inline std::pair<int, int> decode_directed_simple_edge(int n, long long idx) {5093 int u = idx / (n - 1);5094 int rem = idx % (n - 1);5095 return {u, rem + (rem >= u)};5096}50975098// Decodes a linear edge index according to graph mode.5099// O(log n) for undirected, O(1) for directed.5100inline std::pair<int, int>5101decode_graph_edge_index(int n, long long idx, bool directed, bool self_loops) {5102 if (directed) {5103 if (self_loops)5104 return {int(idx / n), int(idx % n)};5105 return decode_directed_simple_edge(n, idx);5106 }5107 if (self_loops)5108 return decode_undirected_loops_edge(n, idx);5109 return decode_undirected_simple_edge(n, idx);5110}51115112} // namespace detail51135114/*5115 * Graph generator.5116 *5117 * Graphs of `n` vertices labeled from 0 to n-1 and `m` edges.5118 * These are labeled graphs, that is, isomorphism is not taken into5119 * account. VWeight is the type of vertex weights, and EWeight is the type of5120 * edge weights. Generator does not generate weights. The weights are to be set5121 * in the wgraph::value.5122 */51235124template <typename VWeight, typename EWeight>5125struct wgraph : gen_base<wgraph<VWeight, EWeight>> {5126 int n_, m_; // Number of vertices and edges.5127 std::set<std::pair<int, int>> edges_; // Edges that were set.5128 bool is_directed_; // If graph is directed.5129 bool has_self_loops_; // If self-loops are allowed.51305131 // Creates graph generator with `n` vertices and `m` edges.5132 // Additionally, you can set if the graph is directed and if self loops are5133 // allowed.5134 // O(1).5135 wgraph(int n, int m, bool is_directed = false, bool has_self_loops = false)5136 : n_(n), m_(m), is_directed_(is_directed),5137 has_self_loops_(has_self_loops) {5138 tgen_ensure(n > 0, "wgraph: number of vertices must be positive");5139 }51405141 // Adds edge between u and v (this edge must be generated).5142 // O(log m).5143 wgraph &add_edge(int u, int v) {5144 tgen_ensure(0 <= std::min(u, v) and std::max(u, v) < n_,5145 "wgraph: vertices must be indexed in [0, n)");51465147 if (!is_directed_ and u > v)5148 std::swap(u, v);5149 edges_.emplace(u, v);5150 tgen_ensure(static_cast<int>(edges_.size()) <= m_,5151 "wgraph: too many edges were added");5152 return *this;5153 }51545155 // Graph value.5156 //5157 // Edges are stored in both directions (if undirected) in adjacency list,5158 // but only u < v in edge list.5159 // Optimized for performance (lazy adjacency list; edge-list constructor5160 // stores edges only).5161 struct value : gen_value_base<value> {5162 using std_type = std::tuple<int, int, std::vector<std::set<int>>>;51635164 int n_; // Number of vertices.5165 std::vector<std::set<int>> adj_; // Adjacency list.5166 std::vector<std::pair<int, int>> edges_; // Edge list.5167 bool is_directed_; // If graph is directed.5168 bool add_1_; // If should add 1 for printing vertex ids.5169 bool print_nm_; // If should print n and m.5170 mutable bool adj_built_{5171 false}; // Lazy cache: true once adj_ is built from edges_; mutable5172 // so const adj() can populate it.5173 std::optional<std::vector<VWeight>> vertex_weights_; // Vertex weights.5174 std::optional<std::vector<EWeight>>5175 edge_weights_; // Edge weights (in same order as edges_ ).51765177 // Creates value from adjacency list. The edges5178 // are considered to be directed.5179 // O(n + m).5180 value(const std::vector<std::set<int>> &adj, bool is_directed = false)5181 : n_(static_cast<int>(adj.size())), adj_(adj),5182 is_directed_(is_directed), add_1_(false), print_nm_(false),5183 adj_built_(true) {5184 for (int u = 0; u < n_; ++u)5185 for (auto v : adj[u]) {5186 tgen_ensure(5187 0 <= v and v < n_,5188 "wgraph: value: vertices must be indexed in [0, n)");5189 // Undirected adjacency is symmetric: count each edge once5190 // (canonical u <= v). Directed: every out-edge appears5191 // once.5192 if (is_directed_ or u <= v)5193 edges_.emplace_back(u, v);5194 }5195 }51965197 // Creates value from `n`, `m`, and edge list. The edges are5198 // considered to be directed.5199 // Optimized for performance (lazy adjacency list; unordered_set dedup).5200 // O(m log m).5201 value(int n, const std::vector<std::pair<int, int>> &edges = {},5202 bool is_directed = false)5203 : n_(n), edges_(), is_directed_(is_directed), add_1_(false),5204 print_nm_(false), adj_built_(false) {5205 edges_.reserve(edges.size());5206 std::unordered_set<uint64_t> seen;5207 seen.reserve(edges.size() * 2 + 1);5208 for (auto [u, v] : edges) {5209 tgen_ensure(5210 0 <= std::min(u, v) and std::max(u, v) < n,5211 "wgraph: value: vertices must be indexed in [0, n)");5212 if (!is_directed_ and u > v)5213 std::swap(u, v);5214 uint64_t key = is_directed_ ? detail::directed_edge_key(u, v)5215 : detail::undirected_edge_key(u, v);5216 if (seen.insert(key).second)5217 edges_.emplace_back(u, v);5218 }5219 }5220 value(int n, const std::set<std::pair<int, int>> &edges,5221 bool is_directed = false)5222 : value(5223 n,5224 std::vector<std::pair<int, int>>(edges.begin(), edges.end()),5225 is_directed) {}5226 value(int n, const std::initializer_list<std::pair<int, int>> &edges,5227 bool is_directed = false)5228 : value(n, std::vector<std::pair<int, int>>(edges), is_directed) {}52295230 // Creates graph from tree (undirected, same edges).5231 // O(n).5232 value(const typename wtree<VWeight, EWeight>::value &t)5233 : value(t.n(), t.edges(), false) {5234 if (t.vertex_weights().has_value()) {5235 vertex_weights_ = *t.vertex_weights();5236 }5237 if (t.edge_weights().has_value()) {5238 edge_weights_ = *t.edge_weights();5239 }5240 }52415242 // Weight type conversion.5243 // O(n + m).5244 template <typename NewVWeight, typename NewEWeight>5245 typename wgraph<NewVWeight, NewEWeight>::value5246 convert_weight_types() const {5247 tgen_ensure(!vertex_weights_.has_value() and5248 !edge_weights_.has_value(),5249 "wgraph: value: cannot convert weight type after "5250 "assigning weights");52515252 ensure_adj_built();5253 typename wgraph<NewVWeight, NewEWeight>::value new_graph(5254 adj_, is_directed_);5255 new_graph.is_directed_ = is_directed_;5256 new_graph.add_1_ = add_1_;5257 new_graph.print_nm_ = print_nm_;5258 return new_graph;5259 }52605261 // Fetches number of vertices.5262 int n() const { return n_; }52635264 // Fetches number of edges.5265 int m() const { return edges_.size(); }52665267 // Fetches if graph is directed;5268 bool is_directed() const { return is_directed_; }52695270 // Fetches a const ref. to adjacency list.5271 const std::vector<std::set<int>> &adj() const {5272 ensure_adj_built();5273 return adj_;5274 }52755276 // Fetches a const ref. to edge set.5277 const std::vector<std::pair<int, int>> &edges() const { return edges_; }52785279 // Fetches vertex weights.5280 const std::optional<std::vector<VWeight>> &vertex_weights() const {5281 return vertex_weights_;5282 }52835284 // Fetches edge weights.5285 const std::optional<std::vector<EWeight>> &edge_weights() const {5286 return edge_weights_;5287 }52885289 // Sets vertex weights.5290 // O(n + m).5291 template <typename NewVWeight = VWeight>5292 typename wgraph<NewVWeight, EWeight>::value set_vertex_weights(5293 const std::vector<NewVWeight> &vertex_weights) const {5294 tgen_ensure(static_cast<int>(vertex_weights.size()) == n(),5295 "wgraph: value: must give `n` vertex weights");52965297 auto new_graph = convert_weight_types<NewVWeight, EWeight>();5298 new_graph.vertex_weights_ = vertex_weights;5299 return new_graph;5300 }53015302 // Sets edge weights.5303 // O(n + m).5304 template <typename NewEWeight = EWeight>5305 typename wgraph<VWeight, NewEWeight>::value5306 set_edge_weights(const std::vector<NewEWeight> &edge_weights) const {5307 tgen_ensure(static_cast<int>(edge_weights.size()) == m(),5308 "wgraph: value: must give `m` edge weights");53095310 auto new_graph = convert_weight_types<VWeight, NewEWeight>();5311 new_graph.edge_weights_ = edge_weights;5312 return new_graph;5313 }53145315 // Enables edge-weighted mode before adding weighted edges5316 // incrementally. The graph must have no edges yet. O(1).5317 value &edge_weighted() {5318 tgen_ensure(m() == 0,5319 "wgraph: value: edge_weighted requires a graph with no "5320 "edges");5321 tgen_ensure(!edge_weights_.has_value(),5322 "wgraph: value: graph is already edge-weighted");53235324 edge_weights_ = std::vector<EWeight>();5325 return *this;5326 }53275328 // Adds 1 to vertex ids, for printing.5329 // O(1).5330 value &add_1() {5331 add_1_ = true;5332 return *this;5333 }53345335 // Prints `n m` on a new line before printing the edges.5336 // O(1).5337 value &print_nm() {5338 print_nm_ = true;5339 return *this;5340 }53415342 // Shuffles the graph's vertex labels (except those in `indices`,5343 // which keep their current label) and edge order. The change is5344 // applied eagerly to the underlying adjacency list, edge list,5345 // vertex weights and edge weights.5346 // O(n + m).5347 value &shuffle_except(std::set<int> indices) {5348 ensure_adj_built();5349 // Builds the relabeling: for each vertex `i`, `new_label[i]` is5350 // its new id. Vertices in `indices` keep their label; the others5351 // are permuted among themselves.5352 std::vector<int> new_label(n());5353 std::vector<int> shuffled;5354 for (int i = 0; i < n(); ++i) {5355 if (indices.count(i))5356 new_label[i] = i;5357 else5358 shuffled.push_back(i);5359 }5360 std::vector<int> targets = shuffled;5361 tgen::shuffle(targets.begin(), targets.end());5362 for (size_t k = 0; k < shuffled.size(); ++k)5363 new_label[shuffled[k]] = targets[k];53645365 // Rewrites adjacency list with new labels.5366 std::vector<std::set<int>> new_adj(n());5367 for (int u = 0; u < n(); ++u)5368 for (int v : adj_[u])5369 new_adj[new_label[u]].insert(new_label[v]);5370 adj_ = new_adj;53715372 // Rewrites edges with new labels (canonical undirected order).5373 for (auto &[u, v] : edges_) {5374 u = new_label[u];5375 v = new_label[v];5376 if (!is_directed_ and u > v)5377 std::swap(u, v);5378 }53795380 // Permutes vertex weights to match the new labels.5381 if (vertex_weights_.has_value()) {5382 std::vector<VWeight> new_vw(n());5383 for (int i = 0; i < n(); ++i)5384 new_vw[new_label[i]] = (*vertex_weights_)[i];5385 vertex_weights_ = new_vw;5386 }53875388 // Shuffles edge order, keeping edge weights aligned.53895390 std::vector<int> perm(edges_.size());5391 std::iota(perm.begin(), perm.end(), 0);5392 tgen::shuffle(perm.begin(), perm.end());53935394 std::vector<std::pair<int, int>> new_edges;5395 std::optional<std::vector<EWeight>> new_ew;5396 if (edge_weights_.has_value())5397 new_ew = std::vector<EWeight>();5398 for (int i : perm) {5399 new_edges.push_back(edges_[i]);5400 if (new_ew.has_value())5401 new_ew->push_back((*edge_weights_)[i]);5402 }54035404 edges_ = new_edges;5405 if (new_ew.has_value())5406 edge_weights_ = new_ew;54075408 return *this;5409 }54105411 // Shuffles the graph's vertices and edge order.5412 // O(n + m).5413 value &shuffle() { return shuffle_except({}); }54145415 // Adds `k` vertices to the graph (labeled n, n+1, ...n+k-1). Updates5416 // `n` accordingly.5417 // O(k) amortized.5418 value &add_vertices(int k, std::optional<std::vector<VWeight>>5419 new_vertex_weights = std::nullopt) {5420 ensure_adj_built();5421 n_ += k;5422 adj_.resize(n());5423 if (new_vertex_weights.has_value()) {5424 tgen_ensure(vertex_weights().has_value(),5425 "wgraph: value: cannot add weighted vertices to "5426 "vertex-unweighted graph");5427 tgen_ensure(5428 static_cast<int>(new_vertex_weights->size()) == k,5429 "wgraph: value: number of vertex weights must be equal "5430 "to number of added vertices");54315432 vertex_weights_->insert(vertex_weights_->end(),5433 new_vertex_weights->begin(),5434 new_vertex_weights->end());5435 } else5436 tgen_ensure(!vertex_weights().has_value(),5437 "wgraph: value: cannot add unweighted vertices to "5438 "vertex-weighted graph");54395440 return *this;5441 }54425443 // Adds edge (u, v).5444 // O(log n) amortized.5445 value &add_edge(int u, int v, std::optional<EWeight> w = std::nullopt) {5446 ensure_adj_built();5447 tgen_ensure(0 <= std::min(u, v) and std::max(u, v) < n(),5448 "wgraph: value: vertex ids must be valid");54495450 if (!is_directed() and u > v)5451 std::swap(u, v);54525453 if (adj_[u].count(v))5454 return *this;54555456 adj_[u].insert(v);5457 if (!is_directed())5458 adj_[v].insert(u);5459 edges_.emplace_back(u, v);54605461 if (w.has_value()) {5462 tgen_ensure(edge_weights().has_value(),5463 "wgraph: value: cannot add weighted edge to "5464 "edge-unweighted graph");54655466 edge_weights_->push_back(*w);5467 } else5468 tgen_ensure(!edge_weights().has_value(),5469 "wgraph: value: cannot add unweighted edge to "5470 "edge-weighted graph");54715472 return *this;5473 }54745475 // Links graph with another `rhs`, adding the edge between u (in left5476 // graph) and v (in right graph). Ids for added vertices are updated5477 // accordingly.5478 // O(rhs.n + rhs.m * log n) amortized.5479 value &link(const value &rhs, int new_u, int new_v,5480 std::optional<EWeight> new_w = std::nullopt) {5481 tgen_ensure(0 <= new_u and new_u < n() and 0 <= new_v and5482 new_v < rhs.n(),5483 "wgraph: value: vertex ids must be valid");54845485 // Edges from right-hand side.5486 int shift = n();5487 add_vertices(rhs.n(), rhs.vertex_weights());5488 for (int i = 0; i < rhs.m(); ++i) {5489 auto [u, v] = rhs.edges()[i];5490 add_edge(shift + u, shift + v,5491 rhs.edge_weights().has_value()5492 ? std::optional<EWeight>((*rhs.edge_weights())[i])5493 : std::nullopt);5494 }54955496 // New edge.5497 add_edge(new_u, shift + new_v, new_w);54985499 return *this;5500 }55015502 // Glues the graph with another `rhs` such that index_pairs[i].first is5503 // considered to be the same as index_pairs[i].second. Ids for added5504 // vertices are updated accordingly.5505 // O(rhs.n + rhs.m * log n) amortized.5506 value &glue(const value &rhs,5507 std::set<std::pair<int, int>> index_pairs) {5508 tgen_ensure(5509 is_directed() == rhs.is_directed(),5510 "wgraph: value: graphs must have the same is_directed value");55115512 // Checks validity of indices.5513 std::set<int> idx_left, idx_right;5514 std::vector<int> right_id_to_left(rhs.n(), -1);5515 for (auto [l, r] : index_pairs) {5516 tgen_ensure(5517 0 <= l and l < n() and 0 <= r and r < rhs.n(),5518 "wgraph: value: vertex indices to glue must be valid");5519 tgen_ensure(idx_left.count(l) == 0 and idx_right.count(r) == 0,5520 "wgraph: value: must not have repeated indices "5521 "on the same side to glue");55225523 idx_left.insert(l);5524 idx_right.insert(r);5525 right_id_to_left[r] = l;5526 }55275528 // Computes new ids of right vertices.5529 std::vector<int> new_right_id(rhs.n(), -1);5530 int intersection_lt = 0;5531 std::optional<std::vector<VWeight>> rhs_vertex_weights;5532 for (int i = 0; i < rhs.n(); ++i) {5533 if (right_id_to_left[i] != -1) {5534 // Is in intersection.5535 ++intersection_lt;5536 new_right_id[i] = right_id_to_left[i];5537 } else {5538 // New id.5539 new_right_id[i] = n() + i - intersection_lt;5540 if (rhs.vertex_weights().has_value()) {5541 if (!rhs_vertex_weights.has_value())5542 rhs_vertex_weights = std::vector<VWeight>();5543 rhs_vertex_weights->push_back(5544 (*rhs.vertex_weights())[i]);5545 }5546 }5547 }55485549 // Adds new vertices and edges.5550 add_vertices(rhs.n() - intersection_lt, rhs_vertex_weights);5551 for (int i = 0; i < rhs.m(); ++i) {5552 auto [u, v] = rhs.edges()[i];5553 add_edge(new_right_id[u], new_right_id[v],5554 rhs.edge_weights().has_value()5555 ? std::optional<EWeight>((*rhs.edge_weights())[i])5556 : std::nullopt);5557 }55585559 return *this;5560 }5561 value &glue(const value &rhs,5562 std::initializer_list<std::pair<int, int>> il) {5563 return glue(rhs, std::set<std::pair<int, int>>(il));5564 }55655566 // Glues the graph with another `rhs` at `indices`. That is, idx in5567 // `indices` are considered to be the same vertex. Ids for added5568 // vertices are updated accordingly.5569 // O(rhs.n + rhs.m * log n) amortized.5570 value &glue(const value &rhs, std::set<int> indices) {5571 std::set<std::pair<int, int>> index_pairs;5572 for (auto i : indices)5573 index_pairs.emplace(i, i);5574 return glue(rhs, index_pairs);5575 }5576 value &glue(const value &rhs, const std::initializer_list<int> &il) {5577 return glue(rhs, std::set<int>(il));5578 }55795580 // Disjoint union.5581 // Shifts ids from `rhs` graph by n().5582 // O(rhs.n + rhs.m * log n) amortized.5583 value &disjoint_union(const value &rhs) {5584 return glue(rhs, std::set<int>());5585 }55865587 // Computes uniformly random subgraph of graph with num_edges edges.5588 // O(n + m).5589 value &random_subgraph(int num_edges) {5590 tgen_ensure(5591 num_edges <= m(),5592 "wgraph: value: can choose at most `m` edges from graph");55935594 std::vector<std::pair<int, int>> new_edges;5595 std::optional<std::vector<EWeight>> new_edge_weights;55965597 int left = m();5598 for (int i = 0; i < m(); ++i) {5599 if (next(1, left--) <= num_edges) {5600 new_edges.push_back(edges()[i]);5601 if (edge_weights_.has_value()) {5602 if (!new_edge_weights.has_value())5603 new_edge_weights = std::vector<EWeight>();5604 new_edge_weights->push_back((*edge_weights())[i]);5605 }5606 --num_edges;5607 }5608 }56095610 edges_ = new_edges;5611 edge_weights_ = new_edge_weights;5612 rebuild_adj_from_edge_list();5613 return *this;5614 }56155616 // Computes a random (not uniform) subgraph with `num_edges` edges that5617 // keeps every connected component connected (does not increase the5618 // number of connected components).5619 // 1. Picks a spanning forest via randomized Prim.5620 // 2. Adds additional edges uniformly at random.5621 // O(n + m).5622 value &random_connected_subgraph(int num_edges) {5623 tgen_ensure(!is_directed_,5624 "wgraph: value: random_connected_subgraph is only for "5625 "undirected graphs");5626 tgen_ensure(5627 num_edges <= m(),5628 "wgraph: value: can choose at most `m` edges from graph");56295630 // Builds an incidence list: for each vertex, the (neighbor, edge5631 // index) pairs.5632 std::vector<std::vector<std::pair<int, int>>> incident(n());5633 for (int i = 0; i < m(); ++i) {5634 auto [u, v] = edges_[i];5635 incident[u].emplace_back(v, i);5636 incident[v].emplace_back(u, i);5637 }56385639 // Randomized Prim.5640 std::vector<bool> vis(n(), false);5641 std::vector<int> queue;5642 std::vector<bool> in_tree(m(), false);5643 int forest_edges = 0;56445645 for (int start = 0; start < n(); ++start) {5646 if (vis[start])5647 continue;5648 vis[start] = true;5649 queue.push_back(start);56505651 while (!queue.empty()) {5652 int i = tgen::next<int>(0, queue.size() - 1);5653 int u = queue[i];5654 std::swap(queue[i], queue.back());5655 queue.pop_back();56565657 for (auto [v, edge_idx] : incident[u]) {5658 if (!vis[v]) {5659 vis[v] = true;5660 queue.push_back(v);5661 in_tree[edge_idx] = true;5662 ++forest_edges;5663 }5664 }5665 }5666 }5667 tgen_ensure(5668 num_edges >= forest_edges,5669 "wgraph: value: random_connected_subgraph needs at least "5670 "`n - c` edges, where `c` is the number of connected "5671 "components");56725673 // Splits edge indices into forest edges and the rest.5674 std::vector<int> tree_idx, rest_idx;5675 for (int i = 0; i < m(); ++i) {5676 if (in_tree[i])5677 tree_idx.push_back(i);5678 else5679 rest_idx.push_back(i);5680 }56815682 tgen::shuffle(rest_idx.begin(), rest_idx.end());56835684 std::vector<int> chosen_idx;5685 chosen_idx.insert(chosen_idx.end(), tree_idx.begin(),5686 tree_idx.end());5687 chosen_idx.insert(chosen_idx.end(), rest_idx.begin(),5688 rest_idx.begin() + num_edges - forest_edges);56895690 detail::tgen_ensure_against_bug(5691 static_cast<int>(chosen_idx.size()) == num_edges,5692 "wgraph: value: chose a wrong number of edges");56935694 std::vector<std::pair<int, int>> new_edges;5695 std::optional<std::vector<EWeight>> new_edge_weights;5696 if (edge_weights_.has_value())5697 new_edge_weights = std::vector<EWeight>();5698 for (int i : chosen_idx) {5699 new_edges.push_back(edges_[i]);5700 if (new_edge_weights.has_value())5701 new_edge_weights->push_back((*edge_weights_)[i]);5702 }57035704 edges_ = new_edges;5705 edge_weights_ = new_edge_weights;5706 rebuild_adj_from_edge_list();5707 return *this;5708 }57095710 // Complement. Self loops are maintained.5711 // O(n^2).5712 value operator!() const {5713 tgen_ensure(!edge_weights_.has_value(),5714 "wgraph: value: cannot compute complement of "5715 "edge-weighted graph");57165717 value complement = *this;5718 complement.ensure_adj_built();5719 std::vector<std::pair<int, int>> compl_edges;5720 for (int i = 0; i < complement.n_; ++i) {5721 std::set<int> complement_adj;5722 for (int j = 0; j < complement.n_; ++j) {5723 bool add_j = false;5724 if (j == i and complement.adj_[i].count(j))5725 add_j = true;5726 if (j != i and !complement.adj_[i].count(j))5727 add_j = true;57285729 if (add_j) {5730 complement_adj.insert(j);5731 // If i > j and !is_directed(), we don't add the edge.5732 if (i <= j or complement.is_directed_) {5733 compl_edges.emplace_back(i, j);5734 }5735 }5736 }5737 std::swap(complement.adj_[i], complement_adj);5738 }5739 std::swap(complement.edges_, compl_edges);57405741 return complement;5742 }57435744 // Concatenates two values.5745 // O(N + M log N), N = n + rhs.n, M = m + rhs.m.5746 value operator+(const value &rhs) const {5747 tgen_ensure(is_directed() == rhs.is_directed(),5748 "wgraph: value: graphs must have the same "5749 "is_directed value");57505751 tgen_ensure(vertex_weights().has_value() ==5752 rhs.vertex_weights().has_value(),5753 "wgraph: value: cannot concatenate vertex-weighted "5754 "wgraph to unweighted");5755 tgen_ensure(edge_weights().has_value() ==5756 rhs.edge_weights().has_value(),5757 "wgraph: value: cannot concatenate edge-weighted "5758 "wgraph to unweighted");57595760 value concat = *this;5761 concat.glue(rhs, std::set<std::pair<int, int>>());5762 concat.add_1_ = add_1_ | rhs.add_1_;5763 concat.print_nm_ = print_nm_ | rhs.print_nm_;57645765 return concat;5766 }57675768 // Prints to std::ostream.5769 // O(n + m).5770 friend std::ostream &operator<<(std::ostream &out, const value &val) {5771 // Prints `n` and `m`.5772 if (val.print_nm_)5773 out << val.n() << " " << val.m() << '\n';57745775 // Prints vertex weights.5776 if (val.vertex_weights()) {5777 for (int i = 0; i < val.n(); ++i) {5778 if (i > 0)5779 out << " ";5780 out << (*val.vertex_weights())[i];5781 }5782 out << '\n';5783 }57845785 // Prints edges.5786 for (int i = 0; i < val.m(); ++i) {5787 auto [u, v] = val.edges()[i];5788 out << (u + val.add_1_) << " " << (v + val.add_1_);57895790 // Edge weight.5791 if (val.edge_weights().has_value())5792 out << " " << (*val.edge_weights())[i];57935794 out << '\n';5795 }57965797 return out;5798 }57995800 // Gets a std::tuple<n, m, adj> representing the value.5801 std::tuple<int, int, std::vector<std::set<int>>> to_std() const {5802 ensure_adj_built();5803 return std_type(n_, m(), adj_);5804 }58055806 private:5807 // Rebuilds adjacency from edges_ after replacing the edge list (e.g.5808 // subgraph operations).5809 // O(m log n).5810 void rebuild_adj_from_edge_list() {5811 adj_.assign(n_, {});5812 for (auto [u, v] : edges_) {5813 adj_[u].insert(v);5814 if (!is_directed_)5815 adj_[v].insert(u);5816 }5817 adj_built_ = true;5818 }58195820 // Builds adj_ from edges_ on first use.5821 // O(1) if already built; O(m log n) otherwise.5822 void ensure_adj_built() const {5823 if (adj_built_)5824 return;5825 const_cast<value *>(this)->rebuild_adj_from_edge_list();5826 }5827 };58285829 // Adds all edges from `rhs` as preset edges.5830 // O(rhs.m * log m).5831 wgraph &add_edges_from(const value &rhs) {5832 tgen_ensure(is_directed_ == rhs.is_directed(),5833 "wgraph: graphs must have the same is_directed value");58345835 for (auto [u, v] : rhs.edges())5836 add_edge(u, v);5837 return *this;5838 }58395840 // Generates graph value.5841 // Optimized for performance: dense no-preset graphs use index sampling;5842 // otherwise gen_remaining_edges.5843 // O(n + m log^2 n) expected.5844 value gen() const {5845 detail::tgen_ensure_against_bug(static_cast<int>(edges_.size()) <= m_,5846 "wgraph: too many edges were added");58475848 // All edges already added.5849 if (static_cast<int>(edges_.size()) == m_)5850 return value(n_, edges_, is_directed_);58515852 // Splits into two cases to optimize performance.58535854 // No presets and m > max_edges / 2: sample m distinct edge indices.5855 if (auto indexed = try_gen_by_edge_index())5856 return *indexed;58575858 // Otherwise: fill preset edges up to m_ with uniform random edges.5859 return gen_remaining_edges(5860 std::vector<std::pair<int, int>>(edges_.begin(), edges_.end()));5861 }58625863 // Gets a (not uniformly) random connected undirected graph.5864 // 1. Preset edges induce a spanning forest on their components.5865 // 2. Then, uniformly random edges between components are added.5866 // 3. Remaining edges are added uniformly at random.5867 // O(n + m log^2 n) expected.5868 value get_connected() const {5869 tgen_ensure(!is_directed_,5870 "wgraph: get_connected is only for undirected graphs");5871 tgen_ensure(m_ >= n_ - 1,5872 "wgraph: connected graph needs at least n - 1 edges");58735874 std::vector<std::pair<int, int>> edges;5875 edges.reserve(m_);58765877 if (edges_.empty()) {5878 if (n_ > 1) {5879 std::vector<int> prufer(n_ - 2);5880 for (int i = 0; i < n_ - 2; ++i)5881 prufer[i] = next<int>(0, n_ - 1);5882 for (auto [u, v] : detail::edges_from_prufer(std::move(prufer)))5883 edges.emplace_back(u, v);5884 }5885 } else {5886 edges.assign(edges_.begin(), edges_.end());58875888 std::vector<std::vector<int>> adj(n_);5889 for (auto [u, v] : edges_) {5890 adj[u].push_back(v);5891 adj[v].push_back(u);5892 }58935894 std::vector<int> comp_size;5895 std::vector<std::vector<int>> component_ids;5896 std::vector<bool> vis(n_, false);5897 std::queue<int> q;58985899 for (int i = 0; i < n_; ++i) {5900 if (vis[i])5901 continue;59025903 vis[i] = true;5904 q.push(i);5905 comp_size.push_back(0);5906 component_ids.emplace_back();5907 while (q.size()) {5908 int u = q.front();5909 q.pop();5910 ++comp_size.back();5911 component_ids.back().push_back(u);5912 for (int v : adj[u]) {5913 if (!vis[v]) {5914 vis[v] = true;5915 q.push(v);5916 }5917 }5918 }5919 }59205921 if (component_ids.size() > 1) {5922 std::vector<int> prufer_values =5923 many_by_distribution(component_ids.size() - 2, comp_size);5924 for (auto [u, v] :5925 detail::edges_from_prufer(std::move(prufer_values)))5926 edges.emplace_back(pick(component_ids[u]),5927 pick(component_ids[v]));5928 }5929 }59305931 return gen_remaining_edges(std::move(edges));5932 }59335934 // Gets a (not uniformly) random directed acyclic graph.5935 // 1. Randomized Kahn (uniform choice among indegree-0 vertices) yields a5936 // random topological order of the preset edges (which must be acyclic).5937 // 2. Extra edges are sampled randomly using the order.5938 // With no preset edges: sample a random graph then orient acyclically.5939 // Optimized for performance (distinct upper-triangle edge-index sampling;5940 // rejection instead of pair::distinct for preset edges).5941 // O(n + m log^2 n) expected.5942 value get_acyclic() const {5943 tgen_ensure(is_directed_,5944 "wgraph: get_acyclic is only for directed graphs");59455946 if (edges_.empty()) {5947 std::vector<int> order(n_);5948 std::iota(order.begin(), order.end(), 0);5949 for (int i = n_ - 1; i > 0; --i)5950 std::swap(order[i], order[next(0, i)]);59515952 const long long max_pairs =5953 static_cast<long long>(n_) * (n_ - 1) / 2;5954 tgen_ensure(m_ <= max_pairs,5955 "wgraph: not enough edges to generate");59565957 std::vector<std::pair<int, int>> edges;5958 edges.reserve(m_);5959 for (long long idx : distinct_range<long long>(0, max_pairs - 1)5960 .gen_list(m_)5961 .to_std()) {5962 auto [i, j] = detail::decode_undirected_simple_edge(n_, idx);5963 edges.emplace_back(order[i], order[j]);5964 }5965 return value(n_, edges, true);5966 }59675968 std::vector<std::vector<int>> adj(n_);5969 std::vector<int> indeg(n_, 0);5970 for (auto [u, v] : edges_) {5971 adj[u].push_back(v);5972 ++indeg[v];5973 }59745975 std::vector<int> available;5976 for (int i = 0; i < n_; ++i)5977 if (indeg[i] == 0)5978 available.push_back(i);59795980 // Random topological order using randomized Kahn's algorithm.5981 std::vector<int> order;5982 while (!available.empty()) {5983 int idx = next(0, static_cast<int>(available.size()) - 1);5984 int u = available[idx];5985 std::swap(available[idx], available.back());5986 available.pop_back();59875988 order.push_back(u);5989 for (int v : adj[u])5990 if (--indeg[v] == 0)5991 available.push_back(v);5992 }59935994 tgen_ensure(static_cast<int>(order.size()) == n_,5995 "wgraph: preset edges contain a directed cycle");59965997 value acyclic(n_, edges_, true);59985999 // Generates final edges.60006001 detail::tgen_ensure_against_bug(acyclic.m() <= m_,6002 "wgraph: too many edges were added");60036004 if (acyclic.m() < m_) {6005 std::vector<int> order_pos(n_);6006 for (int i = 0; i < n_; ++i)6007 order_pos[order[i]] = i;60086009 std::unordered_set<uint64_t> seen;6010 seen.reserve(m_ * 2);6011 for (auto [u, v] : acyclic.edges())6012 seen.insert(6013 detail::undirected_edge_key(order_pos[u], order_pos[v]));60146015 const long long max_pairs =6016 static_cast<long long>(n_) * (n_ - 1) / 2;6017 while (acyclic.m() < m_) {6018 std::pair<int, int> edge;6019 if (!detail::try_generate_distinct(seen, [&] {6020 long long idx = next<long long>(0, max_pairs - 1);6021 edge = detail::decode_undirected_simple_edge(n_, idx);6022 return detail::undirected_edge_key(edge.first,6023 edge.second);6024 }))6025 throw detail::error("wgraph: not enough edges to generate");6026 acyclic.add_edge(order[edge.first], order[edge.second]);6027 }6028 }60296030 return acyclic;6031 }60326033 // Generates a (not uniformly) random skewed connected graph.6034 // 1. Builds the same skewed labeled tree as wtree::gen_skewed(n,6035 // elongation)(root 0, parent(i) = wnext(i, elongation) for i >= 1).6036 // If is_directed, tree edges are oriented down the tree.6037 // 2. Adds the remaining edges: pick an endpoint u uniformly;6038 // pick k uniformly in [1, spread]; walk from u toward the root k6039 // times along tree parents to get v; add edge (v, u).6040 // If elongation is small, generates a graph with small diameter.6041 // If elongation is large, generates a graph with large diameter, with6042 // vertices 0 and n-1 being far apart.6043 // O(n + m log n) if spread is O(1);6044 // O(n log n + m log^2 n) expected otherwise.6045 static value gen_skewed(int n, int m, int elongation, int spread,6046 bool is_directed = false) {6047 tgen_ensure(6048 m >= n - 1,6049 "wgraph: skewed graph needs at least n - 1 edges to be connected");6050 tgen_ensure(spread >= 2,6051 "wgraph: gen_skewed spread must be at least 2");60526053 value skewed(n, {}, is_directed);60546055 std::vector<int> parent(n), depth(n, 0);6056 parent[0] = 0;6057 for (int i = 1; i < n; ++i) {6058 int p = wnext<int>(i, elongation);6059 parent[i] = p;6060 depth[i] = depth[p] + 1;6061 skewed.add_edge(p, i);6062 }60636064 const int extra = m - (n - 1);6065 if (extra == 0)6066 return skewed;60676068 // If spread is large, use binary lifting to find the ancestor.6069 // Otherwise, enumerate O(n * spread) ancestor edges and sample6070 // directly.6071 constexpr int naive_ancestor_spread = 20;60726073 if (spread <= naive_ancestor_spread) {6074 std::vector<std::pair<int, int>> candidates;6075 candidates.reserve(n * spread);6076 for (int u = 0; u < n; ++u) {6077 int max_k = std::min(spread, depth[u]);6078 if (max_k < 2)6079 continue;6080 int v = parent[u];6081 for (int k = 2; k <= max_k; ++k) {6082 v = parent[v];6083 candidates.emplace_back(v, u);6084 }6085 }60866087 tgen_ensure(extra <= static_cast<int>(candidates.size()),6088 "wgraph: not enough edges to generate");60896090 for (auto [v, u] : choose(candidates, extra))6091 skewed.add_edge(v, u);6092 } else {6093 // Binary lifting.6094 int lg = 1;6095 while ((1 << lg) <= n)6096 ++lg;60976098 std::vector<std::vector<int>> up(lg, std::vector<int>(n));6099 for (int v = 0; v < n; ++v)6100 up[0][v] = parent[v];6101 for (int j = 1; j < lg; ++j)6102 for (int v = 0; v < n; ++v)6103 up[j][v] = up[j - 1][up[j - 1][v]];61046105 // Creates uniform generator of edges (u, v) such that v is ancestor6106 // of u. For that, every u has depth[u]-1 choices for v, so we6107 // weight u by min(spread - 1, depth[u] - 1). After that we can6108 // just pick the ancestor uniformly.6109 std::vector<int> distribution = depth;6110 for (int &d : distribution)6111 d = std::max(0, std::min(spread - 1, d - 1));6112 weighted_sampler vertex_choice(distribution);6113 distinct extra_edges([&]() -> std::pair<int, int> {6114 int u = vertex_choice.next();6115 int k = next(2, spread);6116 int v = u;6117 for (int j = 0; j < lg; ++j)6118 if (k >> j & 1)6119 v = up[j][v];6120 return {v, u};6121 });61226123 while (skewed.m() < m) {6124 std::pair<int, int> edge;6125 try {6126 edge = extra_edges.gen();6127 } catch (const std::runtime_error &e) {6128 if (std::string(e.what()) ==6129 "tgen: distinct: no more distinct values")6130 throw detail::error(6131 "wgraph: not enough edges to generate");6132 throw e;6133 }61346135 skewed.add_edge(edge.first, edge.second);6136 }6137 }61386139 return skewed;6140 }61416142 // Generates a random bipartite graph. The first side has vertices6143 // 0 .. n1-1, the second n1 .. n1+n2-1.6144 // Uniform when connected is false (distinct cross-edge indices).6145 // When connected, bipartite Prüfer + rejection fill; not uniform over6146 // connected bipartite graphs.6147 // O(n1 + n2 + m log(n1 * n2)) expected.6148 static value gen_bipartite(int n1, int n2, int m, bool connected = false) {6149 tgen_ensure(m >= 0, "wgraph: number of edges must be nonnegative");6150 long long num_edges = 1LL * n1 * n2;6151 tgen_ensure(m <= num_edges,6152 "wgraph: bipartite graph has at most n1 * n2 edges");6153 if (connected)6154 tgen_ensure(6155 m >= n1 + n2 - 1,6156 "wgraph: connected bipartite graph needs at least n1 + n2 - 1 "6157 "edges");61586159 if (!connected) {6160 std::vector<std::pair<int, int>> edges;6161 edges.reserve(m);6162 for (long long idx : distinct_range<long long>(0, num_edges - 1)6163 .gen_list(m)6164 .to_std())6165 edges.emplace_back(static_cast<int>(idx / n2),6166 n1 + static_cast<int>(idx % n2));6167 return value(n1 + n2, std::move(edges), false);6168 }61696170 std::unordered_set<uint64_t> used_edges;6171 used_edges.reserve(m * 2);6172 std::vector<std::pair<int, int>> edges;6173 edges.reserve(m);61746175 auto pack_edge = [](int u, int v) -> uint64_t {6176 if (u > v)6177 std::swap(u, v);6178 return (static_cast<uint64_t>(u) << 32) | static_cast<uint32_t>(v);6179 };61806181 if (n1 > 0 and n2 > 0) {6182 std::vector<int> prufer(n1 + n2 - 2);6183 for (int i = 0; i < n2 - 1; ++i)6184 prufer[i] = next(0, n1 - 1);6185 for (int i = 0; i < n1 - 1; ++i)6186 prufer[n2 - 1 + i] = next(n1, n1 + n2 - 1);6187 shuffle(prufer.begin(), prufer.end());6188 for (auto [u, v] : detail::edges_from_prufer(std::move(prufer))) {6189 if (u > v)6190 std::swap(u, v);6191 if (used_edges.insert(pack_edge(u, v)).second)6192 edges.emplace_back(u, v);6193 }6194 detail::tgen_ensure_against_bug(6195 used_edges.size() == size_t(n1 + n2 - 1),6196 "wgraph: invalid bipartite spanning tree size");6197 }61986199 while (edges.size() < size_t(m)) {6200 int u = next(0, n1 - 1);6201 int v = next(n1, n1 + n2 - 1);6202 if (used_edges.insert(pack_edge(u, v)).second)6203 edges.emplace_back(u, v);6204 }62056206 return value(n1 + n2, std::move(edges), false);6207 }62086209 private:6210 // If this generator has no preset edges and m is large relative to the6211 // maximum edge count, sample by distinct edge index. Otherwise6212 // std::nullopt.6213 // Optimized for performance (index sampling instead of rejection).6214 // O(m log n).6215 std::optional<value> try_gen_by_edge_index() const {6216 if (!edges_.empty())6217 return std::nullopt;62186219 long long max_edges =6220 detail::max_graph_edges(n_, is_directed_, has_self_loops_);6221 if (m_ > max_edges)6222 throw detail::error("wgraph: not enough edges to generate");6223 if (max_edges <= 0 or 2LL * m_ <= max_edges)6224 return std::nullopt;62256226 std::vector<std::pair<int, int>> edges;6227 edges.reserve(m_);6228 for (long long idx :6229 distinct_range<long long>(0, max_edges - 1).gen_list(m_).to_std())6230 edges.push_back(detail::decode_graph_edge_index(6231 n_, idx, is_directed_, has_self_loops_));62326233 return value(n_, edges, is_directed_);6234 }62356236 // Fills `edges` up to m_ with uniform random edges not already present.6237 // Optimized for performance (uint64 edge keys + try_generate_distinct).6238 // O(m log^2 n) expected.6239 value gen_remaining_edges(std::vector<std::pair<int, int>> edges) const {6240 detail::tgen_ensure_against_bug(static_cast<int>(edges.size()) <= m_,6241 "wgraph: too many edges were added");62426243 if (static_cast<int>(edges.size()) == m_)6244 return value(n_, edges, is_directed_);62456246 edges.reserve(m_);62476248 std::unordered_set<uint64_t> seen;6249 seen.reserve(m_ * 2);6250 for (auto [u, v] : edges) {6251 if (!is_directed_ and u > v)6252 std::swap(u, v);6253 seen.insert(is_directed_ ? detail::directed_edge_key(u, v)6254 : detail::undirected_edge_key(u, v));6255 }62566257 while (static_cast<int>(edges.size()) < m_) {6258 std::pair<int, int> edge;6259 if (!detail::try_generate_distinct(seen, [&] {6260 edge = detail::get_random_graph_edge(n_, is_directed_,6261 has_self_loops_);6262 if (!is_directed_ and edge.first > edge.second)6263 std::swap(edge.first, edge.second);6264 return is_directed_ ? detail::directed_edge_key(edge.first,6265 edge.second)6266 : detail::undirected_edge_key(6267 edge.first, edge.second);6268 }))6269 throw detail::error("wgraph: not enough edges to generate");6270 edges.emplace_back(edge);6271 }62726273 return value(n_, edges, is_directed_);6274 }6275};62766277// Implementation of wtree::value constructor from wgraph.6278// O(n + m alpha(n)).6279template <typename VWeight, typename EWeight>6280wtree<VWeight, EWeight>::value::value(6281 const typename wgraph<VWeight, EWeight>::value &g)6282 : n_(g.n()), adj_(g.n()), add_1_(false), print_n_(false), dsu_(g.n()) {6283 tgen_ensure(g.n() > 0, "wtree: value: graph must have at least one vertex");6284 tgen_ensure(!g.is_directed(),6285 "wtree: value: graph must be undirected to form a tree");62866287 if (g.vertex_weights().has_value())6288 vertex_weights_ = *g.vertex_weights();6289 if (g.edge_weights().has_value())6290 edge_weights_ = std::vector<EWeight>();62916292 if (n_ == 1)6293 return;62946295 std::vector<int> order(g.m());6296 std::iota(order.begin(), order.end(), 0);6297 tgen::shuffle(order.begin(), order.end());62986299 std::vector<std::pair<int, int>> tree_edges;6300 tree_edges.reserve(n_ - 1);63016302 for (int i : order) {6303 auto [u, v] = g.edges()[i];6304 if (!dsu_.unite(u, v))6305 continue;6306 if (u > v)6307 std::swap(u, v);63086309 tree_edges.emplace_back(u, v);6310 adj_[u].insert(v);6311 adj_[v].insert(u);6312 if (edge_weights_.has_value())6313 edge_weights_->push_back((*g.edge_weights())[i]);6314 if (static_cast<int>(tree_edges.size()) == n_ - 1)6315 break;6316 }63176318 tgen_ensure(static_cast<int>(tree_edges.size()) == n_ - 1,6319 "wtree: value: graph must be connected to form a tree");63206321 edges_ = std::move(tree_edges);6322}63236324/*6325 * Other types of weighted-ness.6326 */63276328// Vertex weighted graph.6329template <typename VWeight> using vgraph = wgraph<VWeight, int>;63306331// Edge weighted graph.6332template <typename EWeight> using egraph = wgraph<int, EWeight>;63336334// Unweighted graph.6335using graph = wgraph<int, int>;63366337/*6338 * Standard graphs.6339 */63406341// Complete.6342// O(n^2).6343inline graph::value K(int n) { return graph(n, n * (n - 1) / 2).gen(); }63446345// Path.6346// Path with `n` vertices. The edges of the path are 0 and n-1.6347// If directed, edges are i -> i+1 for i in [0, n-2).6348// O(n).6349inline graph::value P(int n, bool is_directed = false) {6350 graph g(n, n - 1, is_directed);6351 for (int i = 0; i + 1 < n; ++i)6352 g.add_edge(i, i + 1);6353 return g.gen();6354}63556356// Cycle.6357// n >= 3.6358// If directed, edges are i -> (i+1) % n.6359// O(n).6360inline graph::value C(int n, bool is_directed = false) {6361 tgen_ensure(n >= 3, "graph: cycle size must be at least 3");63626363 graph g(n, n, is_directed);6364 for (int i = 0; i < n; ++i)6365 g.add_edge(i, (i + 1) % n);6366 return g.gen();6367}63686369// Complete bipartite.6370// The first side has vertices `0` to `n1-1`, the second side has vertices `n1`6371// to `n1+n2-1`.6372// O(n1 * n2).6373inline graph::value K(int n1, int n2) {6374 graph g(n1 + n2, static_cast<long long>(n1) * n2);6375 for (int i = 0; i < n1; ++i)6376 for (int j = 0; j < n2; ++j)6377 g.add_edge(i, n1 + j);6378 return g.gen();6379}63806381// Star.6382// The center is vertex 0.6383// O(n).6384inline graph::value S(int n) { return K(1, n - 1); }63856386/****************6387 * *6388 * GEOMETRY *6389 * *6390 ****************/63916392namespace geometry {63936394// Point on the plane with coordinates of type T.6395template <typename T> struct point {6396 static_assert(std::is_arithmetic_v<T>,6397 "point requires an arithmetic coordinate type");63986399 // Dot/cross product type: __int128 for T = long long, long long for other6400 // integral T, T for floating-point.6401 using product_t = std::conditional_t<6402 std::is_same_v<T, long long>, detail::i128,6403 std::conditional_t<std::is_integral_v<T>, long long, T>>;64046405 // x and y coordinates.6406 T x_, y_;64076408 // Constructs a point with coordinates x and y.6409 point(T x = 0, T y = 0) : x_(x), y_(y) {}64106411 // Returns the x coordinate.6412 T x() const { return x_; }64136414 // Returns the y coordinate.6415 T y() const { return y_; }64166417 // Equality of coordinates, with epsilon-based equality for floating-point6418 // coordinates (tolerance 1e-9).6419 static bool coord_eq(T a, T b) {6420 if constexpr (std::is_integral_v<T>)6421 return a == b;6422 constexpr T eps = T(1e-9);6423 T d = a - b;6424 return d >= -eps and d <= eps;6425 }64266427 // Lexicographic order (by x, then y).6428 bool operator<(const point &p) const {6429 if (!coord_eq(x_, p.x()))6430 return x_ < p.x();6431 return y_ < p.y();6432 }64336434 // Equality of coordinates.6435 bool operator==(const point &p) const {6436 return coord_eq(x_, p.x()) and coord_eq(y_, p.y());6437 }64386439 // Vector addition.6440 point operator+(const point &p) const {6441 return point(x_ + p.x(), y_ + p.y());6442 }64436444 // Vector subtraction.6445 point operator-(const point &p) const {6446 return point(x_ - p.x(), y_ - p.y());6447 }64486449 // Scalar multiplication.6450 point operator*(T c) const { return point(x_ * c, y_ * c); }64516452 // Dot product.6453 product_t operator*(const point &p) const {6454 if constexpr (std::is_floating_point_v<T>)6455 return x_ * p.x() + y_ * p.y();6456 return product_t(x_) * p.x() + product_t(y_) * p.y();6457 }64586459 // Cross product (signed area of the parallelogram).6460 product_t operator^(const point &p) const {6461 if constexpr (std::is_floating_point_v<T>)6462 return x_ * p.y() - y_ * p.x();6463 return product_t(x_) * p.y() - product_t(y_) * p.x();6464 }64656466 // Prints the point as "x y".6467 friend std::ostream &operator<<(std::ostream &out, const point &p) {6468 return out << p.x() << ' ' << p.y();6469 }6470};64716472// Generates n distinct integer points in [min_coord, max_coord]^2 with no three6473// collinear.6474// O(n).6475inline std::vector<point<long long>>6476random_points_general_position(int n, long long min_coord,6477 long long max_coord) {6478 tgen_ensure(n > 0,6479 "geometry: random_points_general_position: n must be positive");6480 tgen_ensure(max_coord >= min_coord,6481 "geometry: random_points_general_position: min_coord must be "6482 "at most max_coord");6483 tgen_ensure(6484 static_cast<detail::i128>(max_coord) - min_coord <=6485 std::numeric_limits<long long>::max(),6486 "geometry: random_points_general_position: coordinate range too large");6487 uint64_t width = max_coord - min_coord;6488 uint64_t p = math::prime_from(2 * n);64896490 // Requires width >= p - 1 because sheared coordinates lie in [0, p - 1].6491 tgen_ensure(width >= p - 1,6492 "geometry: random_points_general_position: coordinate range "6493 "too small for n");64946495 // Base set: (x, x^-1 mod p) for x = 1, ..., p - 1.6496 //6497 // For a line ax + by + c = 0, substituting y = x^-1 gives ax^2 + cx + b = 06498 // (for x != 0), a quadratic with at most two roots in F_p. So at most two6499 // base points lie on any line. x |-> x^-1 is bijective on {1, ..., p - 1},6500 // so all points are distinct and no three are collinear.6501 std::vector<uint64_t> x_range(p - 1);6502 std::iota(x_range.begin(), x_range.end(), 1);6503 shuffle(x_range.begin(), x_range.end());6504 std::vector<detail::i128> bx(n), by(n);6505 for (int i = 0; i < n; ++i) {6506 uint64_t x = x_range[i];6507 bx[i] = x;6508 by[i] = math::modular_inverse(x, p);6509 }65106511 // Randomize placement without breaking general position: compose elementary6512 // shears in SL(2, F_p), each either [1 r; 0 1] or [1 0; r 1] with6513 // r in {-2, -1, 1, 2} (mod p). Every shear has determinant 1, so their6514 // product is invertible. Invertible linear maps preserve collinearity, so6515 // the image still has no three collinear points.6516 const int num_shears = 8;6517 std::vector<detail::i128> lin_x = bx, lin_y = by;65186519 for (int it = 0; it < num_shears; ++it) {6520 bool vertical_shear = next(2) == 0;6521 int shear_r = pick({-2, -1, 1, 2});65226523 for (int i = 0; i < n; ++i) {6524 if (vertical_shear)6525 lin_x[i] = (lin_x[i] + shear_r * lin_y[i]) % p;6526 else6527 lin_y[i] = (lin_y[i] + shear_r * lin_x[i]) % p;65286529 if (lin_x[i] < 0)6530 lin_x[i] += p;6531 if (lin_y[i] < 0)6532 lin_y[i] += p;6533 }6534 }65356536 detail::i128 min_x = lin_x[0], max_x = lin_x[0], min_y = lin_y[0],6537 max_y = lin_y[0];6538 for (int i = 1; i < n; ++i) {6539 min_x = std::min(min_x, lin_x[i]);6540 max_x = std::max(max_x, lin_x[i]);6541 min_y = std::min(min_y, lin_y[i]);6542 max_y = std::max(max_y, lin_y[i]);6543 }65446545 long long x_shift =6546 min_coord - min_x + next<long long>(0, width - (max_x - min_x));6547 long long y_shift =6548 min_coord - min_y + next<long long>(0, width - (max_y - min_y));65496550 std::vector<point<long long>> pts;6551 for (int i = 0; i < n; ++i)6552 pts.emplace_back(lin_x[i] + x_shift, lin_y[i] + y_shift);6553 return pts;6554}65556556namespace detail {65576558using i128 = tgen::detail::i128;65596560// Signed area of triangle (a, b, p); positive iff (a, b, p) are in6561// counterclockwise order. 0 iff (a, b, p) are collinear. O(1).6562inline i128 ccw(const point<long long> &a, const point<long long> &b,6563 const point<long long> &p) {6564 return (static_cast<i128>(b.x()) - a.x()) *6565 (static_cast<i128>(p.y()) - a.y()) -6566 (static_cast<i128>(b.y()) - a.y()) *6567 (static_cast<i128>(p.x()) - a.x());6568}65696570// Integer projection of P onto line AB (A and B need not be distinct).6571inline i128 proj_on_ab(const point<long long> &P, const point<long long> &A,6572 const point<long long> &B) {6573 return (P - A) * (B - A);6574}65756576// In-place Hamiltonian path on points[left..right-1] with points[left]6577// start and points[right-1] end.6578// O(n log n) expected if points are "random", O(n^2) worst case.6579inline void conquer(std::vector<point<long long>> &points, int left,6580 int right) {6581 if (right - left <= 3)6582 return;65836584 point<long long> A = points[left], B = points[right - 1];65856586 // If all points are collinear, sort them properly and return.6587 bool all_collinear = true;6588 for (int k = left + 1; k < right - 1; ++k) {6589 if (ccw(A, B, points[k]) != 0) {6590 all_collinear = false;6591 break;6592 }6593 }6594 if (all_collinear) {6595 std::sort(points.begin() + left, points.begin() + right,6596 [&](const point<long long> &P, const point<long long> &Q) {6597 return proj_on_ab(P, A, B) < proj_on_ab(Q, A, B);6598 });6599 return;6600 }66016602 // Choses a pivot that is not collinear with A and B.6603 std::vector<int> candidates;6604 for (int k = left + 1; k < right - 1; ++k) {6605 if (ccw(A, B, points[k]) != 0)6606 candidates.push_back(k);6607 }6608 int ci = candidates[next(0, static_cast<int>(candidates.size()) - 1)];6609 point<long long> C = points[ci];66106611 uint64_t wa = next<uint64_t>(1, std::numeric_limits<uint64_t>::max());6612 uint64_t wb = next<uint64_t>(1, std::numeric_limits<uint64_t>::max());6613 bool a_on_positive = ccw(C, A, B) < 0;66146615 // Classify interior points into two sides of the wedge A-C-B for partition.6616 // Collinear points on AB are tie-broken along the segment.6617 i128 proj_sum = proj_on_ab(A, A, B) + proj_on_ab(B, A, B);6618 auto is_positive = [&](const point<long long> &P) -> bool {6619 i128 s = wa * ccw(C, A, P) + wb * ccw(C, B, P);6620 // Weighted wedge side of P w.r.t. C, A, B.6621 if (s != 0)6622 return s > 0;6623 // P is on line AB: split by projection past the midpoint.6624 return 2 * proj_on_ab(P, A, B) > proj_sum;6625 };66266627 // Holds C at points[right-2] while classifying interior points in6628 // [left+1, right-3].6629 if (ci != right - 2)6630 std::swap(points[ci], points[right - 2]);66316632 int i = left + 1;6633 int j = right - 3;6634 while (i < j) {6635 if (is_positive(points[i]) == a_on_positive)6636 ++i;6637 else if (is_positive(points[j]) != a_on_positive)6638 --j;6639 else {6640 std::swap(points[i], points[j]);6641 ++i;6642 --j;6643 }6644 }66456646 // After partition:6647 // points[left]=A | (A,C)... | C | (C,B)... | points[right-1]=B.66486649 // After the swap, p is the index of C (pivot between the two subpaths).6650 int p = i;6651 if (i == j and is_positive(points[i]) == a_on_positive)6652 ++p;6653 std::swap(points[p], points[right - 2]);66546655 // Path A -> C.6656 conquer(points, left, p + 1);6657 // Path C -> B.6658 conquer(points, p, right);6659}66606661// Samples k sorted distinct integers from [left, right] uniformly.6662// Optimized for performance (pool partial Fisher–Yates or complement path for6663// modest ranges; sparse-map fallback otherwise).6664// O(k log k); O(right - left) memory when the range is modest.6665inline std::vector<long long>6666sample_sorted_distinct_in_range(int k, long long left, long long right) {6667 long long universe = right - left + 1;6668 std::vector<long long> res;6669 res.reserve(k);6670 if (k == 0)6671 return res;66726673 constexpr long long pool_threshold = 8'000'000;6674 constexpr long long pool_always_below = 500'000;66756676 if (universe <= pool_threshold and6677 (universe <= pool_always_below or k >= universe / 4)) {6678 size_t u = universe;6679 size_t ks = k;6680 std::vector<long long> pool(u);6681 std::iota(pool.begin(), pool.end(), left);6682 size_t m = ks <= u / 2 ? ks : u - ks;6683 for (size_t i = 0; i < m; ++i) {6684 size_t j = next<size_t>(i, u - 1);6685 std::swap(pool[i], pool[j]);6686 }6687 if (ks <= u / 2) {6688 res.assign(pool.begin(), pool.begin() + ks);6689 std::sort(res.begin(), res.end());6690 } else {6691 std::vector<char> excluded(u, 0);6692 for (size_t i = 0; i < m; ++i)6693 excluded[pool[i] - left] = 1;6694 for (long long v = left; v <= right; ++v)6695 if (!excluded[v - left])6696 res.push_back(v);6697 }6698 } else {6699 std::unordered_map<long long, long long> virtual_list;6700 virtual_list.reserve(k * 2);6701 for (long long i = 0; i < k; ++i) {6702 long long j = next<long long>(i, universe - 1);6703 long long vi = virtual_list.count(i) ? virtual_list[i] : i;6704 long long vj = virtual_list.count(j) ? virtual_list[j] : j;6705 virtual_list[j] = vi;6706 virtual_list[i] = vj;6707 res.push_back(virtual_list[i] + left);6708 }6709 std::sort(res.begin(), res.end());6710 }6711 return res;6712}67136714// Valtr-style signed edge components along one axis from n sorted distinct6715// coordinates. The n differences sum to zero.6716inline std::vector<long long>6717valtr_edge_components(const std::vector<long long> &sorted_coords) {6718 int n = sorted_coords.size();6719 std::vector<long long> left, right;6720 left.reserve(n / 2);6721 right.reserve(n / 2);6722 for (int i = 1; i + 1 < n; ++i) {6723 if (next(2) == 0)6724 left.push_back(sorted_coords[i]);6725 else6726 right.push_back(sorted_coords[i]);6727 }6728 long long lo = sorted_coords.front(), hi = sorted_coords.back();6729 std::vector<long long> seq;6730 seq.reserve(n + 1);6731 seq.push_back(lo);6732 for (long long v : left)6733 seq.push_back(v);6734 seq.push_back(hi);6735 for (auto it = right.rbegin(); it != right.rend(); ++it)6736 seq.push_back(*it);6737 seq.push_back(lo);6738 std::vector<long long> comps(n);6739 for (int i = 0; i < n; ++i)6740 comps[i] = seq[i + 1] - seq[i];6741 return comps;6742}67436744// Drops boundary vertices that are collinear with their cyclic neighbors.6745// O(m), m = |points|.6746inline std::vector<point<long long>>6747simplify_strict_boundary(std::vector<point<long long>> points) {6748 int n = points.size();6749 if (n < 3)6750 return points;67516752 std::vector<point<long long>> strict_points;6753 strict_points.reserve(n);6754 for (int i = 0; i < n; ++i) {6755 if (ccw(points[(i + n - 1) % n], points[i], points[(i + 1) % n]) != 0)6756 strict_points.push_back(points[i]);6757 }6758 return strict_points;6759}67606761// Picks k evenly spaced vertices along a longer cyclic boundary.6762// O(k).6763inline std::vector<point<long long>>6764subsample_boundary(const std::vector<point<long long>> &points, int k) {6765 int n = points.size();6766 if (n <= k)6767 return points;67686769 std::vector<point<long long>> sampled_points;6770 sampled_points.reserve(k);6771 for (int i = 0; i < k; ++i)6772 sampled_points.push_back(points[(static_cast<i128>(i) * n) / k]);6773 return sampled_points;6774}67756776// Random translation so the polygon lies in the box.6777// O(|points|).6778inline void place_inside_box(std::vector<point<long long>> &points,6779 long long min_coord, long long max_coord) {6780 long long width = max_coord - min_coord + 1;67816782 i128 min_x = points[0].x(), max_x = points[0].x();6783 i128 min_y = points[0].y(), max_y = points[0].y();6784 for (const point<long long> &p : points) {6785 min_x = std::min(min_x, static_cast<i128>(p.x()));6786 max_x = std::max(max_x, static_cast<i128>(p.x()));6787 min_y = std::min(min_y, static_cast<i128>(p.y()));6788 max_y = std::max(max_y, static_cast<i128>(p.y()));6789 }67906791 i128 span_x = max_x - min_x;6792 i128 span_y = max_y - min_y;6793 // Random slack keeps the polygon inside the box without filling it.6794 i128 shift_x =6795 min_coord - min_x +6796 next<long long>(0, width - 1 - static_cast<long long>(span_x));6797 i128 shift_y =6798 min_coord - min_y +6799 next<long long>(0, width - 1 - static_cast<long long>(span_y));68006801 for (point<long long> &p : points)6802 p = point<long long>(p.x() + shift_x, p.y() + shift_y);6803}68046805// Random cyclic shift.6806// O(|points|).6807inline void randomize_cyclic_shift(std::vector<point<long long>> &points) {6808 int rot = next(points.size());6809 if (rot > 0)6810 std::rotate(points.begin(), points.begin() + rot, points.end());6811}68126813// Valtr walk for m edges; bbox minimum translated to the origin.6814// O(m log m).6815inline std::vector<point<long long>>6816valtr_vertices(int m, const std::vector<long long> &x_comp,6817 std::vector<long long> y_comp) {6818 shuffle(y_comp.begin(), y_comp.end());68196820 std::vector<point<long long>> edges(m);6821 // Upper half-plane (positive y, or y = 0 and x > 0) sorts before lower.6822 auto upper = [](const point<long long> &p) {6823 return p.y() > 0 or (p.y() == 0 and p.x() > 0);6824 };6825 for (int i = 0; i < m; ++i)6826 edges[i] = point<long long>(x_comp[i], y_comp[i]);68276828 std::sort(edges.begin(), edges.end(),6829 [&upper](const point<long long> &a, const point<long long> &b) {6830 bool au = upper(a), bu = upper(b);6831 if (au != bu)6832 return au;6833 auto cross = a ^ b;6834 if (cross != 0)6835 return cross > 0;6836 return (a * a) < (b * b);6837 });68386839 // Prefix-sum the sorted edge vectors to obtain vertex coordinates.6840 i128 cur_x = 0, cur_y = 0;6841 std::vector<i128> px(m), py(m);6842 for (int i = 0; i < m; ++i) {6843 px[i] = cur_x;6844 py[i] = cur_y;6845 cur_x += edges[i].x();6846 cur_y += edges[i].y();6847 }6848 tgen::detail::tgen_ensure_against_bug(6849 cur_x == 0 and cur_y == 0,6850 "geometry: random_convex_polygon: walk did not close");68516852 i128 min_x = px[0], min_y = py[0];6853 for (int i = 1; i < m; ++i) {6854 min_x = std::min(min_x, px[i]);6855 min_y = std::min(min_y, py[i]);6856 }68576858 // Shift so the bbox minimum is at the origin.6859 std::vector<point<long long>> points;6860 points.reserve(m);6861 for (int i = 0; i < m; ++i)6862 points.emplace_back(px[i] - min_x, py[i] - min_y);6863 return points;6864}68656866} // namespace detail68676868// Generates n vertices of a convex integer polygon inside a box.6869// If strict is true, boundary vertices are guaranteed non-collinear when6870// generation succeeds; retry count depends on n and width.6871// Always returns points in counterclockwise order.6872// O(n log n).6873inline std::vector<point<long long>>6874random_convex_polygon(int n, long long min_coord, long long max_coord,6875 bool strict = false) {6876 tgen_ensure(n >= 3,6877 "geometry: random_convex_polygon: n must be at least 3");6878 tgen_ensure(max_coord >= min_coord,6879 "geometry: random_convex_polygon: min_coord must be at most "6880 "max_coord");6881 tgen_ensure(static_cast<detail::i128>(max_coord) - min_coord + 1 <=6882 std::numeric_limits<long long>::max(),6883 "geometry: random_convex_polygon: coordinate range too large");6884 long long width = max_coord - min_coord + 1;6885 tgen_ensure(6886 width >= n,6887 "geometry: random_convex_polygon: coordinate range too small for n");68886889 // Valtr walk size: n in weak mode; strict mode uses a larger grid so6890 // collinear removal still leaves at least n vertices to subsample.6891 int num_coords = n;6892 if (strict) {6893 // Extra grid lines beyond n: at least 100 (small-n headroom), about6894 // n/1000 for large n, and never more than width - n.6895 int extra = width <= n ? 06896 : std::min<long long>(std::max(100, n / 1000),6897 width - n);6898 num_coords = n + extra;6899 }69006901 // Strict mode retries coordinate sampling when simplification leaves < n6902 // vertices; weak mode has no failure path, so one attempt always suffices.6903 const int max_attempts = strict ? 32 : 1;6904 for (int i = 0; i < max_attempts; ++i) {6905 // Build a convex lattice polygon on [0, width - 1]^2, then translate.6906 std::vector<long long> x_sorted =6907 detail::sample_sorted_distinct_in_range(num_coords, 0, width - 1);6908 std::vector<long long> y_sorted =6909 detail::sample_sorted_distinct_in_range(num_coords, 0, width - 1);6910 std::vector<long long> x_comp = detail::valtr_edge_components(x_sorted);6911 std::vector<long long> y_comp = detail::valtr_edge_components(y_sorted);69126913 std::vector<point<long long>> points =6914 detail::valtr_vertices(num_coords, x_comp, std::move(y_comp));69156916 if (strict) {6917 std::vector<point<long long>> simplified =6918 detail::simplify_strict_boundary(std::move(points));6919 // Tight boxes can leave too few vertices -> resample coordinates.6920 if (static_cast<int>(simplified.size()) < n)6921 continue;69226923 points = detail::subsample_boundary(simplified, n);6924 }69256926 detail::place_inside_box(points, min_coord, max_coord);6927 detail::randomize_cyclic_shift(points);6928 return points;6929 }69306931 // Generation failed.6932 throw tgen::detail::error(6933 "geometry: random_convex_polygon: generation failed: coordinate "6934 "range too small for n");6935}69366937// Random simple polygon through given distinct points.6938// Collinear triples are allowed; fails if all points are collinear.6939// Always returns vertices in counterclockwise order.6940// O(n log n) expected if points are "random", O(n^2) worst case.6941inline std::vector<point<long long>> random_simple_polygon_through_points(6942 const std::vector<point<long long>> &points) {6943 int n = points.size();6944 tgen_ensure(n >= 3,6945 "geometry: random_simple_polygon_through_points: need at "6946 "least 3 points");6947 tgen_ensure(6948 static_cast<int>(6949 std::set<point<long long>>(points.begin(), points.end()).size()) ==6950 n,6951 "geometry: random_simple_polygon_through_points: points must "6952 "be distinct");69536954 int idx_a = 0, idx_b = 0;6955 for (int i = 1; i < n; ++i) {6956 if (points[i] < points[idx_a])6957 idx_a = i;6958 if (points[idx_b] < points[i])6959 idx_b = i;6960 }6961 point<long long> A = points[idx_a], B = points[idx_b];69626963 bool all_collinear = true;6964 for (int i = 0; i < n; ++i) {6965 if (i == idx_a or i == idx_b)6966 continue;6967 if (detail::ccw(A, B, points[i]) != 0) {6968 all_collinear = false;6969 break;6970 }6971 }6972 tgen_ensure(!all_collinear,6973 "geometry: random_simple_polygon_through_points: all points "6974 "are collinear; no simple polygon exists");69756976 // Keep points collinear with AB on the chain that has no other points on6977 // its side, so AB is split through those vertices instead of crossing them6978 // later.6979 int negative_count = 0;6980 for (int i = 0; i < n; ++i) {6981 if (i == idx_a or i == idx_b)6982 continue;6983 if (detail::ccw(A, B, points[i]) < 0)6984 ++negative_count;6985 }69866987 std::vector<point<long long>> chain;6988 chain.push_back(A);6989 int left_count = 0;6990 for (int i = 0; i < n; ++i) {6991 if (i == idx_a or i == idx_b)6992 continue;6993 detail::i128 side = detail::ccw(A, B, points[i]);6994 if (side < 0 or (side == 0 and negative_count == 0)) {6995 chain.push_back(points[i]);6996 ++left_count;6997 }6998 }6999 chain.push_back(B);7000 for (int i = 0; i < n; ++i) {7001 if (i == idx_a or i == idx_b)7002 continue;7003 detail::i128 side = detail::ccw(A, B, points[i]);7004 if (side > 0 or (side == 0 and negative_count != 0))7005 chain.push_back(points[i]);7006 }7007 chain.push_back(A);70087009 int n1 = 2 + left_count;7010 // Upper chain: A -> B.7011 detail::conquer(chain, 0, n1);7012 // Lower chain: B -> A.7013 detail::conquer(chain, n1 - 1, chain.size());70147015 // Cyclic vertex order: chain[1..n1) then chain[n1..end) (skip each path's7016 // start vertex).7017 std::vector<point<long long>> poly;7018 poly.insert(poly.end(), chain.begin() + 1, chain.begin() + n1);7019 poly.insert(poly.end(), chain.begin() + n1, chain.end());7020 return poly;7021}70227023namespace detail {70247025// Samples n distinct integer points in [min_coord, max_coord]^2.7026// O(n log n).7027inline std::vector<point<long long>>7028random_distinct_points_in_box(int n, long long min_coord, long long max_coord) {7029 long long width = max_coord - min_coord;7030 i128 side_128 = width + 1;7031 i128 universe = side_128 * side_128;7032 tgen_ensure(universe <= std::numeric_limits<long long>::max(),7033 "geometry: random_simple_polygon: coordinate range too large");7034 long long side = side_128;7035 tgen_ensure(universe >= n,7036 "geometry: random_simple_polygon: coordinate range too small "7037 "for n distinct points");70387039 // Decodes a linear grid key (x * side + y) to a point.7040 auto decode = [&](long long key) -> point<long long> {7041 return point<long long>(min_coord + key / side, min_coord + key % side);7042 };70437044 // Repeats until all generated points are not collinear.7045 // Runs O(1) expected times.7046 while (true) {7047 std::vector<long long> keys =7048 distinct_range<long long>(0, universe - 1).gen_list(n).to_std();70497050 std::vector<point<long long>> points;7051 points.reserve(n);7052 for (long long key : keys)7053 points.push_back(decode(key));70547055 // Checks if the points are not all collinear.7056 for (int i = 2; i < n; ++i) {7057 if (ccw(points[0], points[1], points[i]) != 0)7058 return points;7059 }7060 }7061}70627063// Axis-aligned edge data for a CCW polygon (interior on the left).7064// O(1).7065struct ortho_poly_edge {7066 long long len;7067 bool horiz;7068 long long fixed;7069 long long lo, hi;7070 int out_x, out_y;7071};70727073// True if a, b, c lie on one horizontal or vertical line.7074// O(1).7075inline bool ortho_axis_collinear(const point<long long> &a,7076 const point<long long> &b,7077 const point<long long> &c) {7078 return (a.x() == b.x() and b.x() == c.x()) or7079 (a.y() == b.y() and b.y() == c.y());7080}70817082// Edge i of poly: orientation, span, exterior normal.7083// O(1).7084inline ortho_poly_edge7085ortho_analyze_edge(const std::vector<point<long long>> &poly, int i) {7086 int m = poly.size();7087 point<long long> a = poly[i], b = poly[(i + 1) % m];7088 ortho_poly_edge e{};7089 if (a.y() == b.y()) {7090 e.horiz = true;7091 e.fixed = a.y();7092 e.lo = std::min(a.x(), b.x());7093 e.hi = std::max(a.x(), b.x());7094 e.out_x = 0;7095 e.out_y = a.x() < b.x() ? -1 : 1;7096 } else {7097 e.fixed = a.x();7098 e.lo = std::min(a.y(), b.y());7099 e.hi = std::max(a.y(), b.y());7100 e.out_x = a.y() < b.y() ? 1 : -1;7101 }7102 e.len = e.hi - e.lo;7103 return e;7104}71057106// True if open axis-aligned segments (a, b) and (c, d) properly cross or7107// overlap (excluding shared endpoints).7108// O(1).7109inline bool ortho_open_seg_cross(const point<long long> &a,7110 const point<long long> &b,7111 const point<long long> &c,7112 const point<long long> &d) {7113 if (a.y() == b.y() and c.y() == d.y()) {7114 if (a.y() != c.y())7115 return false;7116 long long lo1 = std::min(a.x(), b.x()), hi1 = std::max(a.x(), b.x());7117 long long lo2 = std::min(c.x(), d.x()), hi2 = std::max(c.x(), d.x());7118 return lo1 < hi2 and lo2 < hi1;7119 }7120 if (a.x() == b.x() and c.x() == d.x()) {7121 if (a.x() != c.x())7122 return false;7123 long long lo1 = std::min(a.y(), b.y()), hi1 = std::max(a.y(), b.y());7124 long long lo2 = std::min(c.y(), d.y()), hi2 = std::max(c.y(), d.y());7125 return lo1 < hi2 and lo2 < hi1;7126 }7127 if (a.y() == b.y() and c.x() == d.x()) {7128 long long hx = a.y(), vx = c.x();7129 long long hlo = std::min(a.x(), b.x()), hhi = std::max(a.x(), b.x());7130 long long vlo = std::min(c.y(), d.y()), vhi = std::max(c.y(), d.y());7131 return hlo < vx and vx < hhi and vlo < hx and hx < vhi;7132 }7133 if (a.x() == b.x() and c.y() == d.y()) {7134 long long vx = a.x(), hy = c.y();7135 long long vlo = std::min(a.y(), b.y()), vhi = std::max(a.y(), b.y());7136 long long hlo = std::min(c.x(), d.x()), hhi = std::max(c.x(), d.x());7137 return vlo < hy and hy < vhi and hlo < vx and vx < hhi;7138 }7139 return false;7140}71417142// Integral ray-crossing point-in-polygon.7143// O(|poly|).7144inline bool ortho_point_inside(const std::vector<point<long long>> &poly,7145 point<long long> p) {7146 int m = poly.size();7147 bool inside = false;7148 for (int i = 0, j = m - 1; i < m; j = i++) {7149 point<long long> a = poly[i], b = poly[j];7150 if ((a.y() > p.y()) != (b.y() > p.y())) {7151 i128 x_cross = i128(b.x() - a.x()) * (p.y() - a.y()) -7152 i128(p.x() - a.x()) * (b.y() - a.y());7153 if ((a.y() < b.y()) ? x_cross > 0 : x_cross < 0)7154 inside = !inside;7155 }7156 }7157 return inside;7158}71597160// True if p lies strictly in the open segment (a, b).7161// O(1).7162inline bool ortho_point_strictly_interior(point<long long> p,7163 point<long long> a,7164 point<long long> b) {7165 if (a.y() == b.y()) {7166 if (p.y() != a.y())7167 return false;7168 long long lo = std::min(a.x(), b.x()), hi = std::max(a.x(), b.x());7169 return lo < p.x() and p.x() < hi;7170 }7171 if (a.x() == b.x()) {7172 if (p.x() != a.x())7173 return false;7174 long long lo = std::min(a.y(), b.y()), hi = std::max(a.y(), b.y());7175 return lo < p.y() and p.y() < hi;7176 }7177 return false;7178}71797180// True if p lies on the closed segment [a, b].7181// O(1).7182inline bool ortho_point_on_segment(point<long long> p, point<long long> a,7183 point<long long> b) {7184 return p == a or p == b or ortho_point_strictly_interior(p, a, b);7185}71867187// True if splicing add between A and B on edge_i is valid: no vertex of add7188// coincides with poly, no boundary self-contact (forward or reverse7189// T-junctions, collinear overlaps), new segments do not cross other boundary7190// edges, and inward notches keep all of add inside poly.7191// O(|poly|), assuming |add| is O(1).7192inline bool ortho_bump_valid(const std::vector<point<long long>> &poly,7193 point<long long> A, point<long long> B,7194 const std::vector<point<long long>> &add,7195 int edge_i, bool inward) {7196 int m = poly.size();71977198 for (point<long long> v : add)7199 for (point<long long> q : poly)7200 if (v == q)7201 return false;72027203 // Forward T-junction: new vertex on a non-incident edge interior, or on7204 // edge_i but outside the replaced subsegment [A, B].7205 for (point<long long> v : add) {7206 for (int j = 0; j < m; ++j) {7207 point<long long> c = poly[j], d = poly[(j + 1) % m];7208 if (j == edge_i) {7209 if (ortho_point_on_segment(v, c, d) and7210 !ortho_point_on_segment(v, A, B))7211 return false;7212 } else if (ortho_point_strictly_interior(v, c, d)) {7213 return false;7214 }7215 }7216 }72177218 auto seg_ok = [&](point<long long> s0, point<long long> s1) {7219 for (int j = 0; j < m; ++j) {7220 if (j == edge_i)7221 continue;7222 point<long long> c = poly[j], d = poly[(j + 1) % m];7223 if (ortho_open_seg_cross(s0, s1, c, d))7224 return false;7225 }7226 for (int k = 0; k < m; ++k) {7227 point<long long> q = poly[k];7228 if (q == s0 or q == s1 or q == A or q == B)7229 continue;7230 if (ortho_point_strictly_interior(q, s0, s1))7231 return false;7232 }7233 return true;7234 };72357236 point<long long> prev = A;7237 for (point<long long> v : add) {7238 if (!seg_ok(prev, v))7239 return false;7240 prev = v;7241 }7242 if (!seg_ok(prev, B))7243 return false;72447245 if (inward) {7246 for (point<long long> v : add)7247 if (!ortho_point_inside(poly, v))7248 return false;7249 }72507251 return true;7252}72537254// Splices a rectangular tab or notch on edge edge_i over [lo, hi], extending7255// depth units perpendicular to the edge.7256// O(|poly|).7257inline bool ortho_bump_edge(std::vector<point<long long>> &poly, int edge_i,7258 const ortho_poly_edge &e, long long lo,7259 long long hi, long long depth, bool inward) {7260 int m = poly.size();7261 point<long long> A = poly[edge_i], B = poly[(edge_i + 1) % m];72627263 int step_x = inward ? -e.out_x : e.out_x;7264 int step_y = inward ? -e.out_y : e.out_y;72657266 std::vector<point<long long>> add;7267 if (e.horiz) {7268 long long y = e.fixed, y2 = y + step_y * depth;7269 if (A.x() < B.x()) {7270 if (lo > A.x())7271 add.emplace_back(lo, y);7272 add.emplace_back(lo, y2);7273 add.emplace_back(hi, y2);7274 if (hi < B.x())7275 add.emplace_back(hi, y);7276 } else {7277 if (hi < A.x())7278 add.emplace_back(hi, y);7279 add.emplace_back(hi, y2);7280 add.emplace_back(lo, y2);7281 if (lo > B.x())7282 add.emplace_back(lo, y);7283 }7284 } else {7285 long long x = e.fixed, x2 = x + step_x * depth;7286 if (A.y() < B.y()) {7287 if (lo > A.y())7288 add.emplace_back(x, lo);7289 add.emplace_back(x2, lo);7290 add.emplace_back(x2, hi);7291 if (hi < B.y())7292 add.emplace_back(x, hi);7293 } else {7294 if (hi < A.y())7295 add.emplace_back(x, hi);7296 add.emplace_back(x2, hi);7297 add.emplace_back(x2, lo);7298 if (lo > B.y())7299 add.emplace_back(x, lo);7300 }7301 }7302 if (!ortho_bump_valid(poly, A, B, add, edge_i, inward))7303 return false;73047305 poly.insert(poly.begin() + edge_i + 1, add.begin(), add.end());7306 return true;7307}73087309// Picks edge i with probability proportional to7310// e.len * (4 + min(global_timestamp - last_used[i], 8)).7311// O(|poly|).7312inline int ortho_pick_poly_edge(const std::vector<point<long long>> &poly,7313 std::vector<int> &last_used, int &time_stamp) {7314 int m = poly.size();7315 if (last_used.size() != static_cast<size_t>(m)) {7316 last_used.assign(m, 0);7317 time_stamp = 0;7318 }7319 std::vector<long long> weights(m);7320 long long total = 0;7321 for (int i = 0; i < m; ++i) {7322 ortho_poly_edge e = ortho_analyze_edge(poly, i);7323 weights[i] = e.len * (4 + std::min(time_stamp - last_used[i], 8));7324 total += weights[i];7325 }7326 long long pick = next<long long>(0, total - 1);7327 for (int i = 0; i < m; ++i) {7328 pick -= weights[i];7329 if (pick < 0) {7330 last_used[i] = ++time_stamp;7331 return i;7332 }7333 }7334 last_used[m - 1] = ++time_stamp;7335 return m - 1;7336}73377338// One random inflate/cut attempt.7339// O(|poly|).7340inline bool ortho_try_bump(std::vector<point<long long>> &poly, int n,7341 std::vector<int> &last_used, int &time_stamp,7342 bool outward_only = false) {7343 if (poly.size() < 3)7344 return false;73457346 int ei = ortho_pick_poly_edge(poly, last_used, time_stamp);7347 ortho_poly_edge e = ortho_analyze_edge(poly, ei);73487349 // Edge subdivision can leave length-1 segments; a tab needs span >= 2.7350 if (e.len < 2)7351 return false;73527353 // Random subinterval [lo, lo + span] on the edge, with 2 <= span <= e.len.7354 long long span = next<long long>(2, e.len);7355 long long lo = next<long long>(e.lo, e.hi - span);73567357 // max_depth: cap on perpendicular tab/notch height (~sqrt(n), in [2, 12]).7358 // depth: actual height; shallow usually, up to max_depth 10% of the time.7359 long long max_depth =7360 std::clamp<long long>(std::sqrt(n) / 2 + 2, 2LL, 12LL);7361 long long depth =7362 next(10) == 0 ? next<long long>(std::max(2LL, max_depth / 2), max_depth)7363 : next<long long>(1, std::max(2LL, max_depth / 3));73647365 bool inward = !outward_only and next(4) == 0;7366 return ortho_bump_edge(poly, ei, e, lo, lo + span, depth, inward);7367}73687369// Drops axis-aligned collinear vertices.7370// O(n).7371inline std::vector<point<long long>>7372ortho_simplify_collinear(std::vector<point<long long>> poly) {7373 int n = poly.size();7374 if (n < 3)7375 return poly;7376 std::vector<point<long long>> out;7377 out.reserve(n);7378 for (int i = 0; i < n; ++i) {7379 if (!ortho_axis_collinear(poly[(i + n - 1) % n], poly[i],7380 poly[(i + 1) % n]))7381 out.push_back(poly[i]);7382 }7383 return out.size() >= 3 ? out : poly;7384}73857386// Removes one collinear vertex.7387// O(n).7388inline bool ortho_remove_one_collinear(std::vector<point<long long>> &poly) {7389 int n = poly.size();7390 if (n < 4)7391 return false;7392 for (int i = 0; i < n; ++i) {7393 if (!ortho_axis_collinear(poly[(i + n - 1) % n], poly[i],7394 poly[(i + 1) % n]))7395 continue;7396 poly.erase(poly.begin() + i);7397 return true;7398 }7399 return false;7400}74017402// Inserts collinear vertices on straight edges until size reaches `target` (or7403// no edge has spare integer points). Preserves boundary order.7404// O(target).7405inline void ortho_fill_collinear(std::vector<point<long long>> &poly,7406 int target) {7407 int need = target - poly.size();7408 if (need <= 0)7409 return;74107411 std::vector<point<long long>> out;7412 int m = poly.size();7413 out.reserve(poly.size() + need);74147415 for (int i = 0; i < m; ++i) {7416 point<long long> a = poly[i], b = poly[(i + 1) % m];7417 out.push_back(a);7418 if (need <= 0)7419 continue;74207421 ortho_poly_edge e = ortho_analyze_edge(poly, i);7422 long long cap = e.len - 1;7423 if (cap <= 0)7424 continue;74257426 // Adds `take` collinear vertices to the boundary.7427 long long take = std::min<long long>(need, cap);7428 bool forward = e.horiz ? a.x() < b.x() : a.y() < b.y();7429 for (long long k = 0; k < take; ++k) {7430 long long off = (k + 1) * (cap + 1) / (take + 1);7431 long long coord = forward ? e.lo + off : e.hi - off;7432 if (e.horiz)7433 out.push_back({coord, e.fixed});7434 else7435 out.push_back({e.fixed, coord});7436 }7437 need -= take;7438 }7439 poly.swap(out);7440}74417442// Inserts outward depth-1 corrugation tabs on straight edges until size reaches7443// `target` (or no edge has room for another tab). Preserves boundary order.7444// Assumes the scale-up step left every edge on a grid with spacing >= 4.7445// Each tab perturbs a coordinate by 1 and keeps a >= 2 margin from both7446// corners, so tabs can only meet tabs on the same or a directly facing edge.7447// The result is therefore simple and free of collinear triples by construction.7448// O(target).7449inline void ortho_fill_corrugation(std::vector<point<long long>> &poly,7450 size_t target) {7451 size_t n_sz = poly.size();7452 if (n_sz >= target)7453 return;74547455 size_t extra_left = target - n_sz;74567457 std::vector<point<long long>> out;7458 out.reserve(target);7459 int n = poly.size();74607461 for (int i = 0; i < n; ++i) {7462 point<long long> a = poly[i], b = poly[(i + 1) % n];7463 out.push_back(a);74647465 if (extra_left < 4)7466 continue;74677468 ortho_poly_edge e = ortho_analyze_edge(poly, i);7469 if (e.len < 5)7470 continue;74717472 long long dir = (e.horiz ? a.x() < b.x() : a.y() < b.y()) ? 1 : -1;7473 long long start = e.horiz ? a.x() : a.y();7474 long long end = e.horiz ? b.x() : b.y();74757476 for (long long pos = start + 2 * dir;7477 extra_left >= 4 and (pos - end) * dir <= -3; pos += 2 * dir) {7478 if (e.horiz) {7479 long long y2 = e.fixed + e.out_y;7480 out.push_back({pos, e.fixed});7481 out.push_back({pos, y2});7482 out.push_back({pos + dir, y2});7483 out.push_back({pos + dir, e.fixed});7484 } else {7485 long long x2 = e.fixed + e.out_x;7486 out.push_back({e.fixed, pos});7487 out.push_back({x2, pos});7488 out.push_back({x2, pos + dir});7489 out.push_back({e.fixed, pos + dir});7490 }7491 extra_left -= 4;7492 }7493 }74947495 poly.swap(out);7496}74977498// CCW square seed plus boundary inflate/cut to about n vertices.7499// O(n^2) for n <= 1000; O(n) otherwise.7500inline std::vector<point<long long>> build_orthogonal_polygon(int n,7501 bool strict) {7502 bool scale_up = n > 1000;75037504 long long side = std::max<long long>(3, std::sqrt(n));7505 if (scale_up)7506 side = std::clamp(static_cast<long long>(2 * std::sqrt(std::sqrt(n))),7507 8LL, 64LL);7508 std::vector<point<long long>> poly = {7509 {0, 0}, {side, 0}, {side, side}, {0, side}};75107511 int target_ops = scale_up ? std::max(1, static_cast<int>(4 * side - 4) / 2)7512 : std::max(1, (n - 4) / 2);75137514 if (scale_up)7515 target_ops = std::min(7516 target_ops,7517 400 + static_cast<int>(4 * std::sqrt(static_cast<double>(side))));75187519 int failure_limit = std::min(target_ops * 8, 2000);75207521 std::vector<int> last_used;7522 int time_stamp = 0, consecutive_failures = 0;7523 for (int ops = 0; ops < target_ops;) {7524 if (!scale_up and poly.size() + 2 > static_cast<size_t>(n))7525 break;75267527 if (ortho_try_bump(poly, n, last_used, time_stamp, scale_up)) {7528 ++ops;7529 consecutive_failures = 0;7530 } else if (++consecutive_failures >= failure_limit) {7531 break;7532 }7533 }7534 if (strict)7535 poly = ortho_simplify_collinear(std::move(poly));75367537 if (scale_up) {7538 while (poly.size() > static_cast<size_t>(n) and7539 ortho_remove_one_collinear(poly))7540 ;7541 long long upscale = std::max(4LL, ((n + 3) / 4 + side - 1) / side);7542 for (point<long long> &p : poly)7543 p = {p.x() * upscale, p.y() * upscale};7544 if (strict)7545 ortho_fill_corrugation(poly, n);7546 else7547 ortho_fill_collinear(poly, n);7548 } else if (!strict)7549 ortho_fill_collinear(poly, n);75507551 return poly;7552}75537554} // namespace detail75557556// Random simple polygon.7557// If strict, vertex set has no three collinear points7558// (random_points_general_position); otherwise samples distinct grid points7559// (collinear triples allowed), so it might be the case that (polygon[i],7560// polygon[i+1], and polygon[i+2]) are collinear. Polygonizes via7561// random_simple_polygon_through_points. Always counterclockwise.7562// O(n log n) expected.7563inline std::vector<point<long long>>7564random_simple_polygon(int n, long long min_coord, long long max_coord,7565 bool strict = false) {7566 tgen_ensure(n >= 3,7567 "geometry: random_simple_polygon: n must be at least 3");7568 tgen_ensure(max_coord >= min_coord,7569 "geometry: random_simple_polygon: min_coord must be at most "7570 "max_coord");7571 tgen_ensure(static_cast<detail::i128>(max_coord) - min_coord <=7572 std::numeric_limits<long long>::max(),7573 "geometry: random_simple_polygon: coordinate range too large");75747575 std::vector<point<long long>> points =7576 strict ? random_points_general_position(n, min_coord, max_coord)7577 : detail::random_distinct_points_in_box(n, min_coord, max_coord);7578 return random_simple_polygon_through_points(points);7579}75807581// Random orthogonal simple polygon, CCW. Each local bump/scale/fill step7582// preserves full simplicity, so the result is valid by construction.7583// Exactly n vertices when !strict; at most n when strict (near n for n > 1000).7584// O(n^2) for n <= 1000; O(n) otherwise.7585inline std::vector<point<long long>>7586random_orthogonal_polygon(int n, long long min_coord, long long max_coord,7587 bool strict = false) {7588 tgen_ensure(n >= 4,7589 "geometry: random_orthogonal_polygon: n must be at least 4");7590 tgen_ensure(max_coord >= min_coord,7591 "geometry: random_orthogonal_polygon: min_coord must be at "7592 "most max_coord");7593 tgen_ensure(static_cast<detail::i128>(max_coord) - min_coord + 1 <=7594 std::numeric_limits<long long>::max(),7595 "geometry: random_orthogonal_polygon: coordinate range too "7596 "large");7597 long long width = max_coord - min_coord + 1;7598 tgen_ensure(width >= 4,7599 "geometry: random_orthogonal_polygon: coordinate range too "7600 "small");76017602 long long min_side = std::max<long long>(3, std::sqrt(n));7603 if (n > 1000)7604 min_side = std::max<long long>(min_side, (n + 3) / 4);7605 tgen_ensure(min_side < width,7606 "geometry: random_orthogonal_polygon: coordinate range too "7607 "small");76087609 for (int attempt = 0; attempt < 8; ++attempt) {7610 std::vector<point<long long>> poly =7611 detail::build_orthogonal_polygon(n, strict);76127613 if (!strict and poly.size() != static_cast<size_t>(n))7614 continue;76157616 detail::i128 min_x = poly[0].x(), max_x = poly[0].x();7617 detail::i128 min_y = poly[0].y(), max_y = poly[0].y();7618 for (point<long long> p : poly) {7619 min_x = std::min(min_x, detail::i128(p.x()));7620 max_x = std::max(max_x, detail::i128(p.x()));7621 min_y = std::min(min_y, detail::i128(p.y()));7622 max_y = std::max(max_y, detail::i128(p.y()));7623 }7624 if (max_x - min_x >= width or max_y - min_y >= width)7625 continue;76267627 detail::place_inside_box(poly, min_coord, max_coord);7628 detail::randomize_cyclic_shift(poly);7629 return poly;7630 }76317632 throw tgen::detail::error(7633 "geometry: random_orthogonal_polygon: generation failed");7634}76357636} // namespace geometry76377638/************7639 * *7640 * HACK *7641 * *7642 ************/76437644namespace hack {76457646namespace detail {76477648using namespace tgen::detail;76497650// Computes polynomial hash of a string.7651// O(|s|).7652inline int hash_string(const std::string &s, int base, int mod) {7653 long long h = 0;7654 for (char c : s)7655 h = (h * base + c - 'a' + 1) % mod;7656 return h;7657}76587659// Estimates the length of the string to very likely have a collision.7660inline int estimate_length(int alphabet_size, int mod) {7661 // Magic constants.7662 double base_len = 2.5 * std::log(std::sqrt(mod));7663 double scale = std::log(alphabet_size) / std::log(2.0);7664 double adjusted = base_len / std::max(1.0, scale * 0.7);76657666 return static_cast<int>(std::ceil(adjusted));7667}76687669// Collides two strings to have the same polynomial hash.7670// O(sqrt(mod) log(mod)) with high probability.7671inline std::pair<std::string, std::string>7672birthday_attack(const std::vector<std::string> &alphabet, int base, int mod) {7673 tgen_ensure(0 < base and base < mod,7674 "birthday_attack: base must be in (0, mod)");7675 std::map<uint64_t, std::vector<int>> seen;7676 int length = estimate_length(alphabet.size(), mod);76777678 while (true) {7679 std::vector<int> seq(length);76807681 std::string s;76827683 for (int i = 0; i < length; ++i) {7684 seq[i] = next<int>(0, alphabet.size() - 1);7685 s += alphabet[seq[i]];7686 }76877688 int h = hash_string(s, base, mod);76897690 auto it = seen.find(h);7691 if (it != seen.end() and it->second != seq) {7692 std::string a, b;76937694 for (int x : it->second)7695 a += alphabet[x];7696 for (int x : seq)7697 b += alphabet[x];76987699 if (a != b)7700 return {a, b};7701 }77027703 seen[h] = seq;7704 }7705}77067707// Tried to find correct multipliers for unordered_map/set to force7708// collisions. O(1).7709inline std::set<long long> std_hash_multipliers() {7710 std::set<long long> multipliers = {85229};77117712 // Codeforces GCC GNU G++17 7.3.0 case.7713 bool codeforces_gcc_case = true;7714 if (cpp.version_ != 0 and cpp.version_ != 17)7715 codeforces_gcc_case = false;7716 if (compiler.kind_ != compiler_kind::unknown and7717 compiler.kind_ != compiler_kind::gcc)7718 codeforces_gcc_case = false;7719 if (compiler.major_ > 7)7720 codeforces_gcc_case = false;77217722 if (codeforces_gcc_case)7723 multipliers.insert(107897);77247725 return multipliers;7726}77277728} // namespace detail77297730// Fetches prefix of length n of the string "abacabadabacabae...".7731// O(n).7732inline std::string abacaba(int n) {7733 tgen_ensure(n > 0, "str: size must be positive");7734 std::string str = "a";7735 char c = 'a';7736 while (static_cast<int>(str.size()) < n) {7737 int prev_size = str.size();7738 str += ++c;7739 for (int j = 0; j < prev_size and static_cast<int>(str.size()) < n; ++j)7740 str += str[j];7741 }7742 return str;7743}77447745// Two strings that have same polynomial hash for any base, for7746// mod = power of 2 up to 2^64.7747// Thue–Morse.7748// O(1).7749inline std::pair<std::string, std::string> unsigned_polynomial_hash() {7750 std::string a, b;7751 int size = 1 << 10;7752 for (int i = 0; i < size; ++i) {7753 a += 'a' + math::detail::popcount(i) % 2;7754 b += 'a' + ('b' - a[i]);7755 }7756 return {a, b};7757}77587759// Collides two strings to have the same polynomial hash.7760// O(sqrt(mod) log(mod)) with high probability.7761// 0 < base < mod.7762inline std::pair<std::string, std::string> polynomial_hash(int alphabet_size,7763 int base, int mod) {7764 tgen_ensure(alphabet_size > 1,7765 "hack: polynomial_hash: alphabet size must be greater "7766 "than 1");7767 tgen_ensure(0 < base and base < mod,7768 "hack: polynomial_hash: base must be in (0, mod)");77697770 std::vector<std::string> alphabet(alphabet_size);7771 for (int i = 0; i < alphabet_size; ++i)7772 alphabet[i] = std::string(1, 'a' + i);7773 std::iota(alphabet.begin(), alphabet.end(), 'a');7774 return detail::birthday_attack(alphabet, base, mod);7775}77767777// Collides two strings to have the same polynomial hash for multiple bases7778// and mods (up to 2 pairs).7779// O(sqrt(mod) log^2 (mod)) with high probability,7780// with mod = max(mod_1, mod_2).7781inline std::pair<std::string, std::string>7782polynomial_hash(int alphabet_size, std::vector<int> bases,7783 std::vector<int> mods) {7784 tgen_ensure(bases.size() == mods.size(),7785 "hack: polynomial_hash: bases and mods must have the same "7786 "size");7787 tgen_ensure(bases.size() > 0,7788 "hack: polynomial_hash: must have at least one (base, mod) "7789 "pair");7790 tgen_ensure(bases.size() <= 2,7791 "hack: polynomial_hash: multi-hash hack only supported "7792 "for up to 2 (base, mod) pairs");77937794 std::vector<std::string> alphabet(alphabet_size);7795 for (int i = 0; i < alphabet_size; ++i)7796 alphabet[i] = std::string(1, 'a' + i);7797 auto [S1, T1] = detail::birthday_attack(alphabet, bases[0], mods[0]);7798 if (bases.size() == 1)7799 return {S1, T1};7800 return detail::birthday_attack({S1, T1}, bases[1], mods[1]);7801}78027803// Returns a list of integers for unordered_map/set to force collisions.7804// O(size).7805inline std::vector<long long> std_unordered(int size) {7806 tgen_ensure(size > 0, "hack: std_unordered: size must be positive");7807 std::set<long long> multipliers = detail::std_hash_multipliers();7808 long long mult = 1;7809 std::set<long long>::iterator it = multipliers.begin();78107811 std::vector<long long> list;7812 while (static_cast<int>(list.size()) < size) {7813 list.push_back(mult * (*it));7814 ++it;7815 if (it == multipliers.end()) {7816 it = multipliers.begin();7817 ++mult;7818 }7819 }7820 return list;7821}78227823// Returns queries that force \Theta(q sqrt n) asymptotic7824// for Mo algorithm for offline range queries.7825// Forces \Theta(q sqrt n) pointer moves for any ordering.7826// O(n log n + q).7827inline std::vector<std::pair<int, int>> mo_worst_case(int n, int q) {7828 std::set<std::pair<int, int>> queries;78297830 // Adversarial case.7831 int sq = std::sqrt(n);7832 for (int i = 0; i < sq; ++i) {7833 for (int j = i; j < sq; ++j) {7834 if (i * sq < n and j * sq < n)7835 queries.emplace(i * sq, j * sq);7836 }7837 }78387839 // Push extra queries.7840 for (int i = 0; i < n; ++i)7841 if (queries.size() < size_t(q)) {7842 queries.emplace(0, i);7843 queries.emplace(i, i);7844 queries.emplace(i, n - 1);7845 }78467847 std::vector<std::pair<int, int>> pool(queries.begin(), queries.end());7848 while (pool.size() < size_t(q)) {7849 int l = next(0, n - 1);7850 pool.emplace_back(l, next(l, n - 1));7851 }78527853 return choose(shuffled(pool), q);7854}78557856// Returns list of strings that have a high cost to insert in a std::set.7857// Forces cost \Theta(size log(size)).7858// Generates: {b, ab, aab, aaab, ...}.7859// O(size log(size)).7860inline std::vector<std::string> string_set_worst_case(int size) {7861 std::vector<std::string> list;7862 int k = 0, left = size;7863 while (left > 0) {7864 int cur_size = std::min(left, k + 1);7865 left -= cur_size;78667867 char right_char = cur_size == k + 1 ? 'b' : 'c';7868 list.push_back(std::string(cur_size - 1, 'a') + right_char);78697870 ++k;7871 }7872 return tgen::shuffled(list);7873}78747875// Graph for Dijkstra implementations that relax with <= instead of <.7876// Unit-weight layered graph: 0 -> {1,2}, then disjoint7877// 2x2 gadgets (i,i+1) -> {i+2,i+3} for i = 1,3,5,... Many vertices share the7878// same dist from 0; with `d + w <= dist[j]` each pop re-relaxes the whole7879// frontier below it. m = 2(n - 2) edges.7880// O(n).7881inline egraph<int>::value non_strict_relaxation_dijkstra_bug(int n) {7882 tgen_ensure(7883 n >= 3,7884 "hack: non_strict_relaxation_dijkstra_bug: needs at least 3 vertices");78857886 egraph<int>::value g(n, {}, true);7887 g.edge_weighted();7888 g.add_edge(0, 1, 1);7889 g.add_edge(0, 2, 1);7890 for (int i = 1; i + 2 < n; i += 2) {7891 g.add_edge(i, i + 2, 1);7892 if (i + 3 < n)7893 g.add_edge(i, i + 3, 1);78947895 g.add_edge(i + 1, i + 2, 1);7896 if (i + 3 < n)7897 g.add_edge(i + 1, i + 3, 1);7898 }78997900 return g.shuffle_except({0});7901}79027903// Graph for Dijkstra implementations that do not skip stale heap entries7904// (`if (d > dist[i]) continue`).7905// Hub mid = n/2: star 0 -> 1..mid-1 (weights 1..mid-1), funnel i -> mid7906// (weights 1,3,5,...), then mid -> mid+1.. (weight 1).7907// Without a stale-heap check, mid and its in-neighbors are re-popped and7908// re-relax.7909// m = n + mid - 3 edges, mid = floor(n/2).7910// O(n).7911inline egraph<int>::value stale_heap_dijkstra_bug(int n) {7912 tgen_ensure(n >= 4,7913 "hack: stale_heap_dijkstra_bug: needs at least 4 vertices");79147915 int mid = n / 2;7916 egraph<int>::value g(n, {}, true);7917 g.edge_weighted();7918 for (int i = 1; i < mid; ++i)7919 g.add_edge(0, i, i);7920 for (int i = 1; i < mid; ++i)7921 g.add_edge(i, mid, 2 * (mid - i) - 1);7922 for (int i = mid + 1; i < n; ++i)7923 g.add_edge(mid, i, 1);79247925 return g.shuffle_except({0});7926}79277928// Worst-case for FIFO-SPFA.7929// Forces Omega(n^2) from vertex 0 (Theta(n*m), m = 2n - 3).7930// Upper chain ai -> a(i+1) weight 1; lower chain bi -> b(i+1) weight 0;7931// vertical ai -> bi weight 0; cross bi -> a(i+1) weight 1. Upper chain sets7932// loose dist first; cross edges from settled bi then improve a(i+1).7933// m = 2n - 3.7934// O(n).7935inline egraph<int>::value spfa(int n) {7936 tgen_ensure(n >= 2, "hack: spfa: n must be at least 2");7937 tgen_ensure(n % 2 == 0, "hack: spfa: n must be even");79387939 egraph<int>::value g(n, {}, true);7940 g.edge_weighted();79417942 const int k = n / 2;7943 for (int i = 0; i + 1 < k; ++i)7944 g.add_edge(i, i + 1, 1);7945 for (int i = 0; i + 1 < k; ++i)7946 g.add_edge(k + i, k + i + 1, 0);7947 for (int i = 0; i < k; ++i)7948 g.add_edge(i, k + i, 0);7949 for (int i = 0; i + 1 < k; ++i)7950 g.add_edge(k + i, i + 1, 1);79517952 return g.shuffle_except({0});7953}79547955// Zadeh (1972) anti-shortest-paths flow network for Edmonds-Karp and Dinitz.7956// Source is vertex 0; sink is vertex 4l + 2k + 1.7957// n = 4l + 2k + 2, m = 6l + 4k + k^2 - 4.7958// O(l + k^2).7959inline egraph<int>::value dinitz_worst_case(int k, int l) {7960 tgen_ensure(k >= 1, "hack: dinitz_worst_case: k must be at least 1");7961 tgen_ensure(l >= 1, "hack: dinitz_worst_case: l must be at least 1");79627963 const int p1 = 2 * l - 1;7964 const int p2 = 2 * l;7965 const int q1 = 2 * l + 1;7966 const int q2 = 2 * l + 2;7967 const int n = 4 * l + 2 * k + 2;79687969 const int flow_cap = k * k * l;7970 const int layer_cap = k * k;79717972 auto a = [&](int i) { return 2 * l + 3 + 2 * i; };7973 auto b = [&](int i) { return 2 * l + 4 + 2 * i; };7974 auto t = [&](int i) { return 4 * l + 2 * k + 1 - i; };79757976 egraph<int>::value g(n, {}, true);7977 g.edge_weighted();79787979 for (int i = 0; i + 1 < 2 * l - 1; ++i)7980 g.add_edge(i, i + 1, flow_cap);7981 for (int i = 0; i + 1 < 2 * l - 1; ++i)7982 g.add_edge(t(i + 1), t(i), flow_cap);79837984 for (int i = 0; i < 2 * l - 1; i += 2) {7985 g.add_edge(i, i % 4 == 0 ? p1 : p2, layer_cap);7986 g.add_edge(i % 4 == 0 ? q1 : q2, t(i), layer_cap);7987 }79887989 for (int i = 0; i < k; ++i) {7990 g.add_edge(p1, a(i), flow_cap);7991 g.add_edge(p2, b(i), flow_cap);7992 g.add_edge(a(i), q2, flow_cap);7993 g.add_edge(b(i), q1, flow_cap);7994 }79957996 for (int i = 0; i < k; ++i)7997 for (int j = 0; j < k; ++j)7998 g.add_edge(a(i), b(j), 1);79998000 return g;8001}80028003// Returns a mask of length 19938, with weights such that xor-ing with mt199378004// outputs yields 0.8005// O(1).8006template <typename T> std::vector<bool> mt19937_xor_hash() {8007 static_assert(std::is_same_v<T, int> or std::is_same_v<T, long long>,8008 "hack: mt19937_xor_hash: T must be int or long long");80098010 constexpr std::size_t deg = 19937;80118012 std::bitset<deg + 1> a, b, c;8013 b[deg] = c[deg] = 1;8014 std::size_t l = 0, shift = 1;8015 std::mt19937 rng32;8016 std::mt19937_64 rng64;8017 for (std::size_t n = 0; n < deg * 2; ++n) {8018 a >>= 1;8019 if constexpr (std::is_same_v<T, int>)8020 a[deg] = rng32() & 1;8021 else8022 a[deg] = rng64() & 1;80238024 if ((c & a).count() % 2 == 0) {8025 ++shift;8026 continue;8027 }80288029 std::bitset<deg + 1> oc = c;8030 c ^= (b >> shift);8031 if (2 * l <= n) {8032 l = n + 1 - l;8033 b = oc;8034 shift = 1;8035 } else {8036 ++shift;8037 }8038 }80398040 std::vector<bool> mask(deg + 1);8041 for (std::size_t i = 0; i <= deg; ++i)8042 mask[i] = c[i];8043 return mask;8044}80458046// Convex polygon that breaks naive rotating calipers for maximum vertex8047// distance (advances j while dist(i, next(j)) > dist(i, j) instead of using8048// ccw).8049// O(1).8050inline std::vector<geometry::point<double>>8051naive_rotating_calipers_max_dist_bug() {8052 return {8053 {-0.9846, -1.53251}, {0.49946, 1.19525}, {0.79916, 0.98291},8054 {4.02136, -1.57843}, {3.92734, -2.37856}, {3.88558, -2.37188},8055 };8056}80578058namespace detail {80598060// Builds a hack block of order k (length fib(2k+1)).8061// O(fib(2k+1)).8062inline std::vector<int> segment_tree_beats_worst_case_block(int k) {8063 tgen_ensure(k >= 1,8064 "hack: segment_tree_beats_worst_case: k must be at least 1");80658066 std::vector<int> a(k + 1), b(k + 1);8067 std::vector<std::vector<int>> vf(k + 1), vg(k + 1);80688069 a[1] = b[1] = 1;8070 vf[1] = {1};8071 vg[1] = {1, 0};80728073 for (int i = 2; i <= k; ++i) {8074 b[i] = b[i - 1] + a[i - 1];8075 a[i] = b[i] + a[i - 1];8076 for (int x : vf[i - 1])8077 vf[i].push_back(x + a[i] + b[i]);8078 vf[i].push_back(a[i]);8079 for (int x : vg[i - 1])8080 vf[i].push_back(x + a[i]);8081 vg[i] = vf[i];8082 vg[i].push_back(0);8083 for (int x : vg[i - 1])8084 vg[i].push_back(x);8085 }80868087 vf[k].push_back(0);8088 return vf[k];8089}80908091// Appends one update round for the tiled array (offset (round * an) mod L).8092// O(fib(2k+1)).8093inline void8094segment_tree_beats_append_round(std::vector<std::vector<int>> &updates,8095 int block_len, int an, int bn, int n,8096 int round) {8097 const int off = (round * an) % block_len;8098 const int add_off = (off + block_len - bn) % block_len;8099 for (int k = 0; k < block_len; ++k) {8100 const int s = k * block_len * block_len;8101 const int sub_end = off + an;8102 if (sub_end <= block_len)8103 updates.push_back({1, s + off, s + sub_end, bn});8104 else {8105 updates.push_back({1, s + off, s + block_len, bn});8106 updates.push_back({1, s, s + (sub_end - block_len), bn});8107 }8108 const int add_end = add_off + bn;8109 if (add_end <= block_len)8110 updates.push_back({0, s + add_off, s + add_end, an});8111 else {8112 updates.push_back({0, s + add_off, s + block_len, an});8113 updates.push_back({0, s, s + (add_end - block_len), an});8114 }8115 }8116 updates.push_back({2, 0, n, an});8117 for (int k = 0; k < block_len; ++k) {8118 const int s = k * block_len * block_len;8119 updates.push_back({3, s + (off + an - 1) % block_len, 0});8120 }8121}81228123} // namespace detail81248125// Array and updates for worst case of segment tree beats.8126// O(fib(2k+1)^3 + q).8127inline std::pair<std::vector<int>, std::vector<std::vector<int>>>8128segment_tree_beats_worst_case(int k, int q) {8129 tgen_ensure(k >= 1,8130 "hack: segment_tree_beats_worst_case: k must be at least 1");8131 tgen_ensure(k <= 7, "hack: segment_tree_beats_worst_case: k too large");8132 tgen_ensure(q > 0,8133 "hack: segment_tree_beats_worst_case: q must be positive");81348135 const auto &fib = math::fibonacci();8136 const int block_len = fib[k * 2 + 1];8137 const int an = fib[k * 2];8138 const int bn = fib[k * 2 - 1];81398140 const int len = block_len;8141 const int total = len * len * len;81428143 std::vector<int> block = detail::segment_tree_beats_worst_case_block(k);8144 std::vector<int> arr(total, 0);8145 for (int x = 0; x < block_len; ++x) {8146 const int s = x * len * len;8147 for (int i = 0; i < block_len; ++i)8148 arr[s + i] = block[i];8149 }81508151 std::vector<std::vector<int>> updates;8152 updates.reserve(q);8153 const int n = total;8154 for (int round = 0; updates.size() < static_cast<std::size_t>(q); ++round) {8155 detail::segment_tree_beats_append_round(updates, block_len, an, bn, n,8156 round);8157 if (updates.size() > static_cast<std::size_t>(q))8158 updates.resize(q);8159 }8160 return {arr, updates};8161}81628163} // namespace hack81648165/*********************8166 * *8167 * MISCELLANEOUS *8168 * *8169 *********************/81708171namespace misc {81728173// Generates a uniformly random balanced parentheses sequence with k '(' and k8174// ')'. Valid means that for no prefix there are more ')' than '('.8175// O(size).8176inline std::string gen_parenthesis(int size) {8177 tgen_ensure(size > 0 and size % 2 == 0,8178 "misc: parenthesis: size must be a positive even number");81798180 int k = size / 2;8181 std::string s;8182 int open = 0, close = 0;81838184 for (int i = 0; i < size; ++i) {8185 if (open == k) {8186 s += ')';8187 ++close;8188 continue;8189 }8190 if (open == close) {8191 s += '(';8192 ++open;8193 continue;8194 }81958196 long long a = k - open, b = k - close, h = open - close;81978198 // Probability of placing '(':8199 // P('(') = (k - open) * (h + 2) / ((k - open + k - close) * (h + 1))8200 // Derived from ballot numbers ratio.8201 long long num = a * (h + 2);8202 long long den = (a + b) * (h + 1);82038204 if (next<long long>(1, den) <= num) {8205 s += '(';8206 ++open;8207 } else {8208 s += ')';8209 ++close;8210 }8211 }82128213 return s;8214}82158216} // namespace misc82178218} // namespace tgen