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 p = gen_prime(detail::kth_root_floor(left, root),2923						   detail::kth_root_floor(right, root));2924	return *detail::expo(p, root, right);2925}29262927// O(|mods| + log (right)).2928// |rems| = |mods|.2929// rems_i < mods_i.2930inline uint64_t gen_congruent(uint64_t left, uint64_t right,2931							  std::vector<uint64_t> rems,2932							  std::vector<uint64_t> mods) {2933	if (left > right)2934		throw detail::there_is_no_in_range_error("congruent number", left,2935												 right);2936	tgen_ensure(rems.size() == mods.size(),2937				"math: number of remainders and mods must be the same");2938	tgen_ensure(rems.size() > 0, "math: must have at least one congruence");29392940	detail::crt crt;2941	for (int i = 0; i < static_cast<int>(rems.size()); ++i) {2942		tgen_ensure(rems[i] < mods[i],2943					"math: remainder must be smaller than the mod");2944		crt = crt * detail::crt(rems[i], mods[i]);29452946		if (crt.a == -1)2947			throw detail::there_is_no_in_range_error("congruent number", left,2948													 right);2949		if (crt.m > right) {2950			if (!(left <= crt.a and crt.a <= right))2951				throw detail::there_is_no_in_range_error("congruent number",2952														 left, right);29532954			for (int j = 0; j < static_cast<int>(rems.size()); ++j)2955				if (crt.a % mods[j] != rems[j])2956					throw detail::there_is_no_in_range_error("congruent number",2957															 left, right);2958			return crt.a;2959		}2960	}29612962	uint64_t k_min = crt.a >= left ? 0 : ((left - crt.a) + crt.m - 1) / crt.m;2963	uint64_t k_max = (right - crt.a) / crt.m;29642965	if (k_min > k_max)2966		throw detail::there_is_no_in_range_error("congruent number", left,2967												 right);29682969	return crt.a + next(k_min, k_max) * crt.m;2970}29712972// O(log (right)).2973// rem < mod.2974inline uint64_t gen_congruent(uint64_t left, uint64_t right, uint64_t rem,2975							  uint64_t mod) {2976	return gen_congruent(left, right, std::vector<uint64_t>({rem}),2977						 std::vector<uint64_t>({mod}));2978}29792980// First congruent number >= left.2981// O(|mods| + log (left)).2982// |rems| = |mods|.2983// rems_i < mods_i.2984inline uint64_t congruent_from(uint64_t left, std::vector<uint64_t> rems,2985							   std::vector<uint64_t> mods) {2986	tgen_ensure(rems.size() == mods.size(),2987				"math: number of remainders and mods must be the same");2988	tgen_ensure(rems.size() > 0, "math: must have at least one congruence");29892990	detail::crt crt;2991	for (int i = 0; i < static_cast<int>(rems.size()); ++i) {2992		tgen_ensure(rems[i] < mods[i],2993					"math: remainder must be smaller than the mod");2994		crt = crt * detail::crt(rems[i], mods[i]);29952996		if (crt.a == -1)2997			throw detail::there_is_no_from_error("congruent number", left);2998		if (crt.m > std::numeric_limits<uint64_t>::max()) {2999			if (crt.a < left)3000				throw detail::error(3001					"math: congruent number does not exist or is too large");30023003			for (int j = 0; j < static_cast<int>(rems.size()); ++j)3004				if (crt.a % mods[j] != rems[j])3005					throw detail::error("math: congruent number does "3006										"not exist or is too large");3007			return crt.a;3008		}3009	}30103011	uint64_t k = 0;3012	if (crt.a < left)3013		k = ((left - crt.a) + crt.m - 1) / crt.m;3014	detail::i128 result = crt.a + k * crt.m;30153016	if (result > std::numeric_limits<uint64_t>::max())3017		throw detail::error("math: congruent number is too large");3018	return result;3019}30203021// O(log (left))3022// rem < mod.3023inline uint64_t congruent_from(uint64_t left, uint64_t rem, uint64_t mod) {3024	return congruent_from(left, std::vector<uint64_t>{rem},3025						  std::vector<uint64_t>{mod});3026}30273028// Last congruent number <= right.3029// O(|mods| + log (right)).3030// |rems| = |mods|.3031// rems_i < mods_i.3032inline uint64_t congruent_upto(uint64_t right, std::vector<uint64_t> rems,3033							   std::vector<uint64_t> mods) {3034	tgen_ensure(rems.size() == mods.size(),3035				"math: number of remainders and mods must be the same");3036	tgen_ensure(rems.size() > 0, "math: must have at least one congruence");30373038	detail::crt crt;3039	for (int i = 0; i < static_cast<int>(rems.size()); ++i) {3040		tgen_ensure(rems[i] < mods[i],3041					"math: remainder must be smaller than the mod");30423043		crt = crt * detail::crt(rems[i], mods[i]);30443045		if (crt.a == -1)3046			throw detail::there_is_no_upto_error("congruent number", right);3047		if (crt.m > right) {3048			if (!(crt.a <= right))3049				throw detail::there_is_no_upto_error("congruent number", right);30503051			for (int j = 0; j < static_cast<int>(rems.size()); ++j)3052				if (crt.a % mods[j] != rems[j])3053					throw detail::there_is_no_upto_error("congruent number",3054														 right);3055			return crt.a;3056		}3057	}30583059	if (crt.a > right)3060		throw detail::there_is_no_upto_error("congruent number", right);30613062	uint64_t k = (right - crt.a) / crt.m;3063	detail::i128 result = crt.a + k * crt.m;30643065	if (result < 0)3066		throw detail::there_is_no_upto_error("congruent number", right);3067	return result;3068}30693070// O(log r)3071// rem < mod.3072inline uint64_t congruent_upto(uint64_t right, uint64_t rem, uint64_t mod) {3073	return congruent_upto(right, std::vector<uint64_t>{rem},3074						  std::vector<uint64_t>{mod});3075}30763077// Mod used for FFT/NTT.3078inline constexpr int FFT_MOD = 998244353;30793080// Fibonacci sequence up to 2^64.3081inline const std::vector<uint64_t> &fibonacci() {3082	static const std::vector<uint64_t> fib = [] {3083		std::vector<uint64_t> v = {0, 1};3084		while (v.back() <=3085			   std::numeric_limits<uint64_t>::max() - v[v.size() - 2])3086			v.push_back(v.back() + v[v.size() - 2]);3087		return v;3088	}();3089	return fib;3090}30913092// Partition is ordered (composition), that is, (1, 1, 2) != (1, 2, 1).3093// O(n).3094// 0 < n.3095// 0 < part_left.3096inline std::vector<int>3097gen_partition(int n, int part_left = 1,3098			  std::optional<int> part_right = std::nullopt) {3099	if (!part_right.has_value())3100		part_right = n;3101	part_right = std::min(*part_right, n);3102	tgen_ensure(n > 0 and part_left > 0,3103				"math: invalid parameters to gen_partition");3104	tgen_ensure(part_left <= n and *part_right > 0, "math: no such partition");31053106	// dp[i] = log(number of ways to add to i).3107	std::vector<long double> dp(n + 1, detail::LOG_ZERO);3108	dp[0] = detail::LOG_ONE;3109	long double window = detail::LOG_ZERO;3110	for (int i = 1; i <= n; ++i) {3111		if (i >= part_left)3112			window = detail::add_log_space(window, dp[i - part_left]);3113		if (i >= *part_right + 1)3114			window = detail::sub_log_space(window, dp[i - *part_right - 1]);3115		dp[i] = window;3116	}3117	tgen_ensure(dp[n] >= 0, "math: no such partition");31183119	// Crazy math tricks ahead.3120	auto dp_pref = dp;3121	for (int i = 1; i <= n; ++i)3122		dp_pref[i] = detail::add_log_space(dp_pref[i - 1], dp[i]);31233124	std::vector<int> part;3125	int sum = n;3126	while (sum > 0) {3127		// Will generate a number such that what remains is in [l, r].3128		int l = std::max(0, sum - *part_right), r = sum - part_left;3129		detail::tgen_ensure_against_bug(r >= 0, "math: r < 0 in gen_partition");31303131		int nxt_sum = std::min(sum, r);3132		long double random = next<long double>(0, 1);31333134		// We generate a value X (log space), and then choose nxt_sum such3135		// that dp_pref[nxt_sum-1] < X <= dp_pref[nxt_sum].31363137		// Math hack:3138		// Let A = pref[l-1], B = pref[r], U = rand().3139		// X = log[exp(A) + U * (exp(B) - exp(A))]3140		//   = log{exp(B) * [exp(A) / exp(B) + U * (1 - exp(A) / exp(B))]}3141		//   = B + log[exp(A - B) + U - U * exp(A - B))]3142		//   = B + log[U + (1 - U) * exp(A - B)].3143		long double val_l = l ? dp_pref[l - 1] : detail::LOG_ZERO,3144					val_r = dp_pref[r];3145		while (nxt_sum > l and3146			   dp_pref[nxt_sum - 1] >=3147				   val_r + detail::log_space(random +3148											 (1 - random) * exp(val_l - val_r)))3149			--nxt_sum;31503151		part.push_back(sum - nxt_sum);3152		sum = nxt_sum;3153	}31543155	return part;3156}31573158// Partition is ordered (composition), that is, (1, 1, 2) != (1, 2, 1).3159// O(n) time/memory if part_right is not set, O(n * k) time/memory otherwise.3160// 0 < k <= n.3161// 0 <= part_left.3162inline std::vector<int>3163gen_partition_fixed_size(int n, int k, int part_left = 0,3164						 std::optional<int> part_right = std::nullopt) {3165	if (!part_right.has_value())3166		part_right = n;3167	part_right = std::min(*part_right, n);3168	tgen_ensure(0 < k and k <= n and part_left >= 0,3169				"math: invalid parameters to gen_partition_fixed_size");3170	tgen_ensure(static_cast<long long>(k) * part_left <= n and3171					n <= static_cast<long long>(k) * (*part_right),3172				"math: no such partition");31733174	// What we need to distribute to the parts.3175	int s = n - k * part_left;31763177	std::vector<int> part(k);3178	if (*part_right == n) {3179		// Stars and bars - O(n).3180		std::vector<int> cuts = {-1};31813182		int total = s + k - 1, bars = k - 1;3183		for (int i = 0; i < total and bars > 0; ++i)3184			if (next<long double>(0, 1) <3185				static_cast<long double>(bars) / (total - i)) {3186				cuts.push_back(i);3187				--bars;3188			}3189		cuts.push_back(total);31903191		// Recovers parts.3192		for (int i = 0; i < k; ++i)3193			part[i] = cuts[i + 1] - cuts[i] - 1;3194	} else {3195		// DP with log trick - O(nk).3196		int u = *part_right - part_left;31973198		// dp[i][j] = log(#ways to fill i parts with sum j)3199		std::vector<std::vector<long double>> dp(3200			k + 1, std::vector<long double>(s + 1, detail::LOG_ZERO));3201		dp[0][0] = detail::LOG_ONE;32023203		for (int i = 1; i <= k; ++i) {3204			std::vector<long double> pref = dp[i - 1];3205			for (int j = 1; j <= s; ++j)3206				pref[j] = detail::add_log_space(pref[j - 1], dp[i - 1][j]);32073208			for (int j = 0; j <= s; ++j) {3209				dp[i][j] = pref[j];3210				if (j >= u + 1)3211					dp[i][j] = detail::sub_log_space(dp[i][j], pref[j - u - 1]);3212			}3213		}32143215		// Recovers parts backwards.3216		int left_to_distribute = s;3217		for (int i = k; i >= 1; --i) {3218			long double log_total = detail::LOG_ZERO;3219			for (int j = 0; j <= u and j <= left_to_distribute; ++j)3220				log_total = detail::add_log_space(3221					log_total, dp[i - 1][left_to_distribute - j]);3222			detail::tgen_ensure_against_bug(3223				log_total != detail::LOG_ZERO,3224				"math: total == 0 in gen_partition_fixed_size");32253226			// Now we choose a number with probability proportional to3227			// dp[i-1][.].32283229			// log(rand() * total) = log(rand()) + log(total).3230			long double random =3231				detail::log_space(next<long double>(0, 1)) + log_total;32323233			long double cur_prob = detail::LOG_ZERO;3234			int chosen = 0;3235			for (int j = 0; j <= u and j <= left_to_distribute; ++j) {3236				cur_prob = detail::add_log_space(3237					cur_prob, dp[i - 1][left_to_distribute - j]);3238				if (random < cur_prob) {3239					chosen = j;3240					break;3241				}3242			}32433244			part[k - i] = chosen;3245			left_to_distribute -= chosen;3246		}3247	}32483249	for (int &i : part)3250		i += part_left;3251	return part;3252}32533254// Partition is ordered (composition), that is, (1, 1, 2) != (1, 2, 1).3255// Inspired by jngen rndm.partition: random delimiters, sort, gap recovery;3256// omits jngen's part reordering, shuffles, and two-pass redistribution.3257// 0 < k <= n.3258// 0 <= part_left.3259// Not uniformly random; optimized for speed.3260// O(k log k).3261inline std::vector<uint64_t> gen_partition_fixed_size_fast(3262	uint64_t n, int k, uint64_t part_left = 0,3263	std::optional<uint64_t> part_right = std::nullopt) {3264	if (!part_right.has_value())3265		part_right = n;3266	part_right = std::min(*part_right, n);32673268	detail::u128 n128 = n;3269	detail::u128 k128 = k;3270	detail::u128 part_left128 = part_left;3271	detail::u128 part_right128 = *part_right;32723273	tgen_ensure(k > 0 and k128 <= n128,3274				"math: invalid parameters to gen_partition_fixed_size_fast");3275	tgen_ensure(part_right128 >= part_left128 and3276					k128 * part_left128 <= n128 and3277					k128 * part_right128 >= n128,3278				"math: no such partition");32793280	uint64_t slack_total = n128 - k128 * part_left128;3281	uint64_t slack_max = part_right128 - part_left128;32823283	std::vector<uint64_t> part(k);3284	if (k == 1) {3285		part[0] = slack_total;3286	} else {3287		std::vector<uint64_t> cuts(k - 1);3288		for (uint64_t &d : cuts)3289			d = next<uint64_t>(0, slack_total);3290		std::sort(cuts.begin(), cuts.end());32913292		uint64_t prev = 0;3293		for (int i = 0; i + 1 < k; ++i) {3294			part[i] = cuts[i] - prev;3295			prev = cuts[i];3296		}3297		part[k - 1] = slack_total - prev;3298	}32993300	auto add_part_left = [part_left](uint64_t x) -> uint64_t {3301		detail::u128 val = x + part_left;3302		detail::tgen_ensure_against_bug(3303			val <= std::numeric_limits<uint64_t>::max(),3304			"math: part + part_left exceeds uint64_t in "3305			"gen_partition_fixed_size_fast");3306		return val;3307	};33083309	if (slack_max >= slack_total) {3310		for (uint64_t &x : part)3311			x = add_part_left(x);3312		return part;3313	}33143315	detail::u128 remaining = 0;3316	for (uint64_t &x : part) {3317		if (x > slack_max) {3318			remaining += x - slack_max;3319			x = slack_max;3320		}3321		x = add_part_left(x);3322	}33233324	if (remaining > 0) {3325		for (uint64_t &x : part) {3326			if (x < *part_right && remaining > 0) {3327				detail::u128 room = *part_right - x;3328				detail::u128 add = std::min(remaining, room);3329				detail::u128 val = x + add;3330				detail::tgen_ensure_against_bug(3331					val <= *part_right,3332					"math: part exceeds part_right after redistribution in "3333					"gen_partition_fixed_size_fast");3334				x = val;3335				remaining -= add;3336			}3337		}3338		detail::tgen_ensure_against_bug(3339			remaining == 0, "math: remaining mass after redistribution in "3340							"gen_partition_fixed_size_fast");3341	}33423343	return part;3344}33453346// Random partition of elements into k ordered groups (input order preserved).3347// If max_size is unset, part sizes are uniform via gen_partition_fixed_size.3348// If max_size is set, uses gen_partition_fixed_size_fast (not uniform).3349// O(n) if max_size is unset; O(n + k log k) if max_size is set.3350template <typename T>3351std::vector<std::vector<T>>3352partition_elements(std::vector<T> elements, int k, int min_size = 0,3353				   std::optional<uint64_t> max_size = std::nullopt) {3354	size_t n = elements.size();3355	tgen_ensure(k > 0, "math: partition_elements: k must be positive");3356	tgen_ensure(min_size >= 0,3357				"math: partition_elements: min_size must be non-negative");33583359	std::vector<uint64_t> sizes;3360	if (max_size.has_value()) {3361		sizes = gen_partition_fixed_size_fast(n, k, min_size, max_size);3362	} else {3363		for (int sz : gen_partition_fixed_size(n, k, min_size))3364			sizes.push_back(sz);3365	}33663367	std::vector<std::vector<T>> groups;3368	groups.reserve(k);3369	size_t pos = 0;3370	for (uint64_t sz : sizes) {3371		groups.emplace_back(elements.begin() + pos,3372							elements.begin() + pos + sz);3373		pos += sz;3374	}3375	return groups;3376}33773378}; // namespace math33793380/**************3381 *            *3382 *   STRING   *3383 *            *3384 **************/33853386namespace detail {33873388/*3389 * Regex.3390 *3391 * Compatible with testlib's regex.3392 *3393 * Operations:3394 * - A single character yields itself ("a", "3").3395 * - A list of characters inside square braces yields any a random element3396 *   from the list ("[abc123]").3397 * - A range of characters is equivalent to listing them ("[a-z1-9A-Z]").3398 * - A pattern followed by {n} yields the pattern repeated n times ("a{3}").3399 * - A pattern followed by {l,r} yields the pattern repeated between l and r3400 *   times, uniformly at random ("a{3,5}").3401 * - A list of patterns separated by | yields a random pattern from the3402 *   list, uniformly at random ("abc|def|ghi").3403 * - Parentheses can be used for grouping ("a((a|b){3})").3404 *3405 * Examples:3406 * 1. str("[1-9][0-9]{1,2}") generates two- or three-digit numbers.3407 * 2. str("a[b-d]{2}|e") generates "e" or a random string of length 3, with3408 *                       the first character being 'a' and the second and3409 *                       third characters being 'b', 'c', or 'd'.3410 * 3. str("[1-9][0-9]{%d}", n-1) generates n-digit numbers.3411 *3412 * Operations defined by {n} and {l,r} are applied from left to right, and3413 * the pattern that comes before has its delimiters defined either by () or3414 * [] at its end or is taken from the beginning of the pattern (in3415 * "a[bc]{2}", "{2}" is applied to "[bc]", and in "[01]abc{3}", the "{3}" is3416 * applied to "[01]abc").3417 */34183419// If it has children, it is either a SEQ or an OR group, defined by the3420// pattern_ field.3421struct regex_node {3422	// Considered to be repetition of left_bound != -1, pattern if3423	// children_.empty(), otherwise "SEQ" or "OR", defined by the pattern_3424	// field.3425	std::string3426		pattern_; // Either pattern, or "SEQ" or "OR" (if !children_.empty()).3427	std::vector<regex_node> children_; // Children, when SEQ or OR.3428	int left_bound_, right_bound_; // Left and right bounds of the repetition,3429								   // or -1 if not a repetition.3430	double3431		log_space_num_ways_; // Log space number of ways to match the pattern.3432	std::optional<distinct_container<char>>3433		distinct_; // Distinct generator for the pattern, for [chars].34343435	// c or [chars].3436	regex_node(const std::string &pattern)3437		: pattern_(pattern), left_bound_(-1), right_bound_(-1) {3438		if (pattern.size() == 1) {3439			log_space_num_ways_ = math::detail::LOG_ONE;3440			return;3441		}3442		tgen_ensure_against_bug(pattern[0] == '[' and pattern.back() == ']',3443								"str: invalid regex: expected character class");3444		int size = pattern.size() - 2;3445		log_space_num_ways_ = math::detail::log_space(size);3446		distinct_ = distinct_container<char>(pattern.substr(1, size));3447	}3448	// SEQ or OR.3449	regex_node(const std::string &pattern, std::vector<regex_node> &children)3450		: pattern_(pattern), left_bound_(-1), right_bound_(-1) {3451		if (pattern == "SEQ") {3452			// Multiply the number of ways.3453			log_space_num_ways_ = math::detail::LOG_ONE;3454			for (const auto &child : children)3455				log_space_num_ways_ += child.log_space_num_ways_;3456		} else if (pattern == "OR") {3457			// Add the number of ways.3458			log_space_num_ways_ = math::detail::LOG_ZERO;3459			for (const auto &child : children)3460				log_space_num_ways_ = math::detail::add_log_space(3461					log_space_num_ways_, child.log_space_num_ways_);3462		} else3463			tgen_ensure_against_bug("str: invalid regex: expected SEQ or OR");34643465		children_ = std::move(children);3466		children.clear();3467	}3468	// REP.3469	regex_node(int left_bound, int right_bound, regex_node &child)3470		: pattern_("REP"), left_bound_(left_bound), right_bound_(right_bound) {3471		log_space_num_ways_ = math::detail::LOG_ZERO;3472		for (int i = left_bound; i <= right_bound; ++i)3473			log_space_num_ways_ = math::detail::add_log_space(3474				log_space_num_ways_, i * child.log_space_num_ways_);34753476		children_.push_back(std::move(child));3477	}3478};34793480// State of the regex parser.3481struct regex_state {3482	std::vector<regex_node> cur;	  // Current sequence of nodes.3483	std::vector<regex_node> branches; // Branches of the current OR group.3484};34853486// Creates a SEQ node from the current state.3487inline regex_node make_regex_seq(regex_state &st) {3488	return regex_node("SEQ", st.cur);3489}34903491// Finishes current state.3492inline regex_node finish_regex_state(regex_state &st) {3493	// SEQ.3494	if (st.branches.empty())3495		return make_regex_seq(st);34963497	// OR.3498	st.branches.push_back(make_regex_seq(st));3499	return regex_node("OR", st.branches);3500}35013502// Parses a regex pattern into a tree, computing the number of ways to match the3503// pattern.3504inline regex_node parse_regex(std::string regex) {3505	std::string new_regex;3506	for (char c : regex)3507		if (c != ' ')3508			new_regex += c;3509	swap(regex, new_regex);3510	regex_state cur;3511	std::vector<regex_state> stack;35123513	for (size_t i = 0; i < regex.size(); ++i) {3514		char c = regex[i];35153516		if (c == '(') {3517			// Pushes the current state to the stack.3518			stack.push_back(std::move(cur));3519			cur = regex_state();3520		} else if (c == ')') {3521			// Finishes the current state, and adds it to the parent.3522			regex_node node = finish_regex_state(cur);35233524			tgen_ensure(!stack.empty(), "str: invalid regex: unmatched `)`");3525			cur = std::move(stack.back());3526			stack.pop_back();35273528			cur.cur.push_back(std::move(node));3529		} else if (c == '|') {3530			// Starts a new OR group.3531			regex_node node = make_regex_seq(cur);3532			cur.branches.push_back(std::move(node));3533		} else if (c == '[') {3534			// Parses a character class.3535			std::string chars;35363537			for (++i; i < regex.size() and regex[i] != ']'; ++i) {3538				if (i + 2 < regex.size() and regex[i + 1] == '-') {3539					char a = regex[i], b = regex[i + 2];3540					if (a > b)3541						std::swap(a, b);3542					for (char x = a; x <= b; ++x)3543						chars += x;3544					i += 2;3545				} else3546					chars += regex[i];3547			}35483549			tgen_ensure(i < regex.size() and regex[i] == ']',3550						"str: invalid regex: unmatched `[`");3551			cur.cur.emplace_back("[" + chars + "]");3552		} else if (c == '{') {3553			// Parses a repetition.3554			++i;3555			int l = -1, r = -1;35563557			while (i < regex.size() and3558				   isdigit(static_cast<unsigned char>(regex[i]))) {3559				if (l == -1)3560					l = 0;3561				tgen_ensure(l <= static_cast<int>(1e8),3562							"str: invalid regex: number too large inside `{}`");3563				l = 10 * l + (regex[i] - '0');3564				++i;3565			}35663567			if (i < regex.size() and regex[i] == ',') {3568				++i;3569				while (i < regex.size() and3570					   isdigit(static_cast<unsigned char>(regex[i]))) {3571					if (r == -1)3572						r = 0;3573					tgen_ensure(3574						r <= static_cast<int>(1e8),3575						"str: invalid regex: number too large inside `{}`");3576					r = 10 * r + (regex[i] - '0');3577					++i;3578				}3579			} else3580				r = l;35813582			tgen_ensure(i < regex.size() and regex[i] == '}',3583						"str: invalid regex: unmatched `{`");3584			tgen_ensure(l != -1 and r != -1,3585						"str: invalid regex: missing number inside `{}`");3586			tgen_ensure(l <= r,3587						"str: invalid regex: invalid range inside `{}`");35883589			// Creates a REP node from the previous node.3590			tgen_ensure(!cur.cur.empty(),3591						"str: invalid regex: expected expression before `{}`");35923593			regex_node rep(l, r, cur.cur.back());3594			cur.cur.pop_back();3595			cur.cur.push_back(std::move(rep));3596		} else {3597			// Creates a char node.3598			cur.cur.emplace_back(std::string(1, c));3599		}3600	}36013602	tgen_ensure(stack.empty(), "str: invalid regex: unmatched `(`");3603	return finish_regex_state(cur);3604}36053606// Generates a uniformly random string that matches the given regex.3607inline void gen_regex(const regex_node &node, std::string &str) {3608	// For [chars], generate a random character from the list.3609	if (node.pattern_[0] == '[') {3610		str += node.pattern_[1 + next<int>(0, node.pattern_.size() - 3)];3611		return;3612	}36133614	// For REP, generate a random number of times to repeat the pattern.3615	if (node.left_bound_ != -1) {3616		// Generates a random value W from 0 to num_ways.3617		// log(W) = log(random(0, 1) * num_ways)3618		//        = log(random(0, 1)) + log(num_ways).3619		double log_rand = math::detail::log_space(next<double>(0, 1)) +3620						  node.log_space_num_ways_;3621		double cur_prob = math::detail::LOG_ZERO;3622		double child_num_ways = node.children_[0].log_space_num_ways_;36233624		for (int i = node.left_bound_; i <= node.right_bound_; ++i) {3625			cur_prob =3626				math::detail::add_log_space(cur_prob, i * child_num_ways);3627			if (log_rand <= cur_prob) {3628				for (int j = 0; j < i; ++j)3629					gen_regex(node.children_[0], str);3630				return;3631			}3632		}36333634		tgen_ensure_against_bug(false,3635								"str: log_rand > cur_prob in REP gen_regex");3636	}36373638	// For SEQ, generate all children.3639	if (!node.children_.empty() and node.pattern_ == "SEQ") {3640		for (const regex_node &child : node.children_)3641			gen_regex(child, str);3642		return;3643	}36443645	// For OR, generate a random child.3646	if (!node.children_.empty() and node.pattern_ == "OR") {3647		// Generates a random value W from 0 to num_ways.3648		// log(W) = log(random(0, 1) * num_ways)3649		//        = log(random(0, 1)) + log(num_ways).3650		double log_rand = math::detail::log_space(next<double>(0, 1)) +3651						  node.log_space_num_ways_;3652		double cur_prob = math::detail::LOG_ZERO;36533654		for (const regex_node &child : node.children_) {3655			cur_prob = math::detail::add_log_space(cur_prob,3656												   child.log_space_num_ways_);3657			if (log_rand <= cur_prob) {3658				gen_regex(child, str);3659				return;3660			}3661		}36623663		tgen_ensure_against_bug(false,3664								"str: log_rand > cur_prob in OR gen_regex");3665	}36663667	// For char, generate the character.3668	detail::tgen_ensure_against_bug(3669		node.pattern_.size() == 1,3670		"str: invalid regex: expected single character, but got `" +3671			node.pattern_ + "`");3672	str += node.pattern_[0];3673}36743675// Formats a regex string with given arguments.3676template <typename... Args>3677std::string regex_format(const std::string &s, Args &&...args) {3678	if constexpr (sizeof...(Args) == 0) {3679		return s;3680	} else {3681		int size = std::snprintf(nullptr, 0, s.c_str(), args...) + 1;3682		std::string buf(size, '\0');3683		std::snprintf(buf.data(), size, s.c_str(), args...);3684		buf.pop_back(); // remove '\0'3685		return buf;3686	}3687}36883689} // namespace detail36903691/*3692 * String generator.3693 */36943695struct str : gen_base<str> {3696	std::optional<list<char>> list_; // List of characters.3697	std::optional<detail::regex_node>3698		root_; // Root node of the regex tree for the whole string.36993700	// Creates generator for strings of size 'size', with random characters in3701	// [value_left, value_right].3702	str(int size, char value_left = 'a', char value_right = 'z') {3703		tgen_ensure(size > 0, "str: size must be positive");3704		list_ = list<char>(size, value_left, value_right);3705	}37063707	// Creates generator for strings of size 'size', with random characters in3708	// 'chars'.3709	str(int size, std::set<char> chars) {3710		tgen_ensure(size > 0, "str: size must be positive");3711		list_ = list<char>(size, chars);3712	}37133714	// Creates generator for strings that match the given regex.3715	template <typename... Args> str(const std::string &regex, Args &&...args) {3716		tgen_ensure(regex.size() > 0, "str: regex must be non-empty");37173718		root_ = detail::parse_regex(3719			detail::regex_format(regex, std::forward<Args>(args)...));3720	}37213722	// Restricts strings for str[idx] = value.3723	str &fix(int idx, char character) {3724		tgen_ensure(!root_, "str: cannot add restriction for regex");3725		list_->fix(idx, character);3726		return *this;3727	}37283729	// Restricts strings for list[S] to be equal, for given subset S of indices.3730	str &equal(std::set<int> indices) {3731		tgen_ensure(!root_, "str: cannot add restriction for regex");3732		list_->equal(indices);3733		return *this;3734	}37353736	// Restricts strings for str[idx_1] = str[idx_2].3737	str &equal(int idx_1, int idx_2) {3738		tgen_ensure(!root_, "str: cannot add restriction for regex");3739		list_->equal(idx_1, idx_2);3740		return *this;3741	}37423743	// Restricts strings for str[left..right] to have all equal values.3744	str &equal_range(int left, int right) {3745		tgen_ensure(!root_, "str: cannot add restriction for regex");3746		list_->equal_range(left, right);3747		return *this;3748	}37493750	// Restricts strings for all equal chars.3751	str &all_equal() {3752		tgen_ensure(!root_, "str: cannot add restriction for regex");3753		list_->all_equal();3754		return *this;3755	}37563757	// Restricts strings for str[left..right] to be a palindrome.3758	str &palindrome(int left, int right) {3759		tgen_ensure(!root_, "str: cannot add restriction for regex");3760		tgen_ensure(0 <= left and left <= right and right < list_->size_,3761					"str: range indices must be valid");3762		for (int i = left; i < right - (i - left); ++i)3763			equal(i, right - (i - left));3764		return *this;3765	}37663767	// Restricts strings for the entire string to be a palindrome.3768	str &palindrome() {3769		tgen_ensure(!root_, "str: cannot add restriction for regex");3770		return palindrome(0, list_->size_ - 1);3771	}37723773	// Restricts strings for str[S] to be different (distinct), for given subset3774	// S of indices.3775	str &different(std::set<int> indices) {3776		tgen_ensure(!root_, "str: cannot add restriction for regex");3777		list_->different(indices);3778		return *this;3779	}37803781	// Restricts strings for str[idx_1] != str[idx_2].3782	str &different(int idx_1, int idx_2) {3783		tgen_ensure(!root_, "str: cannot add restriction for regex");3784		list_->different(idx_1, idx_2);3785		return *this;3786	}37873788	// Restricts lists for list[left..right] to have all different chars.3789	str &different_range(int left, int right) {3790		tgen_ensure(!root_, "str: cannot add restriction for regex");3791		list_->different_range(left, right);3792		return *this;3793	}37943795	// Restricts strings for all chars to be different.3796	str &all_different() {3797		tgen_ensure(!root_, "str: cannot add restriction for regex");3798		list_->all_different();3799		return *this;3800	}38013802	// str value.3803	struct value : gen_value_base<value> {3804		using tgen_is_sequential_tag = detail::is_sequential_tag;38053806		using value_type = char;3807		using std_type = std::string;3808		std::string str_;38093810		value(const std::string &str) : str_(str) {3811			tgen_ensure(!str_.empty(), "str: value: cannot be empty");3812		}38133814		// Fetches size.3815		int size() const { return str_.size(); }38163817		// Fetches position idx.3818		char &operator[](int idx) {3819			tgen_ensure(0 <= idx and idx < size(),3820						"str: value: index out of bounds");3821			return str_[idx];3822		}3823		const char &operator[](int idx) const {3824			tgen_ensure(0 <= idx and idx < size(),3825						"str: value: index out of bounds");3826			return str_[idx];3827		}38283829		// Sorts characters in non-decreasing order.3830		// O(n log n).3831		value &sort() {3832			std::sort(str_.begin(), str_.end());3833			return *this;3834		}38353836		// Reverses string.3837		// O(n).3838		value &reverse() {3839			std::reverse(str_.begin(), str_.end());3840			return *this;3841		}38423843		// Lowercases all characters.3844		// O(n).3845		value &lowercase() {3846			for (char &c : str_)3847				c = std::tolower(c);3848			return *this;3849		}38503851		// Uppercases all characters.3852		// O(n).3853		value &uppercase() {3854			for (char &c : str_)3855				c = std::toupper(c);3856			return *this;3857		}38583859		// Concatenates two values.3860		// Linear.3861		value operator+(const value &rhs) const {3862			return value(str_ + rhs.str_);3863		}38643865		// Shuffles string uniformly.3866		// O(n).3867		value &shuffle() {3868			for (int i = 0; i < size(); ++i)3869				std::swap(str_[i], str_[next(0, size() - 1)]);3870			return *this;3871		}38723873		// Returns a random character uniformly.3874		// O(1).3875		char pick() const { return str_[next<int>(0, size() - 1)]; }38763877		// Returns str_[i] with probability proportional to distribution[i].3878		// O(1).3879		template <typename Dist>3880		char pick_by_distribution(const std::vector<Dist> &distribution) const {3881			tgen_ensure(static_cast<size_t>(size()) == distribution.size(),3882						"value and distribution must have the same size");3883			return str_[next_by_distribution(distribution)];3884		}3885		template <typename Dist>3886		char pick_by_distribution(3887			const std::initializer_list<Dist> &distribution) const {3888			return pick_by_distribution(std::vector<Dist>(distribution));3889		}38903891		// Chooses k characters uniformly, as in a subsequence of size k.3892		// O(n).3893		value choose(int k) const {3894			tgen_ensure(0 < k and k <= size(),3895						"number of elements to choose must be valid");3896			std::string new_str;3897			int need = k;3898			for (int i = 0; need > 0; ++i) {3899				int left = size() - i;3900				if (next(1, left) <= need) {3901					new_str.push_back(str_[i]);3902					need--;3903				}3904			}3905			return value(new_str);3906		}39073908		// Prints to std::ostream.3909		friend std::ostream &operator<<(std::ostream &out, const value &val) {3910			return out << val.str_;3911		}39123913		// Gets a std::string representing the value.3914		std::string to_std() const { return std_type(str_); }3915	};39163917	// Generates str value.3918	// If created from restrictions: O(n log n).3919	// If created from regex: expected linear.3920	value gen() const {3921		if (root_) {3922			// Regex.3923			std::string ret_str;3924			gen_regex(*root_, ret_str);3925			return value(ret_str);3926		} else {3927			// List.3928			std::vector<char> vec = list_->gen().to_std();3929			return value(std::string(vec.begin(), vec.end()));3930		}3931	}3932};39333934/************3935 *          *3936 *   PAIR   *3937 *          *3938 ************/39393940namespace detail {39413942// Generates pair first == second.3943// O(1).3944template <typename T> std::pair<T, T> gen_eq(T L1, T R1, T L2, T R2) {3945	T L = std::max(L1, L2);3946	T R = std::min(R1, R2);39473948	tgen_ensure(L <= R, "pair: no valid values to generate");3949	T x = next<T>(L, R);3950	return {x, x};3951}39523953// Returns {R1-L1+1, R2-L2+1}.3954template <typename T>3955std::pair<u128, u128> get_n_and_m(T L1, T R1, T L2, T R2) {3956	u128 n = static_cast<i128>(R1) - L1 + 1;3957	u128 m = static_cast<i128>(R2) - L2 + 1;3958	return {n, m};3959}39603961// Returns first + first+1 + ... + last,3962// num_terms terms. Avoids overflow.3963static u128 pos_arith_sum(u128 first, u128 last, u128 num_terms) {3964	u128 x = first + last, y = num_terms;39653966	// x * y / 2, avoiding overflow.3967	if (x % 2 == 0)3968		x /= 2;3969	else3970		y /= 2;39713972	return x * y;3973}39743975// Generates pair first != second.3976// O(1) expected.3977template <typename T> std::pair<T, T> gen_neq(T L1, T R1, T L2, T R2) {3978	auto [n, m] = get_n_and_m(L1, R1, L2, R2);39793980	T L_intersect = std::max(L1, L2);3981	T R_intersect = std::min(R1, R2);3982	u128 inter = static_cast<i128>(R_intersect) - L_intersect + 1;39833984	u128 total = n * m - inter;3985	tgen_ensure(total > 0, "pair: no valid values to generate");39863987	// Runs O(1) expected times in the worst case.3988	T a, b;3989	do {3990		a = next<T>(L1, R1);3991		b = next<T>(L2, R2);3992	} while (a == b);39933994	return {a, b};3995}39963997// For lt, splits 'second' into two regions:3998// 1) second <= R1 -> number of 'first' is (second - L1)3999// 2) second >  R1 -> number of 'first' is (R1 - L1 + 1)4000// Returns {count_region1, count_region2}.4001// O(1).4002template <typename T>4003std::pair<u128, u128> count_lt_regions(T L1, T R1, T L2, T R2) {4004	auto [n, m] = get_n_and_m(L1, R1, L2, R2);40054006	// 'second' must be >= L1 + 1.4007	i128 L_second = std::max<i128>(L2, static_cast<i128>(L1) + 1);4008	i128 R_second = R2;40094010	// Split point for 'second'.4011	i128 split = std::min<i128>(R_second, R1);40124013	// Region 1: b in [L_second, split].4014	u128 len1 = std::max<i128>(0, split - L_second + 1);40154016	u128 count_region1 = 0;4017	if (len1 > 0) {4018		// For b in [L_second, split], there are (b - L1) ways.4019		i128 first = L_second - L1;4020		i128 last = split - L1;40214022		// Arithmetic series first + (first + 1) + ... + last, len1 terms.4023		count_region1 = pos_arith_sum(first, last, len1);4024	}40254026	// Region 2: b > R1.4027	// For b in [R1+1, R_second], there are 'n' ways.4028	i128 L_second_region2 = std::max(L_second, static_cast<i128>(R1) + 1);40294030	u128 len2 = std::max<i128>(0, R_second - L_second_region2 + 1);4031	u128 count_region2 = len2 * n;40324033	return {count_region1, count_region2};4034}40354036// Generates pair first < second.4037// O(log(R1 - L1 + 1) + log(R2 - L2 + 1)).4038template <typename T> std::pair<T, T> gen_lt(T L1, T R1, T L2, T R2) {4039	auto [n, m] = get_n_and_m(L1, R1, L2, R2);40404041	// 'second' needs to be at least L1 + 1 to have a valid value for4042	// 'first'.4043	i128 L_second = std::max<i128>(L2, static_cast<i128>(L1) + 1);4044	i128 R_second = R2;40454046	// Splits 'second' into two regions:4047	// 1) b <= R1 -> number of 'first' is (b - L1);4048	// 2) b >  R1 -> number of 'first' is (R1 - L1 + 1).4049	i128 split = std::min<i128>(R_second, R1);40504051	auto [count_region1, count_region2] = count_lt_regions(L1, R1, L2, R2);4052	u128 total = count_region1 + count_region2;4053	tgen_ensure(total > 0, "pair: no valid values to generate");40544055	u128 k = detail::next128(total);4056	if (k < count_region1) {4057		// Region 1: invert arithmetic series.40584059		// For b in [L_second, split].4060		u128 len1 = std::max<i128>(0, split - L_second + 1);40614062		// We consider b in [L_second, L_second + d].4063		// Each b contributes (b - L1) = base + (b - L_second).4064		// So we sum: base + (base+1) + ... + (base+d)4065		// d in [0, len1).40664067		i128 base = L_second - L1;4068		i128 lo = 0, hi = static_cast<i128>(len1) - 1;40694070		while (lo < hi) {4071			i128 mid = lo + (hi - lo) / 2;40724073			if (pos_arith_sum(base, base + mid, mid + 1) <= k)4074				lo = mid + 1;4075			else4076				hi = mid;4077		}4078		i128 d = lo;40794080		// Subtracts prefix sum with d-1 terms from k.4081		if (d > 0)4082			k -= pos_arith_sum(base, base + d - 1, d);40834084		return {L1 + static_cast<T>(k), L_second + d};4085	} else {4086		// Region 2: uniform block of size n.4087		k -= count_region1;40884089		// For b in [R1+1, R_second], there are 'n' ways.4090		i128 L_second_region2 = std::max(L_second, static_cast<i128>(R1) + 1);40914092		return {L1 + static_cast<T>(k % n),4093				L_second_region2 + static_cast<T>(k / n)};4094	}4095}40964097// Generates pair first > second.4098// O(log(R1 - L1 + 1) + log(R2 - L2 + 1)).4099template <typename T> std::pair<T, T> gen_gt(T L1, T R1, T L2, T R2) {4100	auto [first, second] = gen_lt(L2, R2, L1, R1);4101	return {second, first};4102}41034104// Generates pair first <= second.4105// O(log(R1 - L1 + 1) + log(R2 - L2 + 1)).4106template <typename T> std::pair<T, T> gen_leq(T L1, T R1, T L2, T R2) {4107	// Counts how many pairs are there with first = second.4108	i128 L_intersect = std::max(L1, L2);4109	i128 R_intersect = std::min(R1, R2);4110	u128 eq_count = std::max<i128>(0, R_intersect - L_intersect + 1);41114112	// Counts how many pairs are there with first < second.4113	auto [lt_region1, lt_region2] = count_lt_regions(L1, R1, L2, R2);4114	u128 lt_count = lt_region1 + lt_region2;41154116	u128 total = eq_count + lt_count;4117	tgen_ensure(total > 0, "pair: no valid values to generate");41184119	if (detail::next128(total) < eq_count)4120		return gen_eq(L1, R1, L2, R2);4121	return gen_lt(L1, R1, L2, R2);4122}41234124// Generates pair first >= second.4125// O(log(R1 - L1 + 1) + log(R2 - L2 + 1)).4126template <typename T> std::pair<T, T> gen_geq(T L1, T R1, T L2, T R2) {4127	auto [first, second] = gen_leq(L2, R2, L1, R1);4128	return {second, first};4129}41304131}; // namespace detail41324133/*4134 * Pair generator.4135 *4136 * Pairs of integral types.4137 */41384139template <typename T> struct pair : gen_base<pair<T>> {4140	std::pair<T, T> first_, second_; // Range of first and second values.4141	// Type of restriction.4142	enum class restriction_type { eq, neq, lt, gt, leq, geq, unspecified };4143	restriction_type type_ = restriction_type::unspecified;41444145	// Creates a pair with random values in [first_l, first_r] and [second_l,4146	// second_r].4147	pair(T first_left, T first_right, T second_left, T second_right)4148		: first_(first_left, first_right), second_(second_left, second_right) {4149		tgen_ensure(first_left <= first_right,4150					"pair: first range must be valid");4151		tgen_ensure(second_left <= second_right,4152					"pair: second range must be valid");4153	}41544155	// Creates a pair with random values in [both_l, both_r].4156	pair(T both_left, T both_right)4157		: pair(both_left, both_right, both_left, both_right) {}41584159	// Restricts pair for first = second.4160	pair &eq() {4161		type_ = restriction_type::eq;4162		return *this;4163	}41644165	// Restricts pair for first != second.4166	pair &neq() {4167		type_ = restriction_type::neq;4168		return *this;4169	}41704171	// Restricts pair for first < second.4172	pair &lt() {4173		type_ = restriction_type::lt;4174		return *this;4175	}41764177	// Restricts pair for first > second.4178	pair &gt() {4179		type_ = restriction_type::gt;4180		return *this;4181	}41824183	// Restricts pair for first <= second.4184	pair &leq() {4185		type_ = restriction_type::leq;4186		return *this;4187	}41884189	// Restricts pair for first >= second.4190	pair &geq() {4191		type_ = restriction_type::geq;4192		return *this;4193	}41944195	// Pair value.4196	struct value : gen_value_base<value> {4197		using value_type = T;4198		using std_type = std::pair<T, T>;41994200		std::pair<T, T> pair_;4201		char sep_;42024203		value(const std::pair<T, T> &pair) : pair_(pair), sep_(' ') {}4204		value(const T &first, const T &second)4205			: pair_(first, second), sep_(' ') {}42064207		T first() const { return pair_.first; }4208		T second() const { return pair_.second; }42094210		// Sets the separator for the pair, for printing.4211		value &separator(char sep) {4212			sep_ = sep;4213			return *this;4214		}42154216		// Prints to std::ostream, separated by sep_.4217		friend std::ostream &operator<<(std::ostream &out, const value &val) {4218			return out << val.pair_.first << val.sep_ << val.pair_.second;4219		}42204221		// Gets a std::pair representing the value.4222		auto to_std() const {4223			if constexpr (!detail::is_generator_value<T>::value) {4224				return pair_;4225			} else {4226				std::pair<typename T::std_type, typename T::std_type> pair(4227					pair_.first.to_std(), pair_.second.to_std());4228				return pair;4229			}4230		}4231	};42324233	// Generates a random pair.4234	// O(log(R1 - L1 + 1) + log(R2 - L2 + 1)).4235	value gen() const {4236		T L1 = first_.first, R1 = first_.second;4237		T L2 = second_.first, R2 = second_.second;42384239		switch (type_) {4240		case restriction_type::unspecified:4241			return {next<T>(L1, R1), next<T>(L2, R2)};4242		case restriction_type::eq:4243			return detail::gen_eq<T>(L1, R1, L2, R2);4244		case restriction_type::neq:4245			return detail::gen_neq<T>(L1, R1, L2, R2);4246		case restriction_type::lt:4247			return detail::gen_lt<T>(L1, R1, L2, R2);4248		case restriction_type::gt:4249			return detail::gen_gt<T>(L1, R1, L2, R2);4250		case restriction_type::leq:4251			return detail::gen_leq<T>(L1, R1, L2, R2);4252		case restriction_type::geq:4253			return detail::gen_geq<T>(L1, R1, L2, R2);4254		}4255		throw detail::error("pair: unknown restriction type");4256	}4257};42584259/************4260 *          *4261 *   TREE   *4262 *          *4263 ************/42644265namespace detail {42664267// Generates edges from Prufer sequence.4268// O(n).4269inline std::vector<std::pair<int, int>> edges_from_prufer(std::vector<int> p) {4270	int n = p.size() + 2;42714272	// Degrees.4273	std::vector<int> d(n, 1);4274	for (int i : p)4275		d[i]++;42764277	// Adds last vertex.4278	p.push_back(n - 1);42794280	// Finds first vertex with degree 1.4281	int idx, u;4282	idx = u = find(d.begin(), d.end(), 1) - d.begin();42834284	// Generates edges.4285	std::vector<std::pair<int, int>> edges;4286	for (int v : p) {4287		edges.emplace_back(u, v);4288		if (--d[v] == 1 and v < idx)4289			u = v;4290		else4291			idx = u = find(d.begin() + idx + 1, d.end(), 1) - d.begin();4292	}4293	return edges;4294}42954296// Disjoint set union (union-find) for connectivity queries.4297struct dsu {4298	std::vector<int> parent_;4299	std::vector<unsigned char> rank_;43004301	// Creates a dsu with `n` elements, indexed from 0 to n-1.4302	// Initially every element is in its own set.4303	// O(n).4304	dsu(int n) : parent_(n), rank_(n, 0) {4305		for (int i = 0; i < n; ++i)4306			parent_[i] = i;4307	}43084309	// Adds new elements to the dsu, each in their own new set.4310	// O(k) amortized.4311	void add_elements(int k) {4312		for (int i = 0; i < k; ++i) {4313			int new_id = parent_.size();4314			parent_.push_back(new_id);4315			rank_.push_back(0);4316		}4317	}43184319	// Finds representative of set containing i.4320	// O(alpha(n)) amortized, O(log n) worst case.4321	int find(int i) {4322		return parent_[i] == i ? i : parent_[i] = find(parent_[i]);4323	}43244325	// Merges components of `a` and `b`. Returns if the sets were united, and4326	// false if a and b were in the same set.4327	// O(alpha(n)) amortized, O(log n) worst case.4328	bool unite(int a, int b) {4329		a = find(a);4330		b = find(b);4331		if (a == b)4332			return false;4333		if (rank_[a] > rank_[b])4334			std::swap(a, b);4335		parent_[a] = b;4336		if (rank_[a] == rank_[b])4337			++rank_[b];4338		return true;4339	}4340};43414342} // namespace detail43434344// Forward declaration of wgraph.4345template <typename VWeight, typename EWeight> struct wgraph;43464347/*4348 * Tree generator.4349 *4350 * Unrooted trees with `n` vertices, indexed from 0 to n-1.4351 * These are unrooted undirected labeled trees, that is, isomorphism is not4352 * taken into account. VWeight is the type of vertex weights, and EWeight is4353 * the type of edge weights. Generator does not generate weights. The weights4354 * are to be set in the wtree::value.4355 */43564357template <typename VWeight, typename EWeight>4358struct wtree : gen_base<wtree<VWeight, EWeight>> {4359	int n_;								  // Number of vertices.4360	std::set<std::pair<int, int>> edges_; // Edges that were set.43614362	// Creates tree generator with `n` vertices.4363	// O(1).4364	wtree(int n) : n_(n) {4365		tgen_ensure(n > 0, "wtree: number of vertices must be positive");4366	}43674368	// Adds edge between u and v (this edge must be generated).4369	// O(log n).4370	wtree &add_edge(int u, int v) {4371		tgen_ensure(0 <= std::min(u, v) and std::max(u, v) < n_,4372					"wtree: vertices must be indexed in [0, n)");4373		tgen_ensure(u != v, "wtree: cannot add self loop to tree");43744375		if (u > v)4376			std::swap(u, v);4377		edges_.emplace(u, v);4378		return *this;4379	}43804381	// Tree value.4382	//4383	// Edges are stored in both directions in adjacency list, but only u < v in4384	// edge list.4385	struct value : gen_value_base<value> {4386		using std_type = std::pair<int, std::vector<std::set<int>>>;43874388		int n_;									 // Number of vertices.4389		std::vector<std::set<int>> adj_;		 // Adjacency list.4390		std::vector<std::pair<int, int>> edges_; // Edge list.4391		bool add_1_;   // If should add 1 for printing vertex ids.4392		bool print_n_; // If should print n.4393		std::optional<int> print_parents_; // If should print in parent style4394										   // (stores the root).4395		std::optional<std::vector<VWeight>> vertex_weights_; // Vertex weights.4396		std::optional<std::vector<EWeight>>4397			edge_weights_; // Edge weights (in same order as edges_).4398		detail::dsu dsu_;  // Connectivity of current edges (for cycle checks).43994400		// Creates value from adjacency list.4401		// O(n).4402		value(const std::vector<std::set<int>> &adj)4403			: n_(static_cast<int>(adj.size())), adj_(adj), add_1_(false),4404			  print_n_(false), dsu_(n_) {4405			for (int u = 0; u < n_; ++u)4406				for (auto v : adj[u]) {4407					tgen_ensure(4408						0 <= v and v < n_,4409						"wtree: value: vertices must be indexed in [0, n)");4410					// Symmetric adjacency: count each undirected edge once.4411					if (u < v) {4412						edges_.emplace_back(u, v);4413						tgen_ensure(4414							dsu_.unite(u, v),4415							"wtree: value: initial graph must form a tree");4416					}4417				}4418		}44194420		// Creates value from `n` and edge list.4421		// O(n).4422		value(int n, const std::vector<std::pair<int, int>> &edges)4423			: n_(n), adj_(n), add_1_(false), print_n_(false), dsu_(n) {4424			edges_.reserve(edges.size());4425			for (auto [u, v] : edges) {4426				tgen_ensure(0 <= std::min(u, v) and std::max(u, v) < n,4427							"wtree: value: vertices must be indexed in [0, n)");4428				tgen_ensure(dsu_.unite(u, v),4429							"wtree: value: initial graph must form a tree");4430				if (u > v)4431					std::swap(u, v);4432				edges_.emplace_back(u, v);4433				adj_[u].insert(v);4434				adj_[v].insert(u);4435			}4436		}4437		value(int n, const std::set<std::pair<int, int>> &edges)4438			: value(n, std::vector<std::pair<int, int>>(edges.begin(),4439														edges.end())) {}4440		value(int n, const std::initializer_list<std::pair<int, int>> &edges)4441			: value(n, std::vector<std::pair<int, int>>(edges)) {}44424443		// Creates tree from graph via Kruskal-like random spanning tree.4444		// Implemented after wgraph definition.4445		// O(n + m alpha(n)).4446		value(const typename wgraph<VWeight, EWeight>::value &g);44474448		// Weight type conversion.4449		// O(n).4450		template <typename NewVWeight, typename NewEWeight>4451		typename wtree<NewVWeight, NewEWeight>::value4452		convert_weight_types() const {4453			tgen_ensure(!vertex_weights_.has_value() and4454							!edge_weights_.has_value(),4455						"wtree: value: cannot convert weight type after "4456						"assigning weights");44574458			typename wtree<NewVWeight, NewEWeight>::value new_tree(adj_);4459			new_tree.add_1_ = add_1_;4460			new_tree.print_n_ = print_n_;4461			new_tree.print_parents_ = print_parents_;4462			return new_tree;4463		}44644465		// Fetches number of vertices.4466		int n() const { return n_; }44674468		// Fetches a const ref. to adjacency list.4469		const std::vector<std::set<int>> &adj() const { return adj_; }44704471		// Fetches a const ref. to edge list.4472		const std::vector<std::pair<int, int>> &edges() const { return edges_; }44734474		// Fetches a const ref. to vertex weights.4475		const std::optional<std::vector<VWeight>> &vertex_weights() const {4476			return vertex_weights_;4477		}44784479		// Fetches a const ref. to edge weights.4480		const std::optional<std::vector<EWeight>> &edge_weights() const {4481			return edge_weights_;4482		}44834484		// Sets vertex weights.4485		// O(n).4486		template <typename NewVWeight = VWeight>4487		typename wtree<NewVWeight, EWeight>::value set_vertex_weights(4488			const std::vector<NewVWeight> &vertex_weights) const {4489			tgen_ensure(static_cast<int>(vertex_weights.size()) == n(),4490						"wtree: value: must give `n` vertex weights");44914492			auto new_tree = convert_weight_types<NewVWeight, EWeight>();4493			new_tree.vertex_weights_ = vertex_weights;4494			return new_tree;4495		}44964497		// Sets edge weights.4498		// O(n).4499		template <typename NewEWeight = EWeight>4500		typename wtree<VWeight, NewEWeight>::value4501		set_edge_weights(const std::vector<NewEWeight> &edge_weights) const {4502			tgen_ensure(4503				edge_weights.size() == edges().size(),4504				"wtree: value: must give `edges().size()` edge weights");45054506			auto new_tree = convert_weight_types<VWeight, NewEWeight>();4507			new_tree.edge_weights_ = edge_weights;4508			return new_tree;4509		}45104511		// Enables edge-weighted mode before adding weighted edges4512		// incrementally. The tree must have no edges yet. O(1).4513		value &edge_weighted() {4514			tgen_ensure(edges().size() == 0,4515						"wtree: value: edge_weighted requires a tree with no "4516						"edges");4517			tgen_ensure(!edge_weights_.has_value(),4518						"wtree: value: tree is already edge-weighted");45194520			edge_weights_ = std::vector<EWeight>();4521			return *this;4522		}45234524		// Adds 1 to vertex ids, for printing.4525		// O(1).4526		value &add_1() {4527			add_1_ = true;4528			return *this;4529		}45304531		// Prints `n` on a new line before printing the tree.4532		// O(1).4533		value &print_n() {4534			print_n_ = true;4535			return *this;4536		}45374538		// Prints the tree in parent style.4539		// If root = -1, the root is considered to be 0, and its parent is not4540		// printed. Otherwise, prints the parent of the root as -1. If root = n,4541		// randomizes the root. O(1).4542		value &print_parents(int root = -1) {4543			tgen_ensure(root == -1 or (0 <= root and root < n()) or root == n(),4544						"wtree: value: root must be -1, `n`, or in [0, n)");4545			print_parents_ = root;4546			return *this;4547		}45484549		// Shuffles the tree's vertex labels (except those in `indices`,4550		// which keep their current label) and edge order. The change is4551		// applied eagerly to the underlying adjacency list, edge list,4552		// vertex weights and edge weights.4553		// O(n).4554		value &shuffle_except(std::set<int> indices) {4555			// Builds the relabeling: for each vertex `i`, `new_label[i]` is4556			// its new id. Vertices in `indices` keep their label; the others4557			// are permuted among themselves.4558			std::vector<int> new_label(n());4559			std::vector<int> shuffled;4560			for (int i = 0; i < n(); ++i) {4561				if (indices.count(i))4562					new_label[i] = i;4563				else4564					shuffled.push_back(i);4565			}4566			std::vector<int> targets = shuffled;4567			tgen::shuffle(targets.begin(), targets.end());4568			for (size_t k = 0; k < shuffled.size(); ++k)4569				new_label[shuffled[k]] = targets[k];45704571			// Rewrites adjacency list with new labels.4572			std::vector<std::set<int>> new_adj(n());4573			for (int u = 0; u < n(); ++u)4574				for (int v : adj_[u])4575					new_adj[new_label[u]].insert(new_label[v]);4576			adj_ = std::move(new_adj);45774578			// Rewrites edges with new labels (canonical undirected order).4579			for (auto &[u, v] : edges_) {4580				u = new_label[u];4581				v = new_label[v];4582				if (u > v)4583					std::swap(u, v);4584			}45854586			// Permutes vertex weights to match the new labels.4587			if (vertex_weights_.has_value()) {4588				std::vector<VWeight> new_vw(n());4589				for (int i = 0; i < n(); ++i)4590					new_vw[new_label[i]] = (*vertex_weights_)[i];4591				vertex_weights_ = std::move(new_vw);4592			}45934594			// Rebuilds the dsu so future `add_edge` calls see the new labels.4595			dsu_ = detail::dsu(n());4596			for (auto [u, v] : edges_)4597				dsu_.unite(u, v);45984599			// Shuffles edge order, keeping edge weights aligned.46004601			std::vector<int> perm(edges_.size());4602			std::iota(perm.begin(), perm.end(), 0);4603			tgen::shuffle(perm.begin(), perm.end());46044605			std::vector<std::pair<int, int>> new_edges;4606			std::optional<std::vector<EWeight>> new_ew;4607			if (edge_weights_.has_value())4608				new_ew = std::vector<EWeight>();4609			for (int i : perm) {4610				new_edges.push_back(edges_[i]);4611				if (new_ew.has_value())4612					new_ew->push_back((*edge_weights_)[i]);4613			}4614			edges_ = new_edges;4615			if (new_ew.has_value())4616				edge_weights_ = new_ew;46174618			return *this;4619		}46204621		// Shuffles the tree's vertices and edge order.4622		// O(n).4623		value &shuffle() { return shuffle_except({}); }46244625		// Adds edge (u, v).4626		// O(log n) amortized.4627		value &add_edge(int u, int v, std::optional<EWeight> w = std::nullopt) {4628			tgen_ensure(0 <= std::min(u, v) and std::max(u, v) < n(),4629						"wtree: value: vertex ids must be valid");46304631			if (u > v)4632				std::swap(u, v);46334634			if (adj_[u].count(v))4635				return *this;46364637			adj_[u].insert(v);4638			adj_[v].insert(u);4639			edges_.emplace_back(u, v);4640			tgen_ensure(dsu_.unite(u, v),4641						"wtree: value: added edge must not create a cycle");46424643			if (w.has_value()) {4644				tgen_ensure(edge_weights().has_value(),4645							"wtree: value: cannot add weighted edge to "4646							"edge-unweighted tree");46474648				edge_weights_->push_back(*w);4649			} else4650				tgen_ensure(!edge_weights().has_value(),4651							"wtree: value: cannot add unweighted edge to "4652							"edge-weighted tree");46534654			return *this;4655		}46564657		// Links tree with another `rhs`, adding the edge between u (in left4658		// tree) and v (in right tree). Ids for added vertices are updated4659		// accordingly.4660		// O(rhs.n + rhs.m * log n) amortized.4661		value &link(const value &rhs, int new_u, int new_v,4662					std::optional<EWeight> new_w = std::nullopt) {4663			tgen_ensure(0 <= new_u and new_u < n() and 0 <= new_v and4664							new_v < rhs.n(),4665						"wtree: value: vertex ids must be valid");46664667			// Edges from right-hand side.4668			int shift = n();4669			add_vertices(rhs.n(), rhs.vertex_weights());4670			for (int i = 0; i < static_cast<int>(rhs.edges().size()); ++i) {4671				auto [u, v] = rhs.edges()[i];4672				add_edge(shift + u, shift + v,4673						 rhs.edge_weights().has_value()4674							 ? std::optional<EWeight>((*rhs.edge_weights())[i])4675							 : std::nullopt);4676			}46774678			// New edge.4679			add_edge(new_u, shift + new_v, new_w);46804681			return *this;4682		}46834684		// Glues the tree with another `rhs` such that index_pairs[i].first is4685		// considered to be the same as index_pairs[i].second. Ids for added4686		// vertices are updated accordingly.4687		// O(rhs.n + rhs.m * log n) amortized.4688		value &glue(const value &rhs,4689					std::set<std::pair<int, int>> index_pairs) {4690			// Checks validity of indices.4691			std::set<int> idx_left, idx_right;4692			std::vector<int> right_id_to_left(rhs.n(), -1);4693			for (auto [l, r] : index_pairs) {4694				tgen_ensure(4695					0 <= l and l < n() and 0 <= r and r < rhs.n(),4696					"wtree: value: vertex indices to glue must be valid");4697				tgen_ensure(idx_left.count(l) == 0 and idx_right.count(r) == 0,4698							"wtree: value: must not have repeated indices "4699							"on the same side to glue");47004701				idx_left.insert(l);4702				idx_right.insert(r);4703				right_id_to_left[r] = l;4704			}47054706			// Computes new ids of right vertices.4707			std::vector<int> new_right_id(rhs.n(), -1);4708			int intersection_lt = 0;4709			std::optional<std::vector<VWeight>> rhs_vertex_weights;4710			for (int i = 0; i < rhs.n(); ++i) {4711				if (right_id_to_left[i] != -1) {4712					// Is in intersection.4713					++intersection_lt;4714					new_right_id[i] = right_id_to_left[i];4715				} else {4716					// New id.4717					new_right_id[i] = n() + i - intersection_lt;4718					if (rhs.vertex_weights().has_value()) {4719						if (!rhs_vertex_weights.has_value())4720							rhs_vertex_weights = std::vector<VWeight>();4721						rhs_vertex_weights->push_back(4722							(*rhs.vertex_weights())[i]);4723					}4724				}4725			}47264727			// Adds new vertices and edges.4728			add_vertices(rhs.n() - intersection_lt, rhs_vertex_weights);4729			for (int i = 0; i < static_cast<int>(rhs.edges().size()); ++i) {4730				auto [u, v] = rhs.edges()[i];4731				add_edge(new_right_id[u], new_right_id[v],4732						 rhs.edge_weights().has_value()4733							 ? std::optional<EWeight>((*rhs.edge_weights())[i])4734							 : std::nullopt);4735			}47364737			return *this;4738		}4739		value &glue(const value &rhs,4740					std::initializer_list<std::pair<int, int>> il) {4741			return glue(rhs, std::set<std::pair<int, int>>(il));4742		}47434744		// Glues the tree with another `rhs` at `indices`. That is, idx in4745		// `indices` are considered to be the same vertex. Ids for added4746		// vertices are updated accordingly.4747		// O(rhs.n).4748		value &glue(const value &rhs, std::set<int> indices) {4749			std::set<std::pair<int, int>> index_pairs;4750			for (auto i : indices)4751				index_pairs.emplace(i, i);4752			return glue(rhs, index_pairs);4753		}4754		value &glue(const value &rhs, const std::initializer_list<int> &il) {4755			return glue(rhs, std::set<int>(il));4756		}47574758		// Prints to std::ostream.4759		// O(n).4760		friend std::ostream &operator<<(std::ostream &out, const value &val) {4761			if (val.print_n_)4762				out << val.n() << '\n';47634764			// Prints vertex weights.4765			if (val.vertex_weights()) {4766				for (int i = 0; i < val.n(); ++i) {4767					if (i > 0)4768						out << " ";4769					out << (*val.vertex_weights())[i];4770				}4771				out << '\n';4772			}47734774			tgen_ensure(static_cast<int>(val.edges().size()) == val.n() - 1,4775						"wtree: value: invalid tree to print (number of edges "4776						"must be `n` - 1)");47774778			// Prints in parent style.4779			if (val.print_parents_.has_value()) {4780				tgen_ensure(!val.edge_weights().has_value(),4781							"wtree: value: cannot print parent style if edges "4782							"are weighted");47834784				int root = *val.print_parents_;4785				bool skip_parent_0 = root == -1;4786				if (root == -1)4787					root = 0;4788				if (root == val.n())4789					root = next(0, val.n() - 1);47904791				std::vector<int> parent(val.n(), -1);47924793				std::queue<int> q;4794				std::vector<int> vis(val.n(), false);4795				q.push(root);4796				vis[root] = true;47974798				while (q.size()) {4799					int u = q.front();4800					q.pop();4801					for (int v : val.adj()[u])4802						if (!vis[v]) {4803							vis[v] = true;4804							q.push(v);4805							parent[v] = u;4806						}4807				}48084809				if (skip_parent_0) {4810					for (int i = 1; i < val.n(); ++i) {4811						tgen_ensure(4812							parent[i] < i,4813							"wtree: value: parent of i must be less than i for "4814							"printing in parent style if root is -1");48154816						if (i > 1)4817							out << " ";4818						out << parent[i] + val.add_1_;4819					}4820				} else {4821					for (int i = 0; i < val.n(); ++i) {4822						if (i > 0)4823							out << " ";4824						out << (parent[i] == -1 ? -1 : parent[i]) + val.add_1_;4825					}4826				}48274828				out << '\n';4829				return out;4830			}48314832			// Prints edges.4833			for (int i = 0; i < static_cast<int>(val.edges().size()); ++i) {4834				auto [u, v] = val.edges()[i];4835				out << (u + val.add_1_) << " " << (v + val.add_1_);48364837				// Edge weight.4838				if (val.edge_weights().has_value())4839					out << " " << (*val.edge_weights())[i];48404841				out << '\n';4842			}48434844			return out;4845		}48464847		// Gets a std::pair<n, adj> representing the value.4848		std::pair<int, std::vector<std::set<int>>> to_std() const {4849			return std_type(n_, adj_);4850		}48514852	  private:4853		// Adds `k` vertices to the tree (labeled n, n+1, ...n+k-1). Updates4854		// `n` accordingly. This makes the tree invalid (not a tree anymore).4855		// O(k) amortized.4856		value &add_vertices(int k, std::optional<std::vector<VWeight>>4857									   new_vertex_weights = std::nullopt) {4858			n_ += k;4859			adj_.resize(n());4860			if (new_vertex_weights.has_value()) {4861				tgen_ensure(vertex_weights().has_value(),4862							"wtree: value: cannot add weighted vertices to "4863							"vertex-unweighted tree");4864				tgen_ensure(4865					static_cast<int>(new_vertex_weights->size()) == k,4866					"wtree: value: number of vertex weights must be equal "4867					"to number of added vertices");48684869				vertex_weights_->insert(vertex_weights_->end(),4870										new_vertex_weights->begin(),4871										new_vertex_weights->end());4872			} else4873				tgen_ensure(!vertex_weights().has_value(),4874							"wtree: value: cannot add unweighted vertices to "4875							"vertex-weighted tree");48764877			dsu_.add_elements(k);48784879			return *this;4880		}4881	};48824883	// Generates tree value.4884	// O(n).4885	value gen() const {4886		// Constructs adjacency list.4887		std::vector<std::vector<int>> adj(n_);4888		for (auto [u, v] : edges_) {4889			adj[u].push_back(v);4890			adj[v].push_back(u);4891		}48924893		std::vector<int> comp_size;4894		std::vector<std::vector<int>> component_ids;4895		std::vector<bool> vis(n_, false);4896		std::queue<int> q;48974898		for (int i = 0; i < n_; ++i) {4899			if (vis[i])4900				continue;49014902			vis[i] = true;4903			q.push(i);4904			comp_size.push_back(0);4905			component_ids.emplace_back();4906			while (q.size()) {4907				int u = q.front();4908				q.pop();4909				++comp_size.back();4910				component_ids.back().push_back(u);4911				for (int v : adj[u]) {4912					if (!vis[v]) {4913						vis[v] = true;4914						q.push(v);4915					}4916				}4917			}4918		}49194920		// Creates edges connecting the connected components by treating them as4921		// vertices.4922		std::vector<std::pair<int, int>> new_edges(edges_.begin(),4923												   edges_.end());4924		if (comp_size.size() > 1) {4925			std::vector<int> prufer_values =4926				many_by_distribution(comp_size.size() - 2, comp_size);4927			for (auto [u, v] : detail::edges_from_prufer(prufer_values))4928				new_edges.emplace_back(pick(component_ids[u]),4929									   pick(component_ids[v]));4930		}49314932		return value(n_, new_edges);4933	}49344935	// Generates a (not uniformly) random skewed tree.4936	// Vertex 0 is the root. For each i in 1 .. n-1, parent(i) is4937	// wnext(i, elongation), i.e. a value in [0, i) with skew controlled by4938	// elongation (see wnext).4939	// If elongation is small enough, generates a star (center 0).4940	// If elongation is large enough, generates a path (endpoints 0 and n-1).4941	// O(n).4942	static value gen_skewed(int n, int elongation) {4943		std::vector<std::pair<int, int>> edges;4944		for (int i = 1; i < n; ++i)4945			edges.emplace_back(i, wnext<int>(i, elongation));4946		return value(n, edges);4947	}49484949	// Kruskal-like random tree: random vertex pairs until connected.4950	// Not uniformly random.4951	// O(n log(n) alpha(n)) expected.4952	static value gen_kruskal(int n) {4953		tgen_ensure(n > 0, "wtree: gen_kruskal: n must be positive");4954		if (n == 1)4955			return value(1, {});49564957		detail::dsu components(n);4958		std::vector<std::pair<int, int>> edges;4959		edges.reserve(n - 1);4960		while (edges.size() < size_t(n - 1)) {4961			int u = next(0, n - 1);4962			int v = next(0, n - 1);4963			if (u == v)4964				continue;4965			if (components.unite(u, v))4966				edges.emplace_back(u, v);4967		}4968		return value(n, edges);4969	}4970};49714972/*4973 * Other types of weighted-ness.4974 */49754976// Vertex weighted tree.4977template <typename VWeight> using vtree = wtree<VWeight, int>;49784979// Edge weighted tree.4980template <typename EWeight> using etree = wtree<int, EWeight>;49814982// Unweighted tree.4983using tree = wtree<int, int>;49844985/*************4986 *           *4987 *   GRAPH   *4988 *           *4989 *************/49904991namespace detail {49924993// Canonical undirected edge key for duplicate detection; stores (min(u, v),4994// max(u, v)). O(1).4995inline uint64_t undirected_edge_key(int u, int v) {4996	if (u > v)4997		std::swap(u, v);4998	return (static_cast<uint64_t>(u) << 32) |4999		   static_cast<uint64_t>(static_cast<uint32_t>(v));5000}50015002// Directed edge key for duplicate detection; stores (u, v).5003// O(1).5004inline uint64_t directed_edge_key(int u, int v) {5005	return (static_cast<uint64_t>(u) << 32) |5006		   static_cast<uint64_t>(static_cast<uint32_t>(v));5007}50085009// Maximum number of edges in a simple graph on n vertices.5010// O(1).5011inline long long max_graph_edges(int n, bool directed, bool self_loops) {5012	if (n <= 0)5013		return 0;5014	if (directed)5015		return self_loops ? static_cast<long long>(n) * n5016						  : static_cast<long long>(n) * (n - 1);5017	return self_loops ? static_cast<long long>(n) * (n + 1) / 25018					  : static_cast<long long>(n) * (n - 1) / 2;5019}50205021// Uniform random edge for rejection sampling.5022// O(1) expected.5023inline std::pair<int, int> get_random_graph_edge(int n, bool directed,5024												 bool self_loops) {5025	if (directed) {5026		if (self_loops)5027			return {next<int>(0, n - 1), next<int>(0, n - 1)};5028		int u = next<int>(0, n - 1);5029		int v = next<int>(0, n - 1);5030		while (u == v)5031			v = next<int>(0, n - 1);5032		return {u, v};5033	}5034	if (self_loops) {5035		int u = next<int>(0, n - 1);5036		int v = next<int>(0, n - 1);5037		if (u > v)5038			std::swap(u, v);5039		return {u, v};5040	}5041	int u = next<int>(0, n - 1);5042	int v = next<int>(0, n - 1);5043	while (u == v)5044		v = next<int>(0, n - 1);5045	if (u > v)5046		std::swap(u, v);5047	return {u, v};5048}50495050// Decodes a linear edge index to (u, v) for an undirected simple graph,5051// with u < v.5052// O(log n).5053inline std::pair<int, int> decode_undirected_simple_edge(int n, long long idx) {5054	auto base = [&](int u) -> long long {5055		return static_cast<long long>(u) * (n - 1) -5056			   static_cast<long long>(u) * (u - 1) / 2;5057	};5058	int lo = 0, hi = n - 2;5059	while (lo < hi) {5060		int mid = (lo + hi + 1) / 2;5061		if (base(mid) <= idx)5062			lo = mid;5063		else5064			hi = mid - 1;5065	}5066	return {lo, lo + 1 + int(idx - base(lo))};5067}50685069// Decodes a linear edge index to (u, v) for an undirected graph with loops,5070// with u <= v.5071// O(log n).5072inline std::pair<int, int> decode_undirected_loops_edge(int n, long long idx) {5073	auto base = [&](int u) -> long long {5074		return static_cast<long long>(u) * n -5075			   static_cast<long long>(u) * (u - 1) / 2;5076	};5077	int lo = 0, hi = n - 1;5078	while (lo < hi) {5079		int mid = (lo + hi + 1) / 2;5080		if (base(mid) <= idx)5081			lo = mid;5082		else5083			hi = mid - 1;5084	}5085	return {lo, lo + int(idx - base(lo))};5086}50875088// Decodes a linear edge index to (u, v) for a directed simple graph (no loops).5089// O(1).5090inline std::pair<int, int> decode_directed_simple_edge(int n, long long idx) {5091	int u = idx / (n - 1);5092	int rem = idx % (n - 1);5093	return {u, rem + (rem >= u)};5094}50955096// Decodes a linear edge index according to graph mode.5097// O(log n) for undirected, O(1) for directed.5098inline std::pair<int, int>5099decode_graph_edge_index(int n, long long idx, bool directed, bool self_loops) {5100	if (directed) {5101		if (self_loops)5102			return {int(idx / n), int(idx % n)};5103		return decode_directed_simple_edge(n, idx);5104	}5105	if (self_loops)5106		return decode_undirected_loops_edge(n, idx);5107	return decode_undirected_simple_edge(n, idx);5108}51095110} // namespace detail51115112/*5113 * Graph generator.5114 *5115 * Graphs of `n` vertices labeled from 0 to n-1 and `m` edges.5116 * These are labeled graphs, that is, isomorphism is not taken into5117 * account. VWeight is the type of vertex weights, and EWeight is the type of5118 * edge weights. Generator does not generate weights. The weights are to be set5119 * in the wgraph::value.5120 */51215122template <typename VWeight, typename EWeight>5123struct wgraph : gen_base<wgraph<VWeight, EWeight>> {5124	int n_, m_;							  // Number of vertices and edges.5125	std::set<std::pair<int, int>> edges_; // Edges that were set.5126	bool is_directed_;					  // If graph is directed.5127	bool has_self_loops_;				  // If self-loops are allowed.51285129	// Creates graph generator with `n` vertices and `m` edges.5130	// Additionally, you can set if the graph is directed and if self loops are5131	// allowed.5132	// O(1).5133	wgraph(int n, int m, bool is_directed = false, bool has_self_loops = false)5134		: n_(n), m_(m), is_directed_(is_directed),5135		  has_self_loops_(has_self_loops) {5136		tgen_ensure(n > 0, "wgraph: number of vertices must be positive");5137	}51385139	// Adds edge between u and v (this edge must be generated).5140	// O(log m).5141	wgraph &add_edge(int u, int v) {5142		tgen_ensure(0 <= std::min(u, v) and std::max(u, v) < n_,5143					"wgraph: vertices must be indexed in [0, n)");51445145		if (!is_directed_ and u > v)5146			std::swap(u, v);5147		edges_.emplace(u, v);5148		tgen_ensure(static_cast<int>(edges_.size()) <= m_,5149					"wgraph: too many edges were added");5150		return *this;5151	}51525153	// Graph value.5154	//5155	// Edges are stored in both directions (if undirected) in adjacency list,5156	// but only u < v in edge list.5157	// Optimized for performance (lazy adjacency list; edge-list constructor5158	// stores edges only).5159	struct value : gen_value_base<value> {5160		using std_type = std::tuple<int, int, std::vector<std::set<int>>>;51615162		int n_;									 // Number of vertices.5163		std::vector<std::set<int>> adj_;		 // Adjacency list.5164		std::vector<std::pair<int, int>> edges_; // Edge list.5165		bool is_directed_;						 // If graph is directed.5166		bool add_1_;	// If should add 1 for printing vertex ids.5167		bool print_nm_; // If should print n and m.5168		mutable bool adj_built_{5169			false}; // Lazy cache: true once adj_ is built from edges_; mutable5170					// so const adj() can populate it.5171		std::optional<std::vector<VWeight>> vertex_weights_; // Vertex weights.5172		std::optional<std::vector<EWeight>>5173			edge_weights_; // Edge weights (in same order as edges_ ).51745175		// Creates value from adjacency list. The edges5176		// are considered to be directed.5177		// O(n + m).5178		value(const std::vector<std::set<int>> &adj, bool is_directed = false)5179			: n_(static_cast<int>(adj.size())), adj_(adj),5180			  is_directed_(is_directed), add_1_(false), print_nm_(false),5181			  adj_built_(true) {5182			for (int u = 0; u < n_; ++u)5183				for (auto v : adj[u]) {5184					tgen_ensure(5185						0 <= v and v < n_,5186						"wgraph: value: vertices must be indexed in [0, n)");5187					// Undirected adjacency is symmetric: count each edge once5188					// (canonical u <= v). Directed: every out-edge appears5189					// once.5190					if (is_directed_ or u <= v)5191						edges_.emplace_back(u, v);5192				}5193		}51945195		// Creates value from `n`, `m`, and edge list. The edges are5196		// considered to be directed.5197		// Optimized for performance (lazy adjacency list; unordered_set dedup).5198		// O(m log m).5199		value(int n, const std::vector<std::pair<int, int>> &edges = {},5200			  bool is_directed = false)5201			: n_(n), edges_(), is_directed_(is_directed), add_1_(false),5202			  print_nm_(false), adj_built_(false) {5203			edges_.reserve(edges.size());5204			std::unordered_set<uint64_t> seen;5205			seen.reserve(edges.size() * 2 + 1);5206			for (auto [u, v] : edges) {5207				tgen_ensure(5208					0 <= std::min(u, v) and std::max(u, v) < n,5209					"wgraph: value: vertices must be indexed in [0, n)");5210				if (!is_directed_ and u > v)5211					std::swap(u, v);5212				uint64_t key = is_directed_ ? detail::directed_edge_key(u, v)5213											: detail::undirected_edge_key(u, v);5214				if (seen.insert(key).second)5215					edges_.emplace_back(u, v);5216			}5217		}5218		value(int n, const std::set<std::pair<int, int>> &edges,5219			  bool is_directed = false)5220			: value(5221				  n,5222				  std::vector<std::pair<int, int>>(edges.begin(), edges.end()),5223				  is_directed) {}5224		value(int n, const std::initializer_list<std::pair<int, int>> &edges,5225			  bool is_directed = false)5226			: value(n, std::vector<std::pair<int, int>>(edges), is_directed) {}52275228		// Creates graph from tree (undirected, same edges).5229		// O(n).5230		value(const typename wtree<VWeight, EWeight>::value &t)5231			: value(t.n(), t.edges(), false) {5232			if (t.vertex_weights().has_value()) {5233				vertex_weights_ = *t.vertex_weights();5234			}5235			if (t.edge_weights().has_value()) {5236				edge_weights_ = *t.edge_weights();5237			}5238		}52395240		// Weight type conversion.5241		// O(n + m).5242		template <typename NewVWeight, typename NewEWeight>5243		typename wgraph<NewVWeight, NewEWeight>::value5244		convert_weight_types() const {5245			tgen_ensure(!vertex_weights_.has_value() and5246							!edge_weights_.has_value(),5247						"wgraph: value: cannot convert weight type after "5248						"assigning weights");52495250			ensure_adj_built();5251			typename wgraph<NewVWeight, NewEWeight>::value new_graph(5252				adj_, is_directed_);5253			new_graph.is_directed_ = is_directed_;5254			new_graph.add_1_ = add_1_;5255			new_graph.print_nm_ = print_nm_;5256			return new_graph;5257		}52585259		// Fetches number of vertices.5260		int n() const { return n_; }52615262		// Fetches number of edges.5263		int m() const { return edges_.size(); }52645265		// Fetches if graph is directed;5266		bool is_directed() const { return is_directed_; }52675268		// Fetches a const ref. to adjacency list.5269		const std::vector<std::set<int>> &adj() const {5270			ensure_adj_built();5271			return adj_;5272		}52735274		// Fetches a const ref. to edge set.5275		const std::vector<std::pair<int, int>> &edges() const { return edges_; }52765277		// Fetches vertex weights.5278		const std::optional<std::vector<VWeight>> &vertex_weights() const {5279			return vertex_weights_;5280		}52815282		// Fetches edge weights.5283		const std::optional<std::vector<EWeight>> &edge_weights() const {5284			return edge_weights_;5285		}52865287		// Sets vertex weights.5288		// O(n + m).5289		template <typename NewVWeight = VWeight>5290		typename wgraph<NewVWeight, EWeight>::value set_vertex_weights(5291			const std::vector<NewVWeight> &vertex_weights) const {5292			tgen_ensure(static_cast<int>(vertex_weights.size()) == n(),5293						"wgraph: value: must give `n` vertex weights");52945295			auto new_graph = convert_weight_types<NewVWeight, EWeight>();5296			new_graph.vertex_weights_ = vertex_weights;5297			return new_graph;5298		}52995300		// Sets edge weights.5301		// O(n + m).5302		template <typename NewEWeight = EWeight>5303		typename wgraph<VWeight, NewEWeight>::value5304		set_edge_weights(const std::vector<NewEWeight> &edge_weights) const {5305			tgen_ensure(static_cast<int>(edge_weights.size()) == m(),5306						"wgraph: value: must give `m` edge weights");53075308			auto new_graph = convert_weight_types<VWeight, NewEWeight>();5309			new_graph.edge_weights_ = edge_weights;5310			return new_graph;5311		}53125313		// Enables edge-weighted mode before adding weighted edges5314		// incrementally. The graph must have no edges yet. O(1).5315		value &edge_weighted() {5316			tgen_ensure(m() == 0,5317						"wgraph: value: edge_weighted requires a graph with no "5318						"edges");5319			tgen_ensure(!edge_weights_.has_value(),5320						"wgraph: value: graph is already edge-weighted");53215322			edge_weights_ = std::vector<EWeight>();5323			return *this;5324		}53255326		// Adds 1 to vertex ids, for printing.5327		// O(1).5328		value &add_1() {5329			add_1_ = true;5330			return *this;5331		}53325333		// Prints `n m` on a new line before printing the edges.5334		// O(1).5335		value &print_nm() {5336			print_nm_ = true;5337			return *this;5338		}53395340		// Shuffles the graph's vertex labels (except those in `indices`,5341		// which keep their current label) and edge order. The change is5342		// applied eagerly to the underlying adjacency list, edge list,5343		// vertex weights and edge weights.5344		// O(n + m).5345		value &shuffle_except(std::set<int> indices) {5346			ensure_adj_built();5347			// Builds the relabeling: for each vertex `i`, `new_label[i]` is5348			// its new id. Vertices in `indices` keep their label; the others5349			// are permuted among themselves.5350			std::vector<int> new_label(n());5351			std::vector<int> shuffled;5352			for (int i = 0; i < n(); ++i) {5353				if (indices.count(i))5354					new_label[i] = i;5355				else5356					shuffled.push_back(i);5357			}5358			std::vector<int> targets = shuffled;5359			tgen::shuffle(targets.begin(), targets.end());5360			for (size_t k = 0; k < shuffled.size(); ++k)5361				new_label[shuffled[k]] = targets[k];53625363			// Rewrites adjacency list with new labels.5364			std::vector<std::set<int>> new_adj(n());5365			for (int u = 0; u < n(); ++u)5366				for (int v : adj_[u])5367					new_adj[new_label[u]].insert(new_label[v]);5368			adj_ = new_adj;53695370			// Rewrites edges with new labels (canonical undirected order).5371			for (auto &[u, v] : edges_) {5372				u = new_label[u];5373				v = new_label[v];5374				if (!is_directed_ and u > v)5375					std::swap(u, v);5376			}53775378			// Permutes vertex weights to match the new labels.5379			if (vertex_weights_.has_value()) {5380				std::vector<VWeight> new_vw(n());5381				for (int i = 0; i < n(); ++i)5382					new_vw[new_label[i]] = (*vertex_weights_)[i];5383				vertex_weights_ = new_vw;5384			}53855386			// Shuffles edge order, keeping edge weights aligned.53875388			std::vector<int> perm(edges_.size());5389			std::iota(perm.begin(), perm.end(), 0);5390			tgen::shuffle(perm.begin(), perm.end());53915392			std::vector<std::pair<int, int>> new_edges;5393			std::optional<std::vector<EWeight>> new_ew;5394			if (edge_weights_.has_value())5395				new_ew = std::vector<EWeight>();5396			for (int i : perm) {5397				new_edges.push_back(edges_[i]);5398				if (new_ew.has_value())5399					new_ew->push_back((*edge_weights_)[i]);5400			}54015402			edges_ = new_edges;5403			if (new_ew.has_value())5404				edge_weights_ = new_ew;54055406			return *this;5407		}54085409		// Shuffles the graph's vertices and edge order.5410		// O(n + m).5411		value &shuffle() { return shuffle_except({}); }54125413		// Adds `k` vertices to the graph (labeled n, n+1, ...n+k-1). Updates5414		// `n` accordingly.5415		// O(k) amortized.5416		value &add_vertices(int k, std::optional<std::vector<VWeight>>5417									   new_vertex_weights = std::nullopt) {5418			ensure_adj_built();5419			n_ += k;5420			adj_.resize(n());5421			if (new_vertex_weights.has_value()) {5422				tgen_ensure(vertex_weights().has_value(),5423							"wgraph: value: cannot add weighted vertices to "5424							"vertex-unweighted graph");5425				tgen_ensure(5426					static_cast<int>(new_vertex_weights->size()) == k,5427					"wgraph: value: number of vertex weights must be equal "5428					"to number of added vertices");54295430				vertex_weights_->insert(vertex_weights_->end(),5431										new_vertex_weights->begin(),5432										new_vertex_weights->end());5433			} else5434				tgen_ensure(!vertex_weights().has_value(),5435							"wgraph: value: cannot add unweighted vertices to "5436							"vertex-weighted graph");54375438			return *this;5439		}54405441		// Adds edge (u, v).5442		// O(log n) amortized.5443		value &add_edge(int u, int v, std::optional<EWeight> w = std::nullopt) {5444			ensure_adj_built();5445			tgen_ensure(0 <= std::min(u, v) and std::max(u, v) < n(),5446						"wgraph: value: vertex ids must be valid");54475448			if (!is_directed() and u > v)5449				std::swap(u, v);54505451			if (adj_[u].count(v))5452				return *this;54535454			adj_[u].insert(v);5455			if (!is_directed())5456				adj_[v].insert(u);5457			edges_.emplace_back(u, v);54585459			if (w.has_value()) {5460				tgen_ensure(edge_weights().has_value(),5461							"wgraph: value: cannot add weighted edge to "5462							"edge-unweighted graph");54635464				edge_weights_->push_back(*w);5465			} else5466				tgen_ensure(!edge_weights().has_value(),5467							"wgraph: value: cannot add unweighted edge to "5468							"edge-weighted graph");54695470			return *this;5471		}54725473		// Links graph with another `rhs`, adding the edge between u (in left5474		// graph) and v (in right graph). Ids for added vertices are updated5475		// accordingly.5476		// O(rhs.n + rhs.m * log n) amortized.5477		value &link(const value &rhs, int new_u, int new_v,5478					std::optional<EWeight> new_w = std::nullopt) {5479			tgen_ensure(0 <= new_u and new_u < n() and 0 <= new_v and5480							new_v < rhs.n(),5481						"wgraph: value: vertex ids must be valid");54825483			// Edges from right-hand side.5484			int shift = n();5485			add_vertices(rhs.n(), rhs.vertex_weights());5486			for (int i = 0; i < rhs.m(); ++i) {5487				auto [u, v] = rhs.edges()[i];5488				add_edge(shift + u, shift + v,5489						 rhs.edge_weights().has_value()5490							 ? std::optional<EWeight>((*rhs.edge_weights())[i])5491							 : std::nullopt);5492			}54935494			// New edge.5495			add_edge(new_u, shift + new_v, new_w);54965497			return *this;5498		}54995500		// Glues the graph with another `rhs` such that index_pairs[i].first is5501		// considered to be the same as index_pairs[i].second. Ids for added5502		// vertices are updated accordingly.5503		// O(rhs.n + rhs.m * log n) amortized.5504		value &glue(const value &rhs,5505					std::set<std::pair<int, int>> index_pairs) {5506			tgen_ensure(5507				is_directed() == rhs.is_directed(),5508				"wgraph: value: graphs must have the same is_directed value");55095510			// Checks validity of indices.5511			std::set<int> idx_left, idx_right;5512			std::vector<int> right_id_to_left(rhs.n(), -1);5513			for (auto [l, r] : index_pairs) {5514				tgen_ensure(5515					0 <= l and l < n() and 0 <= r and r < rhs.n(),5516					"wgraph: value: vertex indices to glue must be valid");5517				tgen_ensure(idx_left.count(l) == 0 and idx_right.count(r) == 0,5518							"wgraph: value: must not have repeated indices "5519							"on the same side to glue");55205521				idx_left.insert(l);5522				idx_right.insert(r);5523				right_id_to_left[r] = l;5524			}55255526			// Computes new ids of right vertices.5527			std::vector<int> new_right_id(rhs.n(), -1);5528			int intersection_lt = 0;5529			std::optional<std::vector<VWeight>> rhs_vertex_weights;5530			for (int i = 0; i < rhs.n(); ++i) {5531				if (right_id_to_left[i] != -1) {5532					// Is in intersection.5533					++intersection_lt;5534					new_right_id[i] = right_id_to_left[i];5535				} else {5536					// New id.5537					new_right_id[i] = n() + i - intersection_lt;5538					if (rhs.vertex_weights().has_value()) {5539						if (!rhs_vertex_weights.has_value())5540							rhs_vertex_weights = std::vector<VWeight>();5541						rhs_vertex_weights->push_back(5542							(*rhs.vertex_weights())[i]);5543					}5544				}5545			}55465547			// Adds new vertices and edges.5548			add_vertices(rhs.n() - intersection_lt, rhs_vertex_weights);5549			for (int i = 0; i < rhs.m(); ++i) {5550				auto [u, v] = rhs.edges()[i];5551				add_edge(new_right_id[u], new_right_id[v],5552						 rhs.edge_weights().has_value()5553							 ? std::optional<EWeight>((*rhs.edge_weights())[i])5554							 : std::nullopt);5555			}55565557			return *this;5558		}5559		value &glue(const value &rhs,5560					std::initializer_list<std::pair<int, int>> il) {5561			return glue(rhs, std::set<std::pair<int, int>>(il));5562		}55635564		// Glues the graph with another `rhs` at `indices`. That is, idx in5565		// `indices` are considered to be the same vertex. Ids for added5566		// vertices are updated accordingly.5567		// O(rhs.n + rhs.m * log n) amortized.5568		value &glue(const value &rhs, std::set<int> indices) {5569			std::set<std::pair<int, int>> index_pairs;5570			for (auto i : indices)5571				index_pairs.emplace(i, i);5572			return glue(rhs, index_pairs);5573		}5574		value &glue(const value &rhs, const std::initializer_list<int> &il) {5575			return glue(rhs, std::set<int>(il));5576		}55775578		// Disjoint union.5579		// Shifts ids from `rhs` graph by n().5580		// O(rhs.n + rhs.m * log n) amortized.5581		value &disjoint_union(const value &rhs) {5582			return glue(rhs, std::set<int>());5583		}55845585		// Computes uniformly random subgraph of graph with num_edges edges.5586		// O(n + m).5587		value &random_subgraph(int num_edges) {5588			tgen_ensure(5589				num_edges <= m(),5590				"wgraph: value: can choose at most `m` edges from graph");55915592			std::vector<std::pair<int, int>> new_edges;5593			std::optional<std::vector<EWeight>> new_edge_weights;55945595			int left = m();5596			for (int i = 0; i < m(); ++i) {5597				if (next(1, left--) <= num_edges) {5598					new_edges.push_back(edges()[i]);5599					if (edge_weights_.has_value()) {5600						if (!new_edge_weights.has_value())5601							new_edge_weights = std::vector<EWeight>();5602						new_edge_weights->push_back((*edge_weights())[i]);5603					}5604					--num_edges;5605				}5606			}56075608			edges_ = new_edges;5609			edge_weights_ = new_edge_weights;5610			rebuild_adj_from_edge_list();5611			return *this;5612		}56135614		// Computes a random (not uniform) subgraph with `num_edges` edges that5615		// keeps every connected component connected (does not increase the5616		// number of connected components).5617		// 1. Picks a spanning forest via randomized Prim.5618		// 2. Adds additional edges uniformly at random.5619		// O(n + m).5620		value &random_connected_subgraph(int num_edges) {5621			tgen_ensure(!is_directed_,5622						"wgraph: value: random_connected_subgraph is only for "5623						"undirected graphs");5624			tgen_ensure(5625				num_edges <= m(),5626				"wgraph: value: can choose at most `m` edges from graph");56275628			// Builds an incidence list: for each vertex, the (neighbor, edge5629			// index) pairs.5630			std::vector<std::vector<std::pair<int, int>>> incident(n());5631			for (int i = 0; i < m(); ++i) {5632				auto [u, v] = edges_[i];5633				incident[u].emplace_back(v, i);5634				incident[v].emplace_back(u, i);5635			}56365637			// Randomized Prim.5638			std::vector<bool> vis(n(), false);5639			std::vector<int> queue;5640			std::vector<bool> in_tree(m(), false);5641			int forest_edges = 0;56425643			for (int start = 0; start < n(); ++start) {5644				if (vis[start])5645					continue;5646				vis[start] = true;5647				queue.push_back(start);56485649				while (!queue.empty()) {5650					int i = tgen::next<int>(0, queue.size() - 1);5651					int u = queue[i];5652					std::swap(queue[i], queue.back());5653					queue.pop_back();56545655					for (auto [v, edge_idx] : incident[u]) {5656						if (!vis[v]) {5657							vis[v] = true;5658							queue.push_back(v);5659							in_tree[edge_idx] = true;5660							++forest_edges;5661						}5662					}5663				}5664			}5665			tgen_ensure(5666				num_edges >= forest_edges,5667				"wgraph: value: random_connected_subgraph needs at least "5668				"`n - c` edges, where `c` is the number of connected "5669				"components");56705671			// Splits edge indices into forest edges and the rest.5672			std::vector<int> tree_idx, rest_idx;5673			for (int i = 0; i < m(); ++i) {5674				if (in_tree[i])5675					tree_idx.push_back(i);5676				else5677					rest_idx.push_back(i);5678			}56795680			tgen::shuffle(rest_idx.begin(), rest_idx.end());56815682			std::vector<int> chosen_idx;5683			chosen_idx.insert(chosen_idx.end(), tree_idx.begin(),5684							  tree_idx.end());5685			chosen_idx.insert(chosen_idx.end(), rest_idx.begin(),5686							  rest_idx.begin() + num_edges - forest_edges);56875688			detail::tgen_ensure_against_bug(5689				static_cast<int>(chosen_idx.size()) == num_edges,5690				"wgraph: value: chose a wrong number of edges");56915692			std::vector<std::pair<int, int>> new_edges;5693			std::optional<std::vector<EWeight>> new_edge_weights;5694			if (edge_weights_.has_value())5695				new_edge_weights = std::vector<EWeight>();5696			for (int i : chosen_idx) {5697				new_edges.push_back(edges_[i]);5698				if (new_edge_weights.has_value())5699					new_edge_weights->push_back((*edge_weights_)[i]);5700			}57015702			edges_ = new_edges;5703			edge_weights_ = new_edge_weights;5704			rebuild_adj_from_edge_list();5705			return *this;5706		}57075708		// Complement. Self loops are maintained.5709		// O(n^2).5710		value operator!() const {5711			tgen_ensure(!edge_weights_.has_value(),5712						"wgraph: value: cannot compute complement of "5713						"edge-weighted graph");57145715			value complement = *this;5716			complement.ensure_adj_built();5717			std::vector<std::pair<int, int>> compl_edges;5718			for (int i = 0; i < complement.n_; ++i) {5719				std::set<int> complement_adj;5720				for (int j = 0; j < complement.n_; ++j) {5721					bool add_j = false;5722					if (j == i and complement.adj_[i].count(j))5723						add_j = true;5724					if (j != i and !complement.adj_[i].count(j))5725						add_j = true;57265727					if (add_j) {5728						complement_adj.insert(j);5729						// If i > j and !is_directed(), we don't add the edge.5730						if (i <= j or complement.is_directed_) {5731							compl_edges.emplace_back(i, j);5732						}5733					}5734				}5735				std::swap(complement.adj_[i], complement_adj);5736			}5737			std::swap(complement.edges_, compl_edges);57385739			return complement;5740		}57415742		// Concatenates two values.5743		// O(N + M log N), N = n + rhs.n, M = m + rhs.m.5744		value operator+(const value &rhs) const {5745			tgen_ensure(is_directed() == rhs.is_directed(),5746						"wgraph: value: graphs must have the same "5747						"is_directed value");57485749			tgen_ensure(vertex_weights().has_value() ==5750							rhs.vertex_weights().has_value(),5751						"wgraph: value: cannot concatenate vertex-weighted "5752						"wgraph to unweighted");5753			tgen_ensure(edge_weights().has_value() ==5754							rhs.edge_weights().has_value(),5755						"wgraph: value: cannot concatenate edge-weighted "5756						"wgraph to unweighted");57575758			value concat = *this;5759			concat.glue(rhs, std::set<std::pair<int, int>>());5760			concat.add_1_ = add_1_ | rhs.add_1_;5761			concat.print_nm_ = print_nm_ | rhs.print_nm_;57625763			return concat;5764		}57655766		// Prints to std::ostream.5767		// O(n + m).5768		friend std::ostream &operator<<(std::ostream &out, const value &val) {5769			// Prints `n` and `m`.5770			if (val.print_nm_)5771				out << val.n() << " " << val.m() << '\n';57725773			// Prints vertex weights.5774			if (val.vertex_weights()) {5775				for (int i = 0; i < val.n(); ++i) {5776					if (i > 0)5777						out << " ";5778					out << (*val.vertex_weights())[i];5779				}5780				out << '\n';5781			}57825783			// Prints edges.5784			for (int i = 0; i < val.m(); ++i) {5785				auto [u, v] = val.edges()[i];5786				out << (u + val.add_1_) << " " << (v + val.add_1_);57875788				// Edge weight.5789				if (val.edge_weights().has_value())5790					out << " " << (*val.edge_weights())[i];57915792				out << '\n';5793			}57945795			return out;5796		}57975798		// Gets a std::tuple<n, m, adj> representing the value.5799		std::tuple<int, int, std::vector<std::set<int>>> to_std() const {5800			ensure_adj_built();5801			return std_type(n_, m(), adj_);5802		}58035804	  private:5805		// Rebuilds adjacency from edges_ after replacing the edge list (e.g.5806		// subgraph operations).5807		// O(m log n).5808		void rebuild_adj_from_edge_list() {5809			adj_.assign(n_, {});5810			for (auto [u, v] : edges_) {5811				adj_[u].insert(v);5812				if (!is_directed_)5813					adj_[v].insert(u);5814			}5815			adj_built_ = true;5816		}58175818		// Builds adj_ from edges_ on first use.5819		// O(1) if already built; O(m log n) otherwise.5820		void ensure_adj_built() const {5821			if (adj_built_)5822				return;5823			const_cast<value *>(this)->rebuild_adj_from_edge_list();5824		}5825	};58265827	// Adds all edges from `rhs` as preset edges.5828	// O(rhs.m * log m).5829	wgraph &add_edges_from(const value &rhs) {5830		tgen_ensure(is_directed_ == rhs.is_directed(),5831					"wgraph: graphs must have the same is_directed value");58325833		for (auto [u, v] : rhs.edges())5834			add_edge(u, v);5835		return *this;5836	}58375838	// Generates graph value.5839	// Optimized for performance: dense no-preset graphs use index sampling;5840	// otherwise gen_remaining_edges.5841	// O(n + m log^2 n) expected.5842	value gen() const {5843		detail::tgen_ensure_against_bug(static_cast<int>(edges_.size()) <= m_,5844										"wgraph: too many edges were added");58455846		// All edges already added.5847		if (static_cast<int>(edges_.size()) == m_)5848			return value(n_, edges_, is_directed_);58495850		// Splits into two cases to optimize performance.58515852		// No presets and m > max_edges / 2: sample m distinct edge indices.5853		if (auto indexed = try_gen_by_edge_index())5854			return *indexed;58555856		// Otherwise: fill preset edges up to m_ with uniform random edges.5857		return gen_remaining_edges(5858			std::vector<std::pair<int, int>>(edges_.begin(), edges_.end()));5859	}58605861	// Gets a (not uniformly) random connected undirected graph.5862	// 1. Preset edges induce a spanning forest on their components.5863	// 2. Then, uniformly random edges between components are added.5864	// 3. Remaining edges are added uniformly at random.5865	// O(n + m log^2 n) expected.5866	value get_connected() const {5867		tgen_ensure(!is_directed_,5868					"wgraph: get_connected is only for undirected graphs");5869		tgen_ensure(m_ >= n_ - 1,5870					"wgraph: connected graph needs at least n - 1 edges");58715872		std::vector<std::pair<int, int>> edges;5873		edges.reserve(m_);58745875		if (edges_.empty()) {5876			if (n_ > 1) {5877				std::vector<int> prufer(n_ - 2);5878				for (int i = 0; i < n_ - 2; ++i)5879					prufer[i] = next<int>(0, n_ - 1);5880				for (auto [u, v] : detail::edges_from_prufer(std::move(prufer)))5881					edges.emplace_back(u, v);5882			}5883		} else {5884			edges.assign(edges_.begin(), edges_.end());58855886			std::vector<std::vector<int>> adj(n_);5887			for (auto [u, v] : edges_) {5888				adj[u].push_back(v);5889				adj[v].push_back(u);5890			}58915892			std::vector<int> comp_size;5893			std::vector<std::vector<int>> component_ids;5894			std::vector<bool> vis(n_, false);5895			std::queue<int> q;58965897			for (int i = 0; i < n_; ++i) {5898				if (vis[i])5899					continue;59005901				vis[i] = true;5902				q.push(i);5903				comp_size.push_back(0);5904				component_ids.emplace_back();5905				while (q.size()) {5906					int u = q.front();5907					q.pop();5908					++comp_size.back();5909					component_ids.back().push_back(u);5910					for (int v : adj[u]) {5911						if (!vis[v]) {5912							vis[v] = true;5913							q.push(v);5914						}5915					}5916				}5917			}59185919			if (component_ids.size() > 1) {5920				std::vector<int> prufer_values =5921					many_by_distribution(component_ids.size() - 2, comp_size);5922				for (auto [u, v] :5923					 detail::edges_from_prufer(std::move(prufer_values)))5924					edges.emplace_back(pick(component_ids[u]),5925									   pick(component_ids[v]));5926			}5927		}59285929		return gen_remaining_edges(std::move(edges));5930	}59315932	// Gets a (not uniformly) random directed acyclic graph.5933	// 1. Randomized Kahn (uniform choice among indegree-0 vertices) yields a5934	//    random topological order of the preset edges (which must be acyclic).5935	// 2. Extra edges are sampled randomly using the order.5936	// With no preset edges: sample a random graph then orient acyclically.5937	// Optimized for performance (distinct upper-triangle edge-index sampling;5938	// rejection instead of pair::distinct for preset edges).5939	// O(n + m log^2 n) expected.5940	value get_acyclic() const {5941		tgen_ensure(is_directed_,5942					"wgraph: get_acyclic is only for directed graphs");59435944		if (edges_.empty()) {5945			std::vector<int> order(n_);5946			std::iota(order.begin(), order.end(), 0);5947			for (int i = n_ - 1; i > 0; --i)5948				std::swap(order[i], order[next(0, i)]);59495950			const long long max_pairs =5951				static_cast<long long>(n_) * (n_ - 1) / 2;5952			tgen_ensure(m_ <= max_pairs,5953						"wgraph: not enough edges to generate");59545955			std::vector<std::pair<int, int>> edges;5956			edges.reserve(m_);5957			for (long long idx : distinct_range<long long>(0, max_pairs - 1)5958									 .gen_list(m_)5959									 .to_std()) {5960				auto [i, j] = detail::decode_undirected_simple_edge(n_, idx);5961				edges.emplace_back(order[i], order[j]);5962			}5963			return value(n_, edges, true);5964		}59655966		std::vector<std::vector<int>> adj(n_);5967		std::vector<int> indeg(n_, 0);5968		for (auto [u, v] : edges_) {5969			adj[u].push_back(v);5970			++indeg[v];5971		}59725973		std::vector<int> available;5974		for (int i = 0; i < n_; ++i)5975			if (indeg[i] == 0)5976				available.push_back(i);59775978		// Random topological order using randomized Kahn's algorithm.5979		std::vector<int> order;5980		while (!available.empty()) {5981			int idx = next(0, static_cast<int>(available.size()) - 1);5982			int u = available[idx];5983			std::swap(available[idx], available.back());5984			available.pop_back();59855986			order.push_back(u);5987			for (int v : adj[u])5988				if (--indeg[v] == 0)5989					available.push_back(v);5990		}59915992		tgen_ensure(static_cast<int>(order.size()) == n_,5993					"wgraph: preset edges contain a directed cycle");59945995		value acyclic(n_, edges_, true);59965997		// Generates final edges.59985999		detail::tgen_ensure_against_bug(acyclic.m() <= m_,6000										"wgraph: too many edges were added");60016002		if (acyclic.m() < m_) {6003			std::vector<int> order_pos(n_);6004			for (int i = 0; i < n_; ++i)6005				order_pos[order[i]] = i;60066007			std::unordered_set<uint64_t> seen;6008			seen.reserve(m_ * 2);6009			for (auto [u, v] : acyclic.edges())6010				seen.insert(6011					detail::undirected_edge_key(order_pos[u], order_pos[v]));60126013			const long long max_pairs =6014				static_cast<long long>(n_) * (n_ - 1) / 2;6015			while (acyclic.m() < m_) {6016				std::pair<int, int> edge;6017				if (!detail::try_generate_distinct(seen, [&] {6018						long long idx = next<long long>(0, max_pairs - 1);6019						edge = detail::decode_undirected_simple_edge(n_, idx);6020						return detail::undirected_edge_key(edge.first,6021														   edge.second);6022					}))6023					throw detail::error("wgraph: not enough edges to generate");6024				acyclic.add_edge(order[edge.first], order[edge.second]);6025			}6026		}60276028		return acyclic;6029	}60306031	// Generates a (not uniformly) random skewed connected graph.6032	// 1. Builds the same skewed labeled tree as wtree::gen_skewed(n,6033	//    elongation)(root 0, parent(i) = wnext(i, elongation) for i >= 1).6034	//    If is_directed, tree edges are oriented down the tree.6035	// 2. Adds the remaining edges: pick an endpoint u uniformly;6036	//    pick k uniformly in [1, spread]; walk from u toward the root k6037	//    times along tree parents to get v; add edge (v, u).6038	// If elongation is small, generates a graph with small diameter.6039	// If elongation is large, generates a graph with large diameter, with6040	// vertices 0 and n-1 being far apart.6041	// O(n + m log n) if spread is O(1);6042	// O(n log n + m log^2 n) expected otherwise.6043	static value gen_skewed(int n, int m, int elongation, int spread,6044							bool is_directed = false) {6045		tgen_ensure(6046			m >= n - 1,6047			"wgraph: skewed graph needs at least n - 1 edges to be connected");6048		tgen_ensure(spread >= 2,6049					"wgraph: gen_skewed spread must be at least 2");60506051		value skewed(n, {}, is_directed);60526053		std::vector<int> parent(n), depth(n, 0);6054		parent[0] = 0;6055		for (int i = 1; i < n; ++i) {6056			int p = wnext<int>(i, elongation);6057			parent[i] = p;6058			depth[i] = depth[p] + 1;6059			skewed.add_edge(p, i);6060		}60616062		const int extra = m - (n - 1);6063		if (extra == 0)6064			return skewed;60656066		// If spread is large, use binary lifting to find the ancestor.6067		// Otherwise, enumerate O(n * spread) ancestor edges and sample6068		// directly.6069		constexpr int naive_ancestor_spread = 20;60706071		if (spread <= naive_ancestor_spread) {6072			std::vector<std::pair<int, int>> candidates;6073			candidates.reserve(n * spread);6074			for (int u = 0; u < n; ++u) {6075				int max_k = std::min(spread, depth[u]);6076				if (max_k < 2)6077					continue;6078				int v = parent[u];6079				for (int k = 2; k <= max_k; ++k) {6080					v = parent[v];6081					candidates.emplace_back(v, u);6082				}6083			}60846085			tgen_ensure(extra <= static_cast<int>(candidates.size()),6086						"wgraph: not enough edges to generate");60876088			for (auto [v, u] : choose(candidates, extra))6089				skewed.add_edge(v, u);6090		} else {6091			// Binary lifting.6092			int lg = 1;6093			while ((1 << lg) <= n)6094				++lg;60956096			std::vector<std::vector<int>> up(lg, std::vector<int>(n));6097			for (int v = 0; v < n; ++v)6098				up[0][v] = parent[v];6099			for (int j = 1; j < lg; ++j)6100				for (int v = 0; v < n; ++v)6101					up[j][v] = up[j - 1][up[j - 1][v]];61026103			// Creates uniform generator of edges (u, v) such that v is ancestor6104			// of u. For that, every u has depth[u]-1 choices for v, so we6105			// weight u by min(spread - 1, depth[u] - 1). After that we can6106			// just pick the ancestor uniformly.6107			std::vector<int> distribution = depth;6108			for (int &d : distribution)6109				d = std::max(0, std::min(spread - 1, d - 1));6110			weighted_sampler vertex_choice(distribution);6111			distinct extra_edges([&]() -> std::pair<int, int> {6112				int u = vertex_choice.next();6113				int k = next(2, spread);6114				int v = u;6115				for (int j = 0; j < lg; ++j)6116					if (k >> j & 1)6117						v = up[j][v];6118				return {v, u};6119			});61206121			while (skewed.m() < m) {6122				std::pair<int, int> edge;6123				try {6124					edge = extra_edges.gen();6125				} catch (const std::runtime_error &e) {6126					if (std::string(e.what()) ==6127						"tgen: distinct: no more distinct values")6128						throw detail::error(6129							"wgraph: not enough edges to generate");6130					throw e;6131				}61326133				skewed.add_edge(edge.first, edge.second);6134			}6135		}61366137		return skewed;6138	}61396140	// Generates a random bipartite graph. The first side has vertices6141	// 0 .. n1-1, the second n1 .. n1+n2-1.6142	// Uniform when connected is false (distinct cross-edge indices).6143	// When connected, bipartite Prüfer + rejection fill; not uniform over6144	// connected bipartite graphs.6145	// O(n1 + n2 + m log(n1 * n2)) expected.6146	static value gen_bipartite(int n1, int n2, int m, bool connected = false) {6147		tgen_ensure(m >= 0, "wgraph: number of edges must be nonnegative");6148		long long num_edges = 1LL * n1 * n2;6149		tgen_ensure(m <= num_edges,6150					"wgraph: bipartite graph has at most n1 * n2 edges");6151		if (connected)6152			tgen_ensure(6153				m >= n1 + n2 - 1,6154				"wgraph: connected bipartite graph needs at least n1 + n2 - 1 "6155				"edges");61566157		if (!connected) {6158			std::vector<std::pair<int, int>> edges;6159			edges.reserve(m);6160			for (long long idx : distinct_range<long long>(0, num_edges - 1)6161									 .gen_list(m)6162									 .to_std())6163				edges.emplace_back(static_cast<int>(idx / n2),6164								   n1 + static_cast<int>(idx % n2));6165			return value(n1 + n2, std::move(edges), false);6166		}61676168		std::unordered_set<uint64_t> used_edges;6169		used_edges.reserve(m * 2);6170		std::vector<std::pair<int, int>> edges;6171		edges.reserve(m);61726173		auto pack_edge = [](int u, int v) -> uint64_t {6174			if (u > v)6175				std::swap(u, v);6176			return (static_cast<uint64_t>(u) << 32) | static_cast<uint32_t>(v);6177		};61786179		if (n1 > 0 and n2 > 0) {6180			std::vector<int> prufer(n1 + n2 - 2);6181			for (int i = 0; i < n2 - 1; ++i)6182				prufer[i] = next(0, n1 - 1);6183			for (int i = 0; i < n1 - 1; ++i)6184				prufer[n2 - 1 + i] = next(n1, n1 + n2 - 1);6185			shuffle(prufer.begin(), prufer.end());6186			for (auto [u, v] : detail::edges_from_prufer(std::move(prufer))) {6187				if (u > v)6188					std::swap(u, v);6189				if (used_edges.insert(pack_edge(u, v)).second)6190					edges.emplace_back(u, v);6191			}6192			detail::tgen_ensure_against_bug(6193				used_edges.size() == size_t(n1 + n2 - 1),6194				"wgraph: invalid bipartite spanning tree size");6195		}61966197		while (edges.size() < size_t(m)) {6198			int u = next(0, n1 - 1);6199			int v = next(n1, n1 + n2 - 1);6200			if (used_edges.insert(pack_edge(u, v)).second)6201				edges.emplace_back(u, v);6202		}62036204		return value(n1 + n2, std::move(edges), false);6205	}62066207  private:6208	// If this generator has no preset edges and m is large relative to the6209	// maximum edge count, sample by distinct edge index. Otherwise6210	// std::nullopt.6211	// Optimized for performance (index sampling instead of rejection).6212	// O(m log n).6213	std::optional<value> try_gen_by_edge_index() const {6214		if (!edges_.empty())6215			return std::nullopt;62166217		long long max_edges =6218			detail::max_graph_edges(n_, is_directed_, has_self_loops_);6219		if (m_ > max_edges)6220			throw detail::error("wgraph: not enough edges to generate");6221		if (max_edges <= 0 or 2LL * m_ <= max_edges)6222			return std::nullopt;62236224		std::vector<std::pair<int, int>> edges;6225		edges.reserve(m_);6226		for (long long idx :6227			 distinct_range<long long>(0, max_edges - 1).gen_list(m_).to_std())6228			edges.push_back(detail::decode_graph_edge_index(6229				n_, idx, is_directed_, has_self_loops_));62306231		return value(n_, edges, is_directed_);6232	}62336234	// Fills `edges` up to m_ with uniform random edges not already present.6235	// Optimized for performance (uint64 edge keys + try_generate_distinct).6236	// O(m log^2 n) expected.6237	value gen_remaining_edges(std::vector<std::pair<int, int>> edges) const {6238		detail::tgen_ensure_against_bug(static_cast<int>(edges.size()) <= m_,6239										"wgraph: too many edges were added");62406241		if (static_cast<int>(edges.size()) == m_)6242			return value(n_, edges, is_directed_);62436244		edges.reserve(m_);62456246		std::unordered_set<uint64_t> seen;6247		seen.reserve(m_ * 2);6248		for (auto [u, v] : edges) {6249			if (!is_directed_ and u > v)6250				std::swap(u, v);6251			seen.insert(is_directed_ ? detail::directed_edge_key(u, v)6252									 : detail::undirected_edge_key(u, v));6253		}62546255		while (static_cast<int>(edges.size()) < m_) {6256			std::pair<int, int> edge;6257			if (!detail::try_generate_distinct(seen, [&] {6258					edge = detail::get_random_graph_edge(n_, is_directed_,6259														 has_self_loops_);6260					if (!is_directed_ and edge.first > edge.second)6261						std::swap(edge.first, edge.second);6262					return is_directed_ ? detail::directed_edge_key(edge.first,6263																	edge.second)6264										: detail::undirected_edge_key(6265											  edge.first, edge.second);6266				}))6267				throw detail::error("wgraph: not enough edges to generate");6268			edges.emplace_back(edge);6269		}62706271		return value(n_, edges, is_directed_);6272	}6273};62746275// Implementation of wtree::value constructor from wgraph.6276// O(n + m alpha(n)).6277template <typename VWeight, typename EWeight>6278wtree<VWeight, EWeight>::value::value(6279	const typename wgraph<VWeight, EWeight>::value &g)6280	: n_(g.n()), adj_(g.n()), add_1_(false), print_n_(false), dsu_(g.n()) {6281	tgen_ensure(g.n() > 0, "wtree: value: graph must have at least one vertex");6282	tgen_ensure(!g.is_directed(),6283				"wtree: value: graph must be undirected to form a tree");62846285	if (g.vertex_weights().has_value())6286		vertex_weights_ = *g.vertex_weights();6287	if (g.edge_weights().has_value())6288		edge_weights_ = std::vector<EWeight>();62896290	if (n_ == 1)6291		return;62926293	std::vector<int> order(g.m());6294	std::iota(order.begin(), order.end(), 0);6295	tgen::shuffle(order.begin(), order.end());62966297	std::vector<std::pair<int, int>> tree_edges;6298	tree_edges.reserve(n_ - 1);62996300	for (int i : order) {6301		auto [u, v] = g.edges()[i];6302		if (!dsu_.unite(u, v))6303			continue;6304		if (u > v)6305			std::swap(u, v);63066307		tree_edges.emplace_back(u, v);6308		adj_[u].insert(v);6309		adj_[v].insert(u);6310		if (edge_weights_.has_value())6311			edge_weights_->push_back((*g.edge_weights())[i]);6312		if (static_cast<int>(tree_edges.size()) == n_ - 1)6313			break;6314	}63156316	tgen_ensure(static_cast<int>(tree_edges.size()) == n_ - 1,6317				"wtree: value: graph must be connected to form a tree");63186319	edges_ = std::move(tree_edges);6320}63216322/*6323 * Other types of weighted-ness.6324 */63256326// Vertex weighted graph.6327template <typename VWeight> using vgraph = wgraph<VWeight, int>;63286329// Edge weighted graph.6330template <typename EWeight> using egraph = wgraph<int, EWeight>;63316332// Unweighted graph.6333using graph = wgraph<int, int>;63346335/*6336 * Standard graphs.6337 */63386339// Complete.6340// O(n^2).6341inline graph::value K(int n) { return graph(n, n * (n - 1) / 2).gen(); }63426343// Path.6344// Path with `n` vertices. The edges of the path are 0 and n-1.6345// If directed, edges are i -> i+1 for i in [0, n-2).6346// O(n).6347inline graph::value P(int n, bool is_directed = false) {6348	graph g(n, n - 1, is_directed);6349	for (int i = 0; i + 1 < n; ++i)6350		g.add_edge(i, i + 1);6351	return g.gen();6352}63536354// Cycle.6355// n >= 3.6356// If directed, edges are i -> (i+1) % n.6357// O(n).6358inline graph::value C(int n, bool is_directed = false) {6359	tgen_ensure(n >= 3, "graph: cycle size must be at least 3");63606361	graph g(n, n, is_directed);6362	for (int i = 0; i < n; ++i)6363		g.add_edge(i, (i + 1) % n);6364	return g.gen();6365}63666367// Complete bipartite.6368// The first side has vertices `0` to `n1-1`, the second side has vertices `n1`6369// to `n1+n2-1`.6370// O(n1 * n2).6371inline graph::value K(int n1, int n2) {6372	graph g(n1 + n2, static_cast<long long>(n1) * n2);6373	for (int i = 0; i < n1; ++i)6374		for (int j = 0; j < n2; ++j)6375			g.add_edge(i, n1 + j);6376	return g.gen();6377}63786379// Star.6380// The center is vertex 0.6381// O(n).6382inline graph::value S(int n) { return K(1, n - 1); }63836384/****************6385 *              *6386 *   GEOMETRY   *6387 *              *6388 ****************/63896390namespace geometry {63916392// Point on the plane with coordinates of type T.6393template <typename T> struct point {6394	static_assert(std::is_arithmetic_v<T>,6395				  "point requires an arithmetic coordinate type");63966397	// Dot/cross product type: __int128 for T = long long, long long for other6398	// integral T, T for floating-point.6399	using product_t = std::conditional_t<6400		std::is_same_v<T, long long>, detail::i128,6401		std::conditional_t<std::is_integral_v<T>, long long, T>>;64026403	// x and y coordinates.6404	T x_, y_;64056406	// Constructs a point with coordinates x and y.6407	point(T x = 0, T y = 0) : x_(x), y_(y) {}64086409	// Returns the x coordinate.6410	T x() const { return x_; }64116412	// Returns the y coordinate.6413	T y() const { return y_; }64146415	// Equality of coordinates, with epsilon-based equality for floating-point6416	// coordinates (tolerance 1e-9).6417	static bool coord_eq(T a, T b) {6418		if constexpr (std::is_integral_v<T>)6419			return a == b;6420		constexpr T eps = T(1e-9);6421		T d = a - b;6422		return d >= -eps and d <= eps;6423	}64246425	// Lexicographic order (by x, then y).6426	bool operator<(const point &p) const {6427		if (!coord_eq(x_, p.x()))6428			return x_ < p.x();6429		return y_ < p.y();6430	}64316432	// Equality of coordinates.6433	bool operator==(const point &p) const {6434		return coord_eq(x_, p.x()) and coord_eq(y_, p.y());6435	}64366437	// Vector addition.6438	point operator+(const point &p) const {6439		return point(x_ + p.x(), y_ + p.y());6440	}64416442	// Vector subtraction.6443	point operator-(const point &p) const {6444		return point(x_ - p.x(), y_ - p.y());6445	}64466447	// Scalar multiplication.6448	point operator*(T c) const { return point(x_ * c, y_ * c); }64496450	// Dot product.6451	product_t operator*(const point &p) const {6452		if constexpr (std::is_floating_point_v<T>)6453			return x_ * p.x() + y_ * p.y();6454		return product_t(x_) * p.x() + product_t(y_) * p.y();6455	}64566457	// Cross product (signed area of the parallelogram).6458	product_t operator^(const point &p) const {6459		if constexpr (std::is_floating_point_v<T>)6460			return x_ * p.y() - y_ * p.x();6461		return product_t(x_) * p.y() - product_t(y_) * p.x();6462	}64636464	// Prints the point as "x y".6465	friend std::ostream &operator<<(std::ostream &out, const point &p) {6466		return out << p.x() << ' ' << p.y();6467	}6468};64696470// Generates n distinct integer points in [min_coord, max_coord]^2 with no three6471// collinear.6472// O(n).6473inline std::vector<point<long long>>6474random_points_general_position(int n, long long min_coord,6475							   long long max_coord) {6476	tgen_ensure(n > 0,6477				"geometry: random_points_general_position: n must be positive");6478	tgen_ensure(max_coord >= min_coord,6479				"geometry: random_points_general_position: min_coord must be "6480				"at most max_coord");6481	tgen_ensure(6482		static_cast<detail::i128>(max_coord) - min_coord <=6483			std::numeric_limits<long long>::max(),6484		"geometry: random_points_general_position: coordinate range too large");6485	uint64_t width = max_coord - min_coord;6486	uint64_t p = math::prime_from(2 * n);64876488	// Requires width >= p - 1 because sheared coordinates lie in [0, p - 1].6489	tgen_ensure(width >= p - 1,6490				"geometry: random_points_general_position: coordinate range "6491				"too small for n");64926493	// Base set {(x, x^-1 mod p) : x \in {1, ..., p - 1}}.6494	//6495	// Over F_p, y = x^-1 is a rational map on nonzero x, so any line meets the6496	// graph in at most two points. Distinct x therefore give distinct points;6497	// no three of them can be collinear in the integer plane.6498	std::vector<uint64_t> x_range(p - 1);6499	std::iota(x_range.begin(), x_range.end(), 1);6500	shuffle(x_range.begin(), x_range.end());6501	std::vector<detail::i128> bx(n), by(n);6502	for (int i = 0; i < n; ++i) {6503		uint64_t x = x_range[i];6504		bx[i] = x;6505		by[i] = math::modular_inverse(x, p);6506	}65076508	// Applies a random element of SL(2, p) by composing several elementary6509	// shear matrices over F_p.6510	//6511	// Each step applies either6512	//     [1 r]       or       [1 0]6513	//     [0 1]                [r 1]6514	//6515	// with r \in {-2, -1, 1, 2}, and all arithmetic performed modulo p.6516	// After all iterations, the resulting transformation is6517	// A = S_k S_{k-1} ... S_16518	//6519	// where every S_i has determinant 1. Therefore det(A) = 1 (mod p),6520	// so A is invertible.6521	//6522	// Invertible affine transformations of F_p^2 preserve collinearity and6523	// non-collinearity. Since the base set {(x, x^-1 mod p)} contains no6524	// three collinear points, the transformed set also contains no three6525	// collinear points.6526	const int num_shears = 8;6527	std::vector<detail::i128> lin_x = bx, lin_y = by;65286529	for (int it = 0; it < num_shears; ++it) {6530		bool vertical_shear = next(2) == 0;6531		int shear_r = pick({-2, -1, 1, 2});65326533		for (int i = 0; i < n; ++i) {6534			if (vertical_shear)6535				lin_x[i] = (lin_x[i] + shear_r * lin_y[i]) % p;6536			else6537				lin_y[i] = (lin_y[i] + shear_r * lin_x[i]) % p;65386539			if (lin_x[i] < 0)6540				lin_x[i] += p;6541			if (lin_y[i] < 0)6542				lin_y[i] += p;6543		}6544	}65456546	detail::i128 min_x = lin_x[0], max_x = lin_x[0], min_y = lin_y[0],6547				 max_y = lin_y[0];6548	for (int i = 1; i < n; ++i) {6549		min_x = std::min(min_x, lin_x[i]);6550		max_x = std::max(max_x, lin_x[i]);6551		min_y = std::min(min_y, lin_y[i]);6552		max_y = std::max(max_y, lin_y[i]);6553	}65546555	long long x_shift =6556		min_coord - min_x + next<long long>(0, width - (max_x - min_x));6557	long long y_shift =6558		min_coord - min_y + next<long long>(0, width - (max_y - min_y));65596560	std::vector<point<long long>> pts;6561	for (int i = 0; i < n; ++i)6562		pts.emplace_back(lin_x[i] + x_shift, lin_y[i] + y_shift);6563	return pts;6564}65656566namespace detail {65676568using i128 = tgen::detail::i128;65696570// Signed area of triangle (a, b, p); positive iff (a, b, p) are in6571// counterclockwise order. 0 iff (a, b, p) are collinear. O(1).6572inline i128 ccw(const point<long long> &a, const point<long long> &b,6573				const point<long long> &p) {6574	return (static_cast<i128>(b.x()) - a.x()) *6575			   (static_cast<i128>(p.y()) - a.y()) -6576		   (static_cast<i128>(b.y()) - a.y()) *6577			   (static_cast<i128>(p.x()) - a.x());6578}65796580// Integer projection of P onto line AB (A and B need not be distinct).6581inline i128 proj_on_ab(const point<long long> &P, const point<long long> &A,6582					   const point<long long> &B) {6583	return (P - A) * (B - A);6584}65856586// In-place Hamiltonian path on points[left..right-1] with points[left]6587// start and points[right-1] end.6588// O(n log n) expected if points are "random", O(n^2) worst case.6589inline void conquer(std::vector<point<long long>> &points, int left,6590					int right) {6591	if (right - left <= 3)6592		return;65936594	point<long long> A = points[left], B = points[right - 1];65956596	// If all points are collinear, sort them properly and return.6597	bool all_collinear = true;6598	for (int k = left + 1; k < right - 1; ++k) {6599		if (ccw(A, B, points[k]) != 0) {6600			all_collinear = false;6601			break;6602		}6603	}6604	if (all_collinear) {6605		std::sort(points.begin() + left, points.begin() + right,6606				  [&](const point<long long> &P, const point<long long> &Q) {6607					  return proj_on_ab(P, A, B) < proj_on_ab(Q, A, B);6608				  });6609		return;6610	}66116612	// Choses a pivot that is not collinear with A and B.6613	std::vector<int> candidates;6614	for (int k = left + 1; k < right - 1; ++k) {6615		if (ccw(A, B, points[k]) != 0)6616			candidates.push_back(k);6617	}6618	int ci = candidates[next(0, static_cast<int>(candidates.size()) - 1)];6619	point<long long> C = points[ci];66206621	uint64_t wa = next<uint64_t>(1, std::numeric_limits<uint64_t>::max());6622	uint64_t wb = next<uint64_t>(1, std::numeric_limits<uint64_t>::max());6623	bool a_on_positive = ccw(C, A, B) < 0;66246625	// Classify interior points into two sides of the wedge A-C-B for partition.6626	// Collinear points on AB are tie-broken along the segment.6627	i128 proj_sum = proj_on_ab(A, A, B) + proj_on_ab(B, A, B);6628	auto is_positive = [&](const point<long long> &P) -> bool {6629		i128 s = wa * ccw(C, A, P) + wb * ccw(C, B, P);6630		// Weighted wedge side of P w.r.t. C, A, B.6631		if (s != 0)6632			return s > 0;6633		// P is on line AB: split by projection past the midpoint.6634		return 2 * proj_on_ab(P, A, B) > proj_sum;6635	};66366637	// Holds C at points[right-2] while classifying interior points in6638	// [left+1, right-3].6639	if (ci != right - 2)6640		std::swap(points[ci], points[right - 2]);66416642	int i = left + 1;6643	int j = right - 3;6644	while (i < j) {6645		if (is_positive(points[i]) == a_on_positive)6646			++i;6647		else if (is_positive(points[j]) != a_on_positive)6648			--j;6649		else {6650			std::swap(points[i], points[j]);6651			++i;6652			--j;6653		}6654	}66556656	// After partition:6657	// points[left]=A | (A,C)... | C | (C,B)... | points[right-1]=B.66586659	// After the swap, p is the index of C (pivot between the two subpaths).6660	int p = i;6661	if (i == j and is_positive(points[i]) == a_on_positive)6662		++p;6663	std::swap(points[p], points[right - 2]);66646665	// Path A -> C.6666	conquer(points, left, p + 1);6667	// Path C -> B.6668	conquer(points, p, right);6669}66706671// Samples k sorted distinct integers from [left, right] uniformly.6672// Optimized for performance (pool partial Fisher–Yates or complement path for6673// modest ranges; sparse-map fallback otherwise).6674// O(k log k); O(right - left) memory when the range is modest.6675inline std::vector<long long>6676sample_sorted_distinct_in_range(int k, long long left, long long right) {6677	long long universe = right - left + 1;6678	std::vector<long long> res;6679	res.reserve(k);6680	if (k == 0)6681		return res;66826683	constexpr long long pool_threshold = 8'000'000;6684	constexpr long long pool_always_below = 500'000;66856686	if (universe <= pool_threshold and6687		(universe <= pool_always_below or k >= universe / 4)) {6688		size_t u = universe;6689		size_t ks = k;6690		std::vector<long long> pool(u);6691		std::iota(pool.begin(), pool.end(), left);6692		size_t m = ks <= u / 2 ? ks : u - ks;6693		for (size_t i = 0; i < m; ++i) {6694			size_t j = next<size_t>(i, u - 1);6695			std::swap(pool[i], pool[j]);6696		}6697		if (ks <= u / 2) {6698			res.assign(pool.begin(), pool.begin() + ks);6699			std::sort(res.begin(), res.end());6700		} else {6701			std::vector<char> excluded(u, 0);6702			for (size_t i = 0; i < m; ++i)6703				excluded[pool[i] - left] = 1;6704			for (long long v = left; v <= right; ++v)6705				if (!excluded[v - left])6706					res.push_back(v);6707		}6708	} else {6709		std::unordered_map<long long, long long> virtual_list;6710		virtual_list.reserve(k * 2);6711		for (long long i = 0; i < k; ++i) {6712			long long j = next<long long>(i, universe - 1);6713			long long vi = virtual_list.count(i) ? virtual_list[i] : i;6714			long long vj = virtual_list.count(j) ? virtual_list[j] : j;6715			virtual_list[j] = vi;6716			virtual_list[i] = vj;6717			res.push_back(virtual_list[i] + left);6718		}6719		std::sort(res.begin(), res.end());6720	}6721	return res;6722}67236724// Valtr-style signed edge components along one axis from n sorted distinct6725// coordinates. The n differences sum to zero.6726inline std::vector<long long>6727valtr_edge_components(const std::vector<long long> &sorted_coords) {6728	int n = sorted_coords.size();6729	std::vector<long long> left, right;6730	left.reserve(n / 2);6731	right.reserve(n / 2);6732	for (int i = 1; i + 1 < n; ++i) {6733		if (next(2) == 0)6734			left.push_back(sorted_coords[i]);6735		else6736			right.push_back(sorted_coords[i]);6737	}6738	long long lo = sorted_coords.front(), hi = sorted_coords.back();6739	std::vector<long long> seq;6740	seq.reserve(n + 1);6741	seq.push_back(lo);6742	for (long long v : left)6743		seq.push_back(v);6744	seq.push_back(hi);6745	for (auto it = right.rbegin(); it != right.rend(); ++it)6746		seq.push_back(*it);6747	seq.push_back(lo);6748	std::vector<long long> comps(n);6749	for (int i = 0; i < n; ++i)6750		comps[i] = seq[i + 1] - seq[i];6751	return comps;6752}67536754// Drops boundary vertices that are collinear with their cyclic neighbors.6755// O(m), m = |points|.6756inline std::vector<point<long long>>6757simplify_strict_boundary(std::vector<point<long long>> points) {6758	int n = points.size();6759	if (n < 3)6760		return points;67616762	std::vector<point<long long>> strict_points;6763	strict_points.reserve(n);6764	for (int i = 0; i < n; ++i) {6765		if (ccw(points[(i + n - 1) % n], points[i], points[(i + 1) % n]) != 0)6766			strict_points.push_back(points[i]);6767	}6768	return strict_points;6769}67706771// Picks k evenly spaced vertices along a longer cyclic boundary.6772// O(k).6773inline std::vector<point<long long>>6774subsample_boundary(const std::vector<point<long long>> &points, int k) {6775	int n = points.size();6776	if (n <= k)6777		return points;67786779	std::vector<point<long long>> sampled_points;6780	sampled_points.reserve(k);6781	for (int i = 0; i < k; ++i)6782		sampled_points.push_back(points[(static_cast<i128>(i) * n) / k]);6783	return sampled_points;6784}67856786// Random translation so the polygon lies in the box.6787// O(|points|).6788inline void place_inside_box(std::vector<point<long long>> &points,6789							 long long min_coord, long long max_coord) {6790	long long width = max_coord - min_coord + 1;67916792	i128 min_x = points[0].x(), max_x = points[0].x();6793	i128 min_y = points[0].y(), max_y = points[0].y();6794	for (const point<long long> &p : points) {6795		min_x = std::min(min_x, static_cast<i128>(p.x()));6796		max_x = std::max(max_x, static_cast<i128>(p.x()));6797		min_y = std::min(min_y, static_cast<i128>(p.y()));6798		max_y = std::max(max_y, static_cast<i128>(p.y()));6799	}68006801	i128 span_x = max_x - min_x;6802	i128 span_y = max_y - min_y;6803	// Random slack keeps the polygon inside the box without filling it.6804	i128 shift_x =6805		min_coord - min_x +6806		next<long long>(0, width - 1 - static_cast<long long>(span_x));6807	i128 shift_y =6808		min_coord - min_y +6809		next<long long>(0, width - 1 - static_cast<long long>(span_y));68106811	for (point<long long> &p : points)6812		p = point<long long>(p.x() + shift_x, p.y() + shift_y);6813}68146815// Random cyclic shift; Valtr walk is already counterclockwise.6816// O(|points|).6817inline void randomize_boundary_order(std::vector<point<long long>> &points) {6818	int rot = next(points.size());6819	if (rot > 0)6820		std::rotate(points.begin(), points.begin() + rot, points.end());6821}68226823// Valtr walk for m edges; bbox minimum translated to the origin.6824// O(m log m).6825inline std::vector<point<long long>>6826valtr_vertices(int m, const std::vector<long long> &x_comp,6827			   std::vector<long long> y_comp) {6828	shuffle(y_comp.begin(), y_comp.end());68296830	std::vector<point<long long>> edges(m);6831	// Upper half-plane (positive y, or y = 0 and x > 0) sorts before lower.6832	auto upper = [](const point<long long> &p) {6833		return p.y() > 0 or (p.y() == 0 and p.x() > 0);6834	};6835	for (int i = 0; i < m; ++i)6836		edges[i] = point<long long>(x_comp[i], y_comp[i]);68376838	std::sort(edges.begin(), edges.end(),6839			  [&upper](const point<long long> &a, const point<long long> &b) {6840				  bool au = upper(a), bu = upper(b);6841				  if (au != bu)6842					  return au;6843				  auto cross = a ^ b;6844				  if (cross != 0)6845					  return cross > 0;6846				  return (a * a) < (b * b);6847			  });68486849	// Prefix-sum the sorted edge vectors to obtain vertex coordinates.6850	i128 cur_x = 0, cur_y = 0;6851	std::vector<i128> px(m), py(m);6852	for (int i = 0; i < m; ++i) {6853		px[i] = cur_x;6854		py[i] = cur_y;6855		cur_x += edges[i].x();6856		cur_y += edges[i].y();6857	}6858	tgen::detail::tgen_ensure_against_bug(6859		cur_x == 0 and cur_y == 0,6860		"geometry: random_convex_polygon: walk did not close");68616862	i128 min_x = px[0], min_y = py[0];6863	for (int i = 1; i < m; ++i) {6864		min_x = std::min(min_x, px[i]);6865		min_y = std::min(min_y, py[i]);6866	}68676868	// Shift so the bbox minimum is at the origin.6869	std::vector<point<long long>> points;6870	points.reserve(m);6871	for (int i = 0; i < m; ++i)6872		points.emplace_back(px[i] - min_x, py[i] - min_y);6873	return points;6874}68756876} // namespace detail68776878// Generates n vertices of a convex integer polygon inside a box.6879// If strict is true, boundary vertices are guaranteed non-collinear when6880// generation succeeds; retry count depends on n and width.6881// Always returns points in counterclockwise order.6882// O(n log n).6883inline std::vector<point<long long>>6884random_convex_polygon(int n, long long min_coord, long long max_coord,6885					  bool strict = false) {6886	tgen_ensure(n >= 3,6887				"geometry: random_convex_polygon: n must be at least 3");6888	tgen_ensure(max_coord >= min_coord,6889				"geometry: random_convex_polygon: min_coord must be at most "6890				"max_coord");6891	tgen_ensure(static_cast<detail::i128>(max_coord) - min_coord + 1 <=6892					std::numeric_limits<long long>::max(),6893				"geometry: random_convex_polygon: coordinate range too large");6894	long long width = max_coord - min_coord + 1;6895	tgen_ensure(6896		width >= n,6897		"geometry: random_convex_polygon: coordinate range too small for n");68986899	// Valtr walk size: n in weak mode; strict mode uses a larger grid so6900	// collinear removal still leaves at least n vertices to subsample.6901	int num_coords = n;6902	if (strict) {6903		// Extra grid lines beyond n: at least 100 (small-n headroom), about6904		// n/1000 for large n, and never more than width - n.6905		int extra = width <= n ? 06906							   : std::min<long long>(std::max(100, n / 1000),6907													 width - n);6908		num_coords = n + extra;6909	}69106911	// Strict mode retries coordinate sampling when simplification leaves < n6912	// vertices; weak mode has no failure path, so one attempt always suffices.6913	const int max_attempts = strict ? 32 : 1;6914	for (int i = 0; i < max_attempts; ++i) {6915		// Build a convex lattice polygon on [0, width - 1]^2, then translate.6916		std::vector<long long> x_sorted =6917			detail::sample_sorted_distinct_in_range(num_coords, 0, width - 1);6918		std::vector<long long> y_sorted =6919			detail::sample_sorted_distinct_in_range(num_coords, 0, width - 1);6920		std::vector<long long> x_comp = detail::valtr_edge_components(x_sorted);6921		std::vector<long long> y_comp = detail::valtr_edge_components(y_sorted);69226923		std::vector<point<long long>> points =6924			detail::valtr_vertices(num_coords, x_comp, std::move(y_comp));69256926		if (strict) {6927			std::vector<point<long long>> simplified =6928				detail::simplify_strict_boundary(std::move(points));6929			// Tight boxes can leave too few vertices -> resample coordinates.6930			if (static_cast<int>(simplified.size()) < n)6931				continue;69326933			points = detail::subsample_boundary(simplified, n);6934		}69356936		detail::place_inside_box(points, min_coord, max_coord);6937		detail::randomize_boundary_order(points);6938		return points;6939	}69406941	// Generation failed.6942	throw tgen::detail::error(6943		"geometry: random_convex_polygon: generation failed: coordinate "6944		"range too small for n");6945}69466947// Random simple polygon through given distinct points.6948// Collinear triples are allowed; fails if all points are collinear.6949// O(n log n) expected if points are "random", O(n^2) worst case.6950inline std::vector<point<long long>> random_simple_polygon_through_points(6951	const std::vector<point<long long>> &points) {6952	int n = points.size();6953	tgen_ensure(n >= 3,6954				"geometry: random_simple_polygon_through_points: need at "6955				"least 3 points");6956	tgen_ensure(6957		static_cast<int>(6958			std::set<point<long long>>(points.begin(), points.end()).size()) ==6959			n,6960		"geometry: random_simple_polygon_through_points: points must "6961		"be distinct");69626963	int idx_a = 0, idx_b = 0;6964	for (int i = 1; i < n; ++i) {6965		if (points[i] < points[idx_a])6966			idx_a = i;6967		if (points[idx_b] < points[i])6968			idx_b = i;6969	}6970	point<long long> A = points[idx_a], B = points[idx_b];69716972	bool all_collinear = true;6973	for (int i = 0; i < n; ++i) {6974		if (i == idx_a or i == idx_b)6975			continue;6976		if (detail::ccw(A, B, points[i]) != 0) {6977			all_collinear = false;6978			break;6979		}6980	}6981	tgen_ensure(!all_collinear,6982				"geometry: random_simple_polygon_through_points: all points "6983				"are collinear; no simple polygon exists");69846985	std::vector<point<long long>> chain;6986	chain.push_back(A);6987	int left_count = 0;6988	for (int i = 0; i < n; ++i) {6989		if (i == idx_a or i == idx_b)6990			continue;6991		if (detail::ccw(A, B, points[i]) > 0) {6992			chain.push_back(points[i]);6993			++left_count;6994		}6995	}6996	chain.push_back(B);6997	for (int i = 0; i < n; ++i) {6998		if (i == idx_a or i == idx_b)6999			continue;7000		if (detail::ccw(A, B, points[i]) <= 0)7001			chain.push_back(points[i]);7002	}7003	chain.push_back(A);70047005	int n1 = 2 + left_count;7006	// Upper chain: A -> B.7007	detail::conquer(chain, 0, n1);7008	// Lower chain: B -> A.7009	detail::conquer(chain, n1 - 1, chain.size());70107011	// Cyclic vertex order: chain[1..n1) then chain[n1..end) (skip each path's7012	// start vertex).7013	std::vector<point<long long>> poly;7014	poly.insert(poly.end(), chain.begin() + 1, chain.begin() + n1);7015	poly.insert(poly.end(), chain.begin() + n1, chain.end());7016	return poly;7017}70187019namespace detail {70207021// Samples n distinct integer points in [min_coord, max_coord]^2.7022// O(n log n).7023inline std::vector<point<long long>>7024random_distinct_points_in_box(int n, long long min_coord, long long max_coord) {7025	long long width = max_coord - min_coord;7026	i128 side_128 = width + 1;7027	i128 universe = side_128 * side_128;7028	tgen_ensure(universe <= std::numeric_limits<long long>::max(),7029				"geometry: random_simple_polygon: coordinate range too large");7030	long long side = side_128;7031	tgen_ensure(universe >= n,7032				"geometry: random_simple_polygon: coordinate range too small "7033				"for n distinct points");70347035	// Decodes a linear grid key (x * side + y) to a point.7036	auto decode = [&](long long key) -> point<long long> {7037		return point<long long>(min_coord + key / side, min_coord + key % side);7038	};70397040	// Repeats until all generated points are not collinear.7041	// Runs O(1) times expected.7042	while (true) {7043		std::vector<long long> keys =7044			distinct_range<long long>(0, universe - 1).gen_list(n).to_std();70457046		std::vector<point<long long>> points;7047		points.reserve(n);7048		for (long long key : keys)7049			points.push_back(decode(key));70507051		// Checks if the points are not all collinear.7052		for (int i = 2; i < n; ++i) {7053			if (ccw(points[0], points[1], points[i]) != 0)7054				return points;7055		}7056	}7057}70587059} // namespace detail70607061// Random simple polygon.7062// If strict, vertex set has no three collinear points7063// (random_points_general_position); otherwise samples distinct grid points7064// (collinear triples allowed), so it might be the case that (polygon[i],7065// polygon[i+1], and polygon[i+2]) are collinear. Polygonizes via7066// random_simple_polygon_through_points.7067// O(n log n) expected.7068inline std::vector<point<long long>>7069random_simple_polygon(int n, long long min_coord, long long max_coord,7070					  bool strict = false) {7071	tgen_ensure(n >= 3,7072				"geometry: random_simple_polygon: n must be at least 3");7073	tgen_ensure(max_coord >= min_coord,7074				"geometry: random_simple_polygon: min_coord must be at most "7075				"max_coord");7076	tgen_ensure(static_cast<detail::i128>(max_coord) - min_coord <=7077					std::numeric_limits<long long>::max(),7078				"geometry: random_simple_polygon: coordinate range too large");70797080	std::vector<point<long long>> points =7081		strict ? random_points_general_position(n, min_coord, max_coord)7082			   : detail::random_distinct_points_in_box(n, min_coord, max_coord);7083	return random_simple_polygon_through_points(points);7084}70857086} // namespace geometry70877088/************7089 *          *7090 *   HACK   *7091 *          *7092 ************/70937094namespace hack {70957096namespace detail {70977098using namespace tgen::detail;70997100// Computes polynomial hash of a string.7101// O(|s|).7102inline int hash_string(const std::string &s, int base, int mod) {7103	long long h = 0;7104	for (char c : s)7105		h = (h * base + c - 'a' + 1) % mod;7106	return h;7107}71087109// Estimates the length of the string to very likely have a collision.7110inline int estimate_length(int alphabet_size, int mod) {7111	// Magic constants.7112	double base_len = 2.5 * std::log(std::sqrt(mod));7113	double scale = std::log(alphabet_size) / std::log(2.0);7114	double adjusted = base_len / std::max(1.0, scale * 0.7);71157116	return static_cast<int>(std::ceil(adjusted));7117}71187119// Collides two strings to have the same polynomial hash.7120// O(sqrt(mod) log(mod)) with high probability.7121inline std::pair<std::string, std::string>7122birthday_attack(const std::vector<std::string> &alphabet, int base, int mod) {7123	tgen_ensure(0 < base and base < mod,7124				"birthday_attack: base must be in (0, mod)");7125	std::map<uint64_t, std::vector<int>> seen;7126	int length = estimate_length(alphabet.size(), mod);71277128	while (true) {7129		std::vector<int> seq(length);71307131		std::string s;71327133		for (int i = 0; i < length; ++i) {7134			seq[i] = next<int>(0, alphabet.size() - 1);7135			s += alphabet[seq[i]];7136		}71377138		int h = hash_string(s, base, mod);71397140		auto it = seen.find(h);7141		if (it != seen.end() and it->second != seq) {7142			std::string a, b;71437144			for (int x : it->second)7145				a += alphabet[x];7146			for (int x : seq)7147				b += alphabet[x];71487149			if (a != b)7150				return {a, b};7151		}71527153		seen[h] = seq;7154	}7155}71567157// Tried to find correct multipliers for unordered_map/set to force7158// collisions. O(1).7159inline std::set<long long> std_hash_multipliers() {7160	std::set<long long> multipliers = {85229};71617162	// Codeforces GCC GNU G++17 7.3.0 case.7163	bool codeforces_gcc_case = true;7164	if (cpp.version_ != 0 and cpp.version_ != 17)7165		codeforces_gcc_case = false;7166	if (compiler.kind_ != compiler_kind::unknown and7167		compiler.kind_ != compiler_kind::gcc)7168		codeforces_gcc_case = false;7169	if (compiler.major_ > 7)7170		codeforces_gcc_case = false;71717172	if (codeforces_gcc_case)7173		multipliers.insert(107897);71747175	return multipliers;7176}71777178} // namespace detail71797180// Fetches prefix of length n of the string "abacabadabacabae...".7181// O(n).7182inline std::string abacaba(int n) {7183	tgen_ensure(n > 0, "str: size must be positive");7184	std::string str = "a";7185	char c = 'a';7186	while (static_cast<int>(str.size()) < n) {7187		int prev_size = str.size();7188		str += ++c;7189		for (int j = 0; j < prev_size and static_cast<int>(str.size()) < n; ++j)7190			str += str[j];7191	}7192	return str;7193}71947195// Two strings that have same polynomial hash for any base, for7196// mod = power of 2 up to 2^64.7197// Thue–Morse.7198// O(1).7199inline std::pair<std::string, std::string> unsigned_polynomial_hash() {7200	std::string a, b;7201	int size = 1 << 10;7202	for (int i = 0; i < size; ++i) {7203		a += 'a' + math::detail::popcount(i) % 2;7204		b += 'a' + ('b' - a[i]);7205	}7206	return {a, b};7207}72087209// Collides two strings to have the same polynomial hash.7210// O(sqrt(mod) log(mod)) with high probability.7211// 0 < base < mod.7212inline std::pair<std::string, std::string> polynomial_hash(int alphabet_size,7213														   int base, int mod) {7214	tgen_ensure(alphabet_size > 1,7215				"hack: polynomial_hash: alphabet size must be greater "7216				"than 1");7217	tgen_ensure(0 < base and base < mod,7218				"hack: polynomial_hash: base must be in (0, mod)");72197220	std::vector<std::string> alphabet(alphabet_size);7221	for (int i = 0; i < alphabet_size; ++i)7222		alphabet[i] = std::string(1, 'a' + i);7223	std::iota(alphabet.begin(), alphabet.end(), 'a');7224	return detail::birthday_attack(alphabet, base, mod);7225}72267227// Collides two strings to have the same polynomial hash for multiple bases7228// and mods (up to 2 pairs).7229// O(sqrt(mod) log^2 (mod)) with high probability,7230// with mod = max(mod_1, mod_2).7231inline std::pair<std::string, std::string>7232polynomial_hash(int alphabet_size, std::vector<int> bases,7233				std::vector<int> mods) {7234	tgen_ensure(bases.size() == mods.size(),7235				"hack: polynomial_hash: bases and mods must have the same "7236				"size");7237	tgen_ensure(bases.size() > 0,7238				"hack: polynomial_hash: must have at least one (base, mod) "7239				"pair");7240	tgen_ensure(bases.size() <= 2,7241				"hack: polynomial_hash: multi-hash hack only supported "7242				"for up to 2 (base, mod) pairs");72437244	std::vector<std::string> alphabet(alphabet_size);7245	for (int i = 0; i < alphabet_size; ++i)7246		alphabet[i] = std::string(1, 'a' + i);7247	auto [S1, T1] = detail::birthday_attack(alphabet, bases[0], mods[0]);7248	if (bases.size() == 1)7249		return {S1, T1};7250	return detail::birthday_attack({S1, T1}, bases[1], mods[1]);7251}72527253// Returns a list of integers for unordered_map/set to force collisions.7254// O(size).7255inline std::vector<long long> std_unordered(int size) {7256	tgen_ensure(size > 0, "hack: std_unordered: size must be positive");7257	std::set<long long> multipliers = detail::std_hash_multipliers();7258	long long mult = 1;7259	std::set<long long>::iterator it = multipliers.begin();72607261	std::vector<long long> list;7262	while (static_cast<int>(list.size()) < size) {7263		list.push_back(mult * (*it));7264		++it;7265		if (it == multipliers.end()) {7266			it = multipliers.begin();7267			++mult;7268		}7269	}7270	return list;7271}72727273// Returns queries that force \Theta(q sqrt n) asymptotic7274// for Mo algorithm for offline range queries.7275// Forces \Theta(q sqrt n) pointer moves for any ordering.7276// O(n log n + q).7277inline std::vector<std::pair<int, int>> mo_worst_case(int n, int q) {7278	std::set<std::pair<int, int>> queries;72797280	// Adversarial case.7281	int sq = std::sqrt(n);7282	for (int i = 0; i < sq; ++i) {7283		for (int j = i; j < sq; ++j) {7284			if (i * sq < n and j * sq < n)7285				queries.emplace(i * sq, j * sq);7286		}7287	}72887289	// Push extra queries.7290	for (int i = 0; i < n; ++i)7291		if (queries.size() < size_t(q)) {7292			queries.emplace(0, i);7293			queries.emplace(i, i);7294			queries.emplace(i, n - 1);7295		}72967297	std::vector<std::pair<int, int>> pool(queries.begin(), queries.end());7298	while (pool.size() < size_t(q)) {7299		int l = next(0, n - 1);7300		pool.emplace_back(l, next(l, n - 1));7301	}73027303	return choose(shuffled(pool), q);7304}73057306// Returns list of strings that have a high cost to insert in a std::set.7307// Forces cost \Theta(size log(size)).7308// Generates: {b, ab, aab, aaab, ...}.7309// O(size log(size)).7310inline std::vector<std::string> string_set_worst_case(int size) {7311	std::vector<std::string> list;7312	int k = 0, left = size;7313	while (left > 0) {7314		int cur_size = std::min(left, k + 1);7315		left -= cur_size;73167317		char right_char = cur_size == k + 1 ? 'b' : 'c';7318		list.push_back(std::string(cur_size - 1, 'a') + right_char);73197320		++k;7321	}7322	return tgen::shuffled(list);7323}73247325// Graph for Dijkstra implementations that relax with <= instead of <.7326// Unit-weight layered graph: 0 -> {1,2}, then disjoint7327// 2x2 gadgets (i,i+1) -> {i+2,i+3} for i = 1,3,5,... Many vertices share the7328// same dist from 0; with `d + w <= dist[j]` each pop re-relaxes the whole7329// frontier below it. m = 2(n - 2) edges.7330// O(n).7331inline egraph<int>::value non_strict_relaxation_dijkstra_bug(int n) {7332	tgen_ensure(7333		n >= 3,7334		"hack: non_strict_relaxation_dijkstra_bug: needs at least 3 vertices");73357336	egraph<int>::value g(n, {}, true);7337	g.edge_weighted();7338	g.add_edge(0, 1, 1);7339	g.add_edge(0, 2, 1);7340	for (int i = 1; i + 2 < n; i += 2) {7341		g.add_edge(i, i + 2, 1);7342		if (i + 3 < n)7343			g.add_edge(i, i + 3, 1);73447345		g.add_edge(i + 1, i + 2, 1);7346		if (i + 3 < n)7347			g.add_edge(i + 1, i + 3, 1);7348	}73497350	return g.shuffle_except({0});7351}73527353// Graph for Dijkstra implementations that do not skip stale heap entries7354// (`if (d > dist[i]) continue`).7355// Hub mid = n/2: star 0 -> 1..mid-1 (weights 1..mid-1), funnel i -> mid7356// (weights 1,3,5,...), then mid -> mid+1.. (weight 1).7357// Without a stale-heap check, mid and its in-neighbors are re-popped and7358// re-relax.7359// m = n + mid - 3 edges, mid = floor(n/2).7360// O(n).7361inline egraph<int>::value stale_heap_dijkstra_bug(int n) {7362	tgen_ensure(n >= 4,7363				"hack: stale_heap_dijkstra_bug: needs at least 4 vertices");73647365	int mid = n / 2;7366	egraph<int>::value g(n, {}, true);7367	g.edge_weighted();7368	for (int i = 1; i < mid; ++i)7369		g.add_edge(0, i, i);7370	for (int i = 1; i < mid; ++i)7371		g.add_edge(i, mid, 2 * (mid - i) - 1);7372	for (int i = mid + 1; i < n; ++i)7373		g.add_edge(mid, i, 1);73747375	return g.shuffle_except({0});7376}73777378// Worst-case for FIFO-SPFA.7379// Forces Omega(n^2) from vertex 0 (Theta(n*m), m = 2n - 3).7380// Upper chain ai -> a(i+1) weight 1; lower chain bi -> b(i+1) weight 0;7381// vertical ai -> bi weight 0; cross bi -> a(i+1) weight 1. Upper chain sets7382// loose dist first; cross edges from settled bi then improve a(i+1).7383// m = 2n - 3.7384// O(n).7385inline egraph<int>::value spfa(int n) {7386	tgen_ensure(n >= 2, "hack: spfa: n must be at least 2");7387	tgen_ensure(n % 2 == 0, "hack: spfa: n must be even");73887389	egraph<int>::value g(n, {}, true);7390	g.edge_weighted();73917392	const int k = n / 2;7393	for (int i = 0; i + 1 < k; ++i)7394		g.add_edge(i, i + 1, 1);7395	for (int i = 0; i + 1 < k; ++i)7396		g.add_edge(k + i, k + i + 1, 0);7397	for (int i = 0; i < k; ++i)7398		g.add_edge(i, k + i, 0);7399	for (int i = 0; i + 1 < k; ++i)7400		g.add_edge(k + i, i + 1, 1);74017402	return g.shuffle_except({0});7403}74047405// Zadeh (1972) anti-shortest-paths flow network for Edmonds-Karp and Dinitz.7406// Source is vertex 0; sink is vertex 4l + 2k + 1.7407// n = 4l + 2k + 2, m = 6l + 4k + k^2 - 4.7408// O(l + k^2).7409inline egraph<int>::value dinitz_worst_case(int k, int l) {7410	tgen_ensure(k >= 1, "hack: dinitz_worst_case: k must be at least 1");7411	tgen_ensure(l >= 1, "hack: dinitz_worst_case: l must be at least 1");74127413	const int p1 = 2 * l - 1;7414	const int p2 = 2 * l;7415	const int q1 = 2 * l + 1;7416	const int q2 = 2 * l + 2;7417	const int n = 4 * l + 2 * k + 2;74187419	const int flow_cap = k * k * l;7420	const int layer_cap = k * k;74217422	auto a = [&](int i) { return 2 * l + 3 + 2 * i; };7423	auto b = [&](int i) { return 2 * l + 4 + 2 * i; };7424	auto t = [&](int i) { return 4 * l + 2 * k + 1 - i; };74257426	egraph<int>::value g(n, {}, true);7427	g.edge_weighted();74287429	for (int i = 0; i + 1 < 2 * l - 1; ++i)7430		g.add_edge(i, i + 1, flow_cap);7431	for (int i = 0; i + 1 < 2 * l - 1; ++i)7432		g.add_edge(t(i + 1), t(i), flow_cap);74337434	for (int i = 0; i < 2 * l - 1; i += 2) {7435		g.add_edge(i, i % 4 == 0 ? p1 : p2, layer_cap);7436		g.add_edge(i % 4 == 0 ? q1 : q2, t(i), layer_cap);7437	}74387439	for (int i = 0; i < k; ++i) {7440		g.add_edge(p1, a(i), flow_cap);7441		g.add_edge(p2, b(i), flow_cap);7442		g.add_edge(a(i), q2, flow_cap);7443		g.add_edge(b(i), q1, flow_cap);7444	}74457446	for (int i = 0; i < k; ++i)7447		for (int j = 0; j < k; ++j)7448			g.add_edge(a(i), b(j), 1);74497450	return g;7451}74527453// Returns a mask of length 19938, with weights such that xor-ing with mt199377454// outputs yields 0.7455// O(1).7456template <typename T> std::vector<bool> mt19937_xor_hash() {7457	static_assert(std::is_same_v<T, int> or std::is_same_v<T, long long>,7458				  "hack: mt19937_xor_hash: T must be int or long long");74597460	constexpr std::size_t deg = 19937;74617462	std::bitset<deg + 1> a, b, c;7463	b[deg] = c[deg] = 1;7464	std::size_t l = 0, shift = 1;7465	std::mt19937 rng32;7466	std::mt19937_64 rng64;7467	for (std::size_t n = 0; n < deg * 2; ++n) {7468		a >>= 1;7469		if constexpr (std::is_same_v<T, int>)7470			a[deg] = rng32() & 1;7471		else7472			a[deg] = rng64() & 1;74737474		if ((c & a).count() % 2 == 0) {7475			++shift;7476			continue;7477		}74787479		std::bitset<deg + 1> oc = c;7480		c ^= (b >> shift);7481		if (2 * l <= n) {7482			l = n + 1 - l;7483			b = oc;7484			shift = 1;7485		} else {7486			++shift;7487		}7488	}74897490	std::vector<bool> mask(deg + 1);7491	for (std::size_t i = 0; i <= deg; ++i)7492		mask[i] = c[i];7493	return mask;7494}74957496// Convex polygon that breaks naive rotating calipers for maximum vertex7497// distance (advances j while dist(i, next(j)) > dist(i, j) instead of using7498// ccw).7499// O(1).7500inline std::vector<geometry::point<double>>7501naive_rotating_calipers_max_dist_bug() {7502	return {7503		{-0.9846, -1.53251}, {0.49946, 1.19525},  {0.79916, 0.98291},7504		{4.02136, -1.57843}, {3.92734, -2.37856}, {3.88558, -2.37188},7505	};7506}75077508namespace detail {75097510// Builds a hack block of order k (length fib(2k+1)).7511// O(fib(2k+1)).7512inline std::vector<int> segment_tree_beats_worst_case_block(int k) {7513	tgen_ensure(k >= 1,7514				"hack: segment_tree_beats_worst_case: k must be at least 1");75157516	std::vector<int> a(k + 1), b(k + 1);7517	std::vector<std::vector<int>> vf(k + 1), vg(k + 1);75187519	a[1] = b[1] = 1;7520	vf[1] = {1};7521	vg[1] = {1, 0};75227523	for (int i = 2; i <= k; ++i) {7524		b[i] = b[i - 1] + a[i - 1];7525		a[i] = b[i] + a[i - 1];7526		for (int x : vf[i - 1])7527			vf[i].push_back(x + a[i] + b[i]);7528		vf[i].push_back(a[i]);7529		for (int x : vg[i - 1])7530			vf[i].push_back(x + a[i]);7531		vg[i] = vf[i];7532		vg[i].push_back(0);7533		for (int x : vg[i - 1])7534			vg[i].push_back(x);7535	}75367537	vf[k].push_back(0);7538	return vf[k];7539}75407541// Appends one update round for the tiled array (offset (round * an) mod L).7542// O(fib(2k+1)).7543inline void7544segment_tree_beats_append_round(std::vector<std::vector<int>> &updates,7545								int block_len, int an, int bn, int n,7546								int round) {7547	const int off = (round * an) % block_len;7548	const int add_off = (off + block_len - bn) % block_len;7549	for (int k = 0; k < block_len; ++k) {7550		const int s = k * block_len * block_len;7551		const int sub_end = off + an;7552		if (sub_end <= block_len)7553			updates.push_back({1, s + off, s + sub_end, bn});7554		else {7555			updates.push_back({1, s + off, s + block_len, bn});7556			updates.push_back({1, s, s + (sub_end - block_len), bn});7557		}7558		const int add_end = add_off + bn;7559		if (add_end <= block_len)7560			updates.push_back({0, s + add_off, s + add_end, an});7561		else {7562			updates.push_back({0, s + add_off, s + block_len, an});7563			updates.push_back({0, s, s + (add_end - block_len), an});7564		}7565	}7566	updates.push_back({2, 0, n, an});7567	for (int k = 0; k < block_len; ++k) {7568		const int s = k * block_len * block_len;7569		updates.push_back({3, s + (off + an - 1) % block_len, 0});7570	}7571}75727573} // namespace detail75747575// Array and updates for worst case of segment tree beats.7576// O(fib(2k+1)^3 + q).7577inline std::pair<std::vector<int>, std::vector<std::vector<int>>>7578segment_tree_beats_worst_case(int k, int q) {7579	tgen_ensure(k >= 1,7580				"hack: segment_tree_beats_worst_case: k must be at least 1");7581	tgen_ensure(k <= 7, "hack: segment_tree_beats_worst_case: k too large");7582	tgen_ensure(q > 0,7583				"hack: segment_tree_beats_worst_case: q must be positive");75847585	const auto &fib = math::fibonacci();7586	const int block_len = fib[k * 2 + 1];7587	const int an = fib[k * 2];7588	const int bn = fib[k * 2 - 1];75897590	const int len = block_len;7591	const int total = len * len * len;75927593	std::vector<int> block = detail::segment_tree_beats_worst_case_block(k);7594	std::vector<int> arr(total, 0);7595	for (int x = 0; x < block_len; ++x) {7596		const int s = x * len * len;7597		for (int i = 0; i < block_len; ++i)7598			arr[s + i] = block[i];7599	}76007601	std::vector<std::vector<int>> updates;7602	updates.reserve(q);7603	const int n = total;7604	for (int round = 0; updates.size() < static_cast<std::size_t>(q); ++round) {7605		detail::segment_tree_beats_append_round(updates, block_len, an, bn, n,7606												round);7607		if (updates.size() > static_cast<std::size_t>(q))7608			updates.resize(q);7609	}7610	return {arr, updates};7611}76127613} // namespace hack76147615/*********************7616 *                   *7617 *   MISCELLANEOUS   *7618 *                   *7619 *********************/76207621namespace misc {76227623// Generates a uniformly random balanced parentheses sequence with k '(' and k7624// ')'. Valid means that for no prefix there are more ')' than '('.7625// O(size).7626inline std::string gen_parenthesis(int size) {7627	tgen_ensure(size > 0 and size % 2 == 0,7628				"misc: parenthesis: size must be a positive even number");76297630	int k = size / 2;7631	std::string s;7632	int open = 0, close = 0;76337634	for (int i = 0; i < size; ++i) {7635		if (open == k) {7636			s += ')';7637			++close;7638			continue;7639		}7640		if (open == close) {7641			s += '(';7642			++open;7643			continue;7644		}76457646		long long a = k - open, b = k - close, h = open - close;76477648		// Probability of placing '(':7649		// P('(') = (k - open) * (h + 2) / ((k - open + k - close) * (h + 1))7650		// Derived from ballot numbers ratio.7651		long long num = a * (h + 2);7652		long long den = (a + b) * (h + 1);76537654		if (next<long long>(1, den) <= num) {7655			s += '(';7656			++open;7657		} else {7658			s += ')';7659			++close;7660		}7661	}76627663	return s;7664}76657666} // namespace misc76677668} // namespace tgen