ライブラリ一式
| 投稿者 | manta1130 |
|---|---|
| 言語 | Rust |
| 更新日時 (JST) | 2026-11-03 08:38:24 |
| サイズ | 182,239 Byte |
| 参照URL | |
日時はサーバーが記録したものです。更新すると更新日時が変わるので、不都合がある場合は新しく投稿してください。
#[allow(unused_imports)]
use proconio::{
fastout, input,
marker::{Bytes, Chars, Isize1, Usize1},
};
fn main() {
println!("ライブラリ用のダミーコードです。");
}
//https://github.com/rust-lang-ja/ac-library-rs
//https://github.com/manta1130/competitive-template-rs
use complex::*;
use convolution::*;
use counter::*;
use dlp::*;
use dpfactorial::*;
use dsu::*;
use dsu_potential::*;
use eratosthenes::*;
use fenwicktree::*;
use fft::*;
use geometry::*;
use global::*;
use graph::*;
use graph_search::*;
use input::*;
use lazysegtree::*;
use lca::*;
use math::*;
use matrix::*;
use maxflow::*;
use mex::*;
use mincostflow::*;
use modint::*;
use ordfloat::*;
use pcg64::*;
use primenumber::*;
use raq_bit::*;
use rng_trait::*;
use rollinghash::*;
use scc::*;
use segtree::*;
use segtree_template::*;
use string::*;
use timer::*;
use twosat::*;
use vectools::*;
use xorshift::*;
pub mod complex {
#![warn(missing_docs)]
use std::fmt;
use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};
#[derive(Clone, PartialEq, Debug, Eq, Copy, Default)]
pub struct Complex<T> {
pub re: T,
pub im: T,
}
impl<T> Complex<T> {
#[must_use]
pub const fn new(re: T, im: T) -> Self {
Self { re, im }
}
}
impl<T: Copy + Neg<Output = T>> Complex<T> {
#[must_use]
pub fn conj(self) -> Self {
Self::new(self.re, -self.im)
}
}
impl<T> Complex<T>
where
T: Copy + Add<Output = T> + Mul<Output = T>,
{
#[must_use]
pub fn norm_sqr(self) -> T {
self.re * self.re + self.im * self.im
}
}
impl<T> Complex<T>
where
T: Copy + Add<Output = T> + Mul<Output = T> + Into<f64>,
{
#[must_use]
pub fn norm(self) -> f64 {
self.norm_sqr().into().sqrt()
}
#[must_use]
pub fn arg(self) -> f64 {
self.im.into().atan2(self.re.into())
}
}
impl<T> From<(T, T)> for Complex<T> {
fn from((re, im): (T, T)) -> Self {
Self::new(re, im)
}
}
impl<T: fmt::Display> fmt::Display for Complex<T> {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(f, "{}{:+}i", self.re, self.im)
}
}
impl<T: Add<Output = T>> Add for Complex<T> {
type Output = Self;
fn add(self, rhs: Self) -> Self {
Self::new(self.re + rhs.re, self.im + rhs.im)
}
}
impl<T: AddAssign> AddAssign for Complex<T> {
fn add_assign(&mut self, rhs: Self) {
self.re += rhs.re;
self.im += rhs.im;
}
}
impl<T: Sub<Output = T>> Sub for Complex<T> {
type Output = Self;
fn sub(self, rhs: Self) -> Self {
Self::new(self.re - rhs.re, self.im - rhs.im)
}
}
impl<T: SubAssign> SubAssign for Complex<T> {
fn sub_assign(&mut self, rhs: Self) {
self.re -= rhs.re;
self.im -= rhs.im;
}
}
impl<T> Mul for Complex<T>
where
T: Copy + Add<Output = T> + Sub<Output = T> + Mul<Output = T>,
{
type Output = Self;
fn mul(self, rhs: Self) -> Self {
Self::new(
self.re * rhs.re - self.im * rhs.im,
self.re * rhs.im + self.im * rhs.re,
)
}
}
impl<T> MulAssign for Complex<T>
where
T: Copy + Add<Output = T> + Sub<Output = T> + Mul<Output = T>,
{
fn mul_assign(&mut self, rhs: Self) {
*self = *self * rhs;
}
}
impl<T> Div for Complex<T>
where
T: Copy
+ Add<Output = T>
+ Sub<Output = T>
+ Mul<Output = T>
+ Div<Output = T>
+ Neg<Output = T>,
{
type Output = Self;
fn div(self, rhs: Self) -> Self {
let denominator = rhs.norm_sqr();
let numerator = self * rhs.conj();
Self::new(numerator.re / denominator, numerator.im / denominator)
}
}
impl<T> DivAssign for Complex<T>
where
T: Copy
+ Add<Output = T>
+ Sub<Output = T>
+ Mul<Output = T>
+ Div<Output = T>
+ Neg<Output = T>,
{
fn div_assign(&mut self, rhs: Self) {
*self = *self / rhs;
}
}
}
pub mod convolution {
macro_rules! modulus {
($($name:ident),*) => {
$(
#[derive(Copy, Clone, Eq, PartialEq)]
enum $name {}
impl Modulus for $name {
const VALUE: u32 = $name as _;
const HINT_VALUE_IS_PRIME: bool = true;
fn butterfly_cache() -> &'static ::std::thread::LocalKey<::std::cell::RefCell<::std::option::Option<$crate::modint::ButterflyCache<Self>>>> {
thread_local! {
static BUTTERFLY_CACHE: ::std::cell::RefCell<::std::option::Option<$crate::modint::ButterflyCache<$name>>> = ::std::default::Default::default();
}
&BUTTERFLY_CACHE
}
}
)*
};
}
use crate::{
internal_bit, internal_math,
modint::{ButterflyCache, Modulus, RemEuclidU32, StaticModInt},
};
use std::{
cmp,
convert::{TryFrom, TryInto as _},
fmt,
};
#[allow(clippy::many_single_char_names)]
pub fn convolution<M>(a: &[StaticModInt<M>], b: &[StaticModInt<M>]) -> Vec<StaticModInt<M>>
where
M: Modulus,
{
if a.is_empty() || b.is_empty() {
return vec![];
}
let (n, m) = (a.len(), b.len());
if cmp::min(n, m) <= 60 {
let (n, m, a, b) = if n < m { (m, n, b, a) } else { (n, m, a, b) };
let mut ans = vec![StaticModInt::new(0); n + m - 1];
for i in 0..n {
for j in 0..m {
ans[i + j] += a[i] * b[j];
}
}
return ans;
}
let (mut a, mut b) = (a.to_owned(), b.to_owned());
let z = 1 << internal_bit::ceil_pow2((n + m - 1) as _);
a.resize(z, StaticModInt::raw(0));
butterfly(&mut a);
b.resize(z, StaticModInt::raw(0));
butterfly(&mut b);
for (a, b) in a.iter_mut().zip(&b) {
*a *= b;
}
butterfly_inv(&mut a);
a.resize(n + m - 1, StaticModInt::raw(0));
let iz = StaticModInt::new(z).inv();
for a in &mut a {
*a *= iz;
}
a
}
pub fn convolution_raw<T, M>(a: &[T], b: &[T]) -> Vec<T>
where
T: RemEuclidU32 + TryFrom<u32> + Clone,
T::Error: fmt::Debug,
M: Modulus,
{
let a = a.iter().cloned().map(Into::into).collect::<Vec<_>>();
let b = b.iter().cloned().map(Into::into).collect::<Vec<_>>();
convolution::<M>(&a, &b)
.into_iter()
.map(|z| {
z.val()
.try_into()
.expect("the numeric type is smaller than the modulus")
})
.collect()
}
#[allow(clippy::many_single_char_names)]
pub fn convolution_i64(a: &[i64], b: &[i64]) -> Vec<i64> {
const M1: u64 = 754_974_721; // 2^24
const M2: u64 = 167_772_161; // 2^25
const M3: u64 = 469_762_049; // 2^26
const M2M3: u64 = M2 * M3;
const M1M3: u64 = M1 * M3;
const M1M2: u64 = M1 * M2;
const M1M2M3: u64 = M1M2.wrapping_mul(M3);
modulus!(M1, M2, M3);
if a.is_empty() || b.is_empty() {
return vec![];
}
let (_, i1) = internal_math::inv_gcd(M2M3 as _, M1 as _);
let (_, i2) = internal_math::inv_gcd(M1M3 as _, M2 as _);
let (_, i3) = internal_math::inv_gcd(M1M2 as _, M3 as _);
let c1 = convolution_raw::<i64, M1>(a, b);
let c2 = convolution_raw::<i64, M2>(a, b);
let c3 = convolution_raw::<i64, M3>(a, b);
c1.into_iter()
.zip(c2)
.zip(c3)
.map(|((c1, c2), c3)| {
const OFFSET: &[u64] = &[0, 0, M1M2M3, 2 * M1M2M3, 3 * M1M2M3];
let mut x = [(c1, i1, M1, M2M3), (c2, i2, M2, M1M3), (c3, i3, M3, M1M2)]
.iter()
.map(|&(c, i, m1, m2)| {
c.wrapping_mul(i).rem_euclid(m1 as _).wrapping_mul(m2 as _)
})
.fold(0, i64::wrapping_add);
let mut diff = c1 - internal_math::safe_mod(x, M1 as _);
if diff < 0 {
diff += M1 as i64;
}
x = x.wrapping_sub(OFFSET[diff.rem_euclid(5) as usize] as _);
x
})
.collect()
}
#[allow(clippy::many_single_char_names)]
fn butterfly<M: Modulus>(a: &mut [StaticModInt<M>]) {
let n = a.len();
let h = internal_bit::ceil_pow2(n as u32);
M::butterfly_cache().with(|cache| {
let mut cache = cache.borrow_mut();
let ButterflyCache { sum_e, .. } = cache.get_or_insert_with(prepare);
for ph in 1..=h {
let w = 1 << (ph - 1);
let p = 1 << (h - ph);
let mut now = StaticModInt::<M>::new(1);
for s in 0..w {
let offset = s << (h - ph + 1);
for i in 0..p {
let l = a[i + offset];
let r = a[i + offset + p] * now;
a[i + offset] = l + r;
a[i + offset + p] = l - r;
}
now *= sum_e[(!s).trailing_zeros() as usize];
}
}
});
}
#[allow(clippy::many_single_char_names)]
fn butterfly_inv<M: Modulus>(a: &mut [StaticModInt<M>]) {
let n = a.len();
let h = internal_bit::ceil_pow2(n as u32);
M::butterfly_cache().with(|cache| {
let mut cache = cache.borrow_mut();
let ButterflyCache { sum_ie, .. } = cache.get_or_insert_with(prepare);
for ph in (1..=h).rev() {
let w = 1 << (ph - 1);
let p = 1 << (h - ph);
let mut inow = StaticModInt::<M>::new(1);
for s in 0..w {
let offset = s << (h - ph + 1);
for i in 0..p {
let l = a[i + offset];
let r = a[i + offset + p];
a[i + offset] = l + r;
a[i + offset + p] = StaticModInt::new(M::VALUE + l.val() - r.val()) * inow;
}
inow *= sum_ie[(!s).trailing_zeros() as usize];
}
}
});
}
fn prepare<M: Modulus>() -> ButterflyCache<M> {
let g = StaticModInt::<M>::raw(internal_math::primitive_root(M::VALUE as i32) as u32);
let mut es = [StaticModInt::<M>::raw(0); 30]; // es[i]^(2^(2+i)) == 1
let mut ies = [StaticModInt::<M>::raw(0); 30];
let cnt2 = (M::VALUE - 1).trailing_zeros() as usize;
let mut e = g.pow(((M::VALUE - 1) >> cnt2).into());
let mut ie = e.inv();
for i in (2..=cnt2).rev() {
es[i - 2] = e;
ies[i - 2] = ie;
e *= e;
ie *= ie;
}
let sum_e = es
.iter()
.scan(StaticModInt::new(1), |acc, e| {
*acc *= e;
Some(*acc)
})
.collect();
let sum_ie = ies
.iter()
.scan(StaticModInt::new(1), |acc, ie| {
*acc *= ie;
Some(*acc)
})
.collect();
ButterflyCache { sum_e, sum_ie }
}
}
pub mod counter {
#![warn(missing_docs)]
use std::collections::HashMap;
use std::hash::Hash;
pub trait Counter: Iterator + Sized {
fn counter(self) -> HashMap<Self::Item, usize>
where
Self::Item: Eq + Hash,
{
let mut result = HashMap::new();
for item in self {
*result.entry(item).or_default() += 1;
}
result
}
}
impl<I: Iterator> Counter for I {}
}
pub mod dlp {
#![warn(missing_docs)]
use std::collections::HashMap;
fn gcd(mut a: u64, mut b: u64) -> u64 {
while b != 0 {
(a, b) = (b, a % b);
}
a
}
fn mul_mod(a: u64, b: u64, modulus: u64) -> u64 {
(u128::from(a) * u128::from(b) % u128::from(modulus)) as u64
}
fn pow_mod(mut base: u64, mut exponent: u64, modulus: u64) -> u64 {
let mut result = 1 % modulus;
while exponent != 0 {
if exponent & 1 != 0 {
result = mul_mod(result, base, modulus);
}
base = mul_mod(base, base, modulus);
exponent >>= 1;
}
result
}
fn inverse_mod(value: u64, modulus: u64) -> u64 {
fn egcd(a: i128, b: i128) -> (i128, i128, i128) {
if b == 0 {
(a, 1, 0)
} else {
let (g, x, y) = egcd(b, a % b);
(g, y, x - (a / b) * y)
}
}
let (g, x, _) = egcd(i128::from(value), i128::from(modulus));
assert_eq!(g, 1);
x.rem_euclid(i128::from(modulus)) as u64
}
fn bsgs_coprime(base: u64, target: u64, modulus: u64, positive: bool) -> Option<u64> {
if gcd(target, modulus) != 1 {
return None;
}
let block = modulus.isqrt().saturating_add(1);
let mut babies = HashMap::new();
let (first, last) = if positive { (1, block) } else { (0, block - 1) };
let mut value = pow_mod(base, first, modulus);
for exponent in first..=last {
if value == target {
return Some(exponent);
}
if babies.insert(value, exponent).is_some() {
return None;
}
value = mul_mod(value, base, modulus);
}
let giant = pow_mod(inverse_mod(base, modulus), block, modulus);
let mut current = target;
for i in 0..=block {
if let Some(&j) = babies.get(¤t) {
return i.checked_mul(block)?.checked_add(j);
}
current = mul_mod(current, giant, modulus);
}
None
}
#[must_use]
pub fn dlp(base: u64, target: u64, modulus: u64, allow_zero: bool) -> Option<u64> {
assert_ne!(modulus, 0, "the modulus must be positive");
if modulus == 1 {
return Some(u64::from(!allow_zero));
}
let base = base % modulus;
let mut target = target % modulus;
let mut reduced_modulus = modulus;
let mut accumulated = 1_u64;
let mut removed = 0_u64;
loop {
let divisor = gcd(base, reduced_modulus);
if divisor == 1 {
break;
}
if target == accumulated && (allow_zero || removed > 0) {
return Some(removed);
}
if !target.is_multiple_of(divisor) {
return None;
}
target /= divisor;
reduced_modulus /= divisor;
accumulated = mul_mod(accumulated, base / divisor, reduced_modulus);
removed += 1;
}
let adjusted = mul_mod(
target,
inverse_mod(accumulated, reduced_modulus),
reduced_modulus,
);
let positive = !allow_zero && removed == 0;
bsgs_coprime(base % reduced_modulus, adjusted, reduced_modulus, positive)
.and_then(|x| x.checked_add(removed))
}
}
pub mod dpfactorial {
#![warn(missing_docs)]
use crate::modint::{ModInt, ModIntBase};
#[derive(Clone, Debug)]
pub struct DpFactorial<M: ModIntBase = ModInt> {
factorials: Vec<M>,
inverse_factorials: Vec<M>,
modulus: u32,
}
impl<M: ModIntBase> Default for DpFactorial<M> {
fn default() -> Self {
Self::new()
}
}
impl<M: ModIntBase> DpFactorial<M> {
#[must_use]
pub fn new() -> Self {
Self {
factorials: vec![M::new(1_u8)],
inverse_factorials: vec![M::new(1_u8)],
modulus: M::modulus(),
}
}
fn ensure(&mut self, n: usize) {
assert_eq!(
M::modulus(),
self.modulus,
"the modulus changed after DpFactorial construction"
);
assert!(
n < self.modulus as usize,
"factorial inverse is not supported at or above the modulus"
);
if n < self.factorials.len() {
return;
}
let old_len = self.factorials.len();
self.factorials.reserve(n + 1 - old_len);
self.inverse_factorials.reserve(n + 1 - old_len);
for i in old_len..=n {
let inverse = if i == 1 {
M::new(1_u8)
} else {
let remainder = self.modulus as usize % i;
assert_ne!(remainder, 0, "factorial inverse does not exist");
-M::new(self.modulus as usize / i)
* self.factorials[remainder - 1]
* self.inverse_factorials[remainder]
};
let factorial = self.factorials[i - 1] * M::new(i);
let inverse_factorial = self.inverse_factorials[i - 1] * inverse;
self.factorials.push(factorial);
self.inverse_factorials.push(inverse_factorial);
}
}
pub fn factorial(&mut self, n: usize) -> M {
self.ensure(n);
self.factorials[n]
}
pub fn factorial_inverse(&mut self, n: usize) -> M {
self.ensure(n);
self.inverse_factorials[n]
}
pub fn combination(&mut self, n: usize, r: usize) -> M {
if r > n {
return M::new(0_u8);
}
self.ensure(n);
self.factorials[n] * self.inverse_factorials[r] * self.inverse_factorials[n - r]
}
pub fn combination_small(&mut self, n: usize, r: usize) -> M {
if r > n {
return M::new(0_u8);
}
let r = r.min(n - r);
self.ensure(r);
if r == 0 {
return M::new(1_u8);
}
((n - r + 1)..=n).fold(self.inverse_factorials[r], |acc, x| acc * M::new(x))
}
pub fn permutation(&mut self, n: usize, r: usize) -> M {
if r > n {
return M::new(0_u8);
}
self.ensure(n);
self.factorials[n] * self.inverse_factorials[n - r]
}
}
}
pub mod dsu {
#[derive(Clone, Debug)]
pub struct Dsu {
n: usize,
parent_or_size: Vec<i32>,
}
impl Dsu {
pub fn new(size: usize) -> Self {
Self {
n: size,
parent_or_size: vec![-1; size],
}
}
pub fn merge(&mut self, a: usize, b: usize) -> usize {
assert!(a < self.n);
assert!(b < self.n);
let (mut x, mut y) = (self.leader(a), self.leader(b));
if x == y {
return x;
}
if -self.parent_or_size[x] < -self.parent_or_size[y] {
std::mem::swap(&mut x, &mut y);
}
self.parent_or_size[x] += self.parent_or_size[y];
self.parent_or_size[y] = x as i32;
x
}
pub fn same(&mut self, a: usize, b: usize) -> bool {
assert!(a < self.n);
assert!(b < self.n);
self.leader(a) == self.leader(b)
}
pub fn leader(&mut self, a: usize) -> usize {
assert!(a < self.n);
if self.parent_or_size[a] < 0 {
return a;
}
self.parent_or_size[a] = self.leader(self.parent_or_size[a] as usize) as i32;
self.parent_or_size[a] as usize
}
pub fn size(&mut self, a: usize) -> usize {
assert!(a < self.n);
let x = self.leader(a);
-self.parent_or_size[x] as usize
}
pub fn groups(&mut self) -> Vec<Vec<usize>> {
let mut leader_buf = vec![0; self.n];
let mut group_size = vec![0; self.n];
for i in 0..self.n {
leader_buf[i] = self.leader(i);
group_size[leader_buf[i]] += 1;
}
let mut result = vec![Vec::new(); self.n];
for i in 0..self.n {
result[i].reserve(group_size[i]);
}
for i in 0..self.n {
result[leader_buf[i]].push(i);
}
result
.into_iter()
.filter(|x| !x.is_empty())
.collect::<Vec<Vec<usize>>>()
}
}
}
pub mod dsu_potential {
#![warn(missing_docs)]
#[derive(Copy, Clone, Debug, Eq, PartialEq)]
pub enum PotentialMerge {
Merged(usize),
AlreadyConsistent(usize),
Contradiction(usize),
}
#[derive(Clone, Debug)]
pub struct DsuPotential {
parent_or_size: Vec<isize>,
weight_to_parent: Vec<i64>,
}
impl DsuPotential {
#[must_use]
pub fn new(size: usize) -> Self {
Self {
parent_or_size: vec![-1; size],
weight_to_parent: vec![0; size],
}
}
#[must_use]
pub fn len(&self) -> usize {
self.parent_or_size.len()
}
#[must_use]
pub fn is_empty(&self) -> bool {
self.len() == 0
}
pub fn merge(&mut self, a: usize, b: usize, difference: i64) -> PotentialMerge {
self.check(a);
self.check(b);
let pa = self.potential(a);
let pb = self.potential(b);
let (mut x, mut y) = (self.leader(a), self.leader(b));
let mut root_difference = difference
.checked_add(pa)
.and_then(|v| v.checked_sub(pb))
.expect("potential overflow");
if x == y {
return if root_difference == 0 {
PotentialMerge::AlreadyConsistent(x)
} else {
PotentialMerge::Contradiction(x)
};
}
if -self.parent_or_size[x] < -self.parent_or_size[y] {
std::mem::swap(&mut x, &mut y);
root_difference = root_difference.checked_neg().expect("potential overflow");
}
self.parent_or_size[x] += self.parent_or_size[y];
self.parent_or_size[y] = x as isize;
self.weight_to_parent[y] = root_difference;
PotentialMerge::Merged(x)
}
pub fn same(&mut self, a: usize, b: usize) -> bool {
self.check(a);
self.check(b);
self.leader(a) == self.leader(b)
}
pub fn leader(&mut self, vertex: usize) -> usize {
self.check(vertex);
let mut root = vertex;
let mut path = Vec::new();
while self.parent_or_size[root] >= 0 {
path.push(root);
root = self.parent_or_size[root] as usize;
}
let mut total = 0_i64;
for node in path.into_iter().rev() {
total = total
.checked_add(self.weight_to_parent[node])
.expect("potential overflow");
self.parent_or_size[node] = root as isize;
self.weight_to_parent[node] = total;
}
root
}
pub fn potential(&mut self, vertex: usize) -> i64 {
self.leader(vertex);
self.weight_to_parent[vertex]
}
pub fn diff(&mut self, a: usize, b: usize) -> Option<i64> {
if !self.same(a, b) {
return None;
}
self.potential(b).checked_sub(self.potential(a))
}
pub fn size(&mut self, vertex: usize) -> usize {
let root = self.leader(vertex);
(-self.parent_or_size[root]) as usize
}
pub fn groups(&mut self) -> Vec<Vec<usize>> {
let mut group_indices = vec![usize::MAX; self.len()];
let mut groups = Vec::<Vec<usize>>::new();
for v in 0..self.len() {
let root = self.leader(v);
let index = if group_indices[root] == usize::MAX {
group_indices[root] = groups.len();
groups.push(Vec::new());
groups.len() - 1
} else {
group_indices[root]
};
groups[index].push(v);
}
groups
}
fn check(&self, v: usize) {
assert!(v < self.len(), "vertex index out of bounds");
}
}
}
pub mod eratosthenes {
#![warn(missing_docs)]
use std::convert::TryFrom;
use std::error::Error;
use std::fmt;
#[derive(Clone, Debug)]
pub struct EratosthenesSieve {
smallest: Vec<usize>,
primes: Vec<usize>,
}
impl EratosthenesSieve {
#[must_use]
pub fn new(limit: usize) -> Self {
let mut smallest = vec![0; limit.checked_add(1).expect("sieve limit overflow")];
let mut primes = Vec::new();
for value in 2..=limit {
if smallest[value] == 0 {
smallest[value] = value;
primes.push(value);
}
for &prime in &primes {
if prime > smallest[value] || prime > limit / value {
break;
}
smallest[value * prime] = prime;
}
}
Self { smallest, primes }
}
#[must_use]
pub fn get_factor(&self, value: usize) -> Option<usize> {
assert!(value < self.smallest.len(), "value exceeds sieve limit");
(value >= 2).then(|| self.smallest[value])
}
#[must_use]
pub fn is_prime(&self, value: usize) -> bool {
self.get_factor(value) == Some(value)
}
pub fn prime_factorization(&self, mut value: usize) -> std::vec::IntoIter<usize> {
assert!(value < self.smallest.len(), "value exceeds sieve limit");
let mut factors = Vec::new();
while value >= 2 {
let p = self.smallest[value];
factors.push(p);
value /= p;
}
factors.into_iter()
}
#[must_use]
pub fn primes(&self) -> &[usize] {
&self.primes
}
}
#[derive(Clone, Debug, Eq, PartialEq)]
pub enum SieveRangeError {
ReversedRange,
RangeTooLarge,
}
impl fmt::Display for SieveRangeError {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(f, "{self:?}")
}
}
impl Error for SieveRangeError {}
#[derive(Clone, Debug)]
pub struct EratosthenesSieveRange {
start: u64,
end: u64,
is_prime: Vec<bool>,
base_primes: Vec<u64>,
}
impl EratosthenesSieveRange {
pub fn new(start: u64, end: u64) -> Result<Self, SieveRangeError> {
if start > end {
return Err(SieveRangeError::ReversedRange);
}
let width = end
.checked_sub(start)
.and_then(|x| x.checked_add(1))
.and_then(|x| usize::try_from(x).ok())
.ok_or(SieveRangeError::RangeTooLarge)?;
let root = end.isqrt();
let base = EratosthenesSieve::new(
usize::try_from(root).map_err(|_| SieveRangeError::RangeTooLarge)?,
);
let base_primes = base.primes().iter().map(|&x| x as u64).collect::<Vec<_>>();
let mut is_prime = vec![true; width];
for &p in &base_primes {
let square = p * p;
let first = start.div_ceil(p).saturating_mul(p).max(square);
if first > end {
continue;
}
let mut value = first;
loop {
is_prime[(value - start) as usize] = false;
match value.checked_add(p) {
Some(next) if next <= end => value = next,
_ => break,
}
}
}
if start == 0 {
is_prime[0] = false;
if end >= 1 {
is_prime[1] = false;
}
} else if start == 1 {
is_prime[0] = false;
}
Ok(Self {
start,
end,
is_prime,
base_primes,
})
}
#[must_use]
pub fn is_prime(&self, value: u64) -> Option<bool> {
(self.start..=self.end)
.contains(&value)
.then(|| self.is_prime[(value - self.start) as usize])
}
pub fn prime_factorization(&self, value: u64) -> Option<Vec<u64>> {
if !(self.start..=self.end).contains(&value) {
return None;
}
let mut value = value;
let mut factors = Vec::new();
for &p in &self.base_primes {
if p > value / p {
break;
}
while value.is_multiple_of(p) {
factors.push(p);
value /= p;
}
}
if value > 1 {
factors.push(value);
}
Some(factors)
}
}
}
pub mod fenwicktree {
use std::ops::{Bound, RangeBounds};
#[derive(Clone)]
pub struct FenwickTree<T> {
n: usize,
ary: Vec<T>,
e: T,
}
impl<T: Clone + std::ops::AddAssign<T>> FenwickTree<T> {
pub fn new(n: usize, e: T) -> Self {
FenwickTree {
n,
ary: vec![e.clone(); n],
e,
}
}
pub fn accum(&self, mut idx: usize) -> T {
let mut sum = self.e.clone();
while idx > 0 {
sum += self.ary[idx - 1].clone();
idx &= idx - 1;
}
sum
}
pub fn add<U: Clone>(&mut self, mut idx: usize, val: U)
where
T: std::ops::AddAssign<U>,
{
let n = self.n;
idx += 1;
while idx <= n {
self.ary[idx - 1] += val.clone();
idx += idx & idx.wrapping_neg();
}
}
pub fn sum<R>(&self, range: R) -> T
where
T: std::ops::Sub<Output = T>,
R: RangeBounds<usize>,
{
let r = match range.end_bound() {
Bound::Included(r) => r + 1,
Bound::Excluded(r) => *r,
Bound::Unbounded => self.n,
};
let l = match range.start_bound() {
Bound::Included(l) => *l,
Bound::Excluded(l) => l + 1,
Bound::Unbounded => return self.accum(r),
};
self.accum(r) - self.accum(l)
}
}
}
pub mod fft {
#![warn(missing_docs)]
use crate::Complex;
use std::error::Error;
use std::f64::consts::TAU;
use std::fmt;
#[derive(Copy, Clone, Debug, Eq, PartialEq)]
pub enum FftDirection {
Forward,
Inverse,
}
#[derive(Copy, Clone, Debug, Eq, PartialEq)]
pub enum FftError {
Empty,
NonPowerOfTwo(usize),
}
impl fmt::Display for FftError {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
match self {
Self::Empty => write!(f, "FFT input must not be empty"),
Self::NonPowerOfTwo(n) => write!(f, "FFT length {n} is not a power of two"),
}
}
}
impl Error for FftError {}
pub fn fft(values: &mut [Complex<f64>], direction: FftDirection) -> Result<(), FftError> {
let n = values.len();
if n == 0 {
return Err(FftError::Empty);
}
if !n.is_power_of_two() {
return Err(FftError::NonPowerOfTwo(n));
}
if n == 1 {
return Ok(());
}
let shift = usize::BITS - n.trailing_zeros();
for i in 0..n {
let j = i.reverse_bits() >> shift;
if i < j {
values.swap(i, j);
}
}
let sign = if direction == FftDirection::Forward {
-1.0
} else {
1.0
};
let mut width = 2;
while width <= n {
let (sin, cos) = (sign * TAU / width as f64).sin_cos();
let root = Complex::new(cos, sin);
for chunk in values.chunks_exact_mut(width) {
let mut factor = Complex::new(1.0, 0.0);
let half = width / 2;
for i in 0..half {
let even = chunk[i];
let odd = chunk[i + half] * factor;
chunk[i] = even + odd;
chunk[i + half] = even - odd;
factor *= root;
}
}
width *= 2;
}
if direction == FftDirection::Inverse {
for value in values {
value.re /= n as f64;
value.im /= n as f64;
}
}
Ok(())
}
pub fn fft_padded(
values: &[Complex<f64>],
direction: FftDirection,
) -> Result<Vec<Complex<f64>>, FftError> {
if values.is_empty() {
return Ok(Vec::new());
}
let mut output = values.to_vec();
output.resize(values.len().next_power_of_two(), Complex::new(0.0, 0.0));
fft(&mut output, direction)?;
Ok(output)
}
}
pub mod geometry {
#![warn(missing_docs)]
use std::cmp::Ordering;
use std::ops::{Add, AddAssign, Mul, Neg, Sub, SubAssign};
#[derive(Debug, Copy, Clone, PartialEq, Eq, Default)]
pub struct Point<T> {
pub x: T,
pub y: T,
}
impl<T> Point<T> {
#[must_use]
pub const fn new(x: T, y: T) -> Self {
Self { x, y }
}
}
impl<T: Add<Output = T>> Add for Point<T> {
type Output = Self;
fn add(self, rhs: Self) -> Self {
Self::new(self.x + rhs.x, self.y + rhs.y)
}
}
impl<T: AddAssign> AddAssign for Point<T> {
fn add_assign(&mut self, rhs: Self) {
self.x += rhs.x;
self.y += rhs.y;
}
}
impl<T: Sub<Output = T>> Sub for Point<T> {
type Output = Self;
fn sub(self, rhs: Self) -> Self {
Self::new(self.x - rhs.x, self.y - rhs.y)
}
}
impl<T: SubAssign> SubAssign for Point<T> {
fn sub_assign(&mut self, rhs: Self) {
self.x -= rhs.x;
self.y -= rhs.y;
}
}
#[must_use]
pub fn cross<T>(a: Point<T>, b: Point<T>) -> T
where
T: Mul<Output = T> + Sub<Output = T>,
{
a.x * b.y - a.y * b.x
}
#[must_use]
pub fn dot<T>(a: Point<T>, b: Point<T>) -> T
where
T: Mul<Output = T> + Add<Output = T>,
{
a.x * b.x + a.y * b.y
}
#[must_use]
pub fn norm<T>(a: Point<T>) -> T
where
T: Copy + Mul<Output = T> + Add<Output = T>,
{
dot(a, a)
}
#[derive(Copy, Clone, Debug, Eq, PartialEq)]
pub enum Ccw {
CounterClockwise,
Clockwise,
OnlineBack,
OnlineFront,
OnSegment,
}
#[must_use]
pub fn ccw<T>(a: Point<T>, b: Point<T>, c: Point<T>) -> Ccw
where
T: Copy + Default + PartialOrd + Sub<Output = T> + Mul<Output = T> + Add<Output = T>,
{
let b = b - a;
let c = c - a;
let zero = T::default();
let area = cross(b, c);
if area > zero {
Ccw::CounterClockwise
} else if area < zero {
Ccw::Clockwise
} else if dot(b, c) < zero {
Ccw::OnlineBack
} else if norm(b) < norm(c) {
Ccw::OnlineFront
} else {
Ccw::OnSegment
}
}
fn upper_half<T>(p: &Point<T>) -> bool
where
T: Default + PartialOrd,
{
p.y > T::default() || p.y == T::default() && p.x >= T::default()
}
pub fn polar_cmp<T>(a: &Point<T>, b: &Point<T>) -> Ordering
where
T: Copy + Default + Ord + Mul<Output = T> + Sub<Output = T> + Add<Output = T>,
{
match upper_half(a).cmp(&upper_half(b)).reverse() {
Ordering::Equal => match cross(*a, *b).cmp(&T::default()).reverse() {
Ordering::Equal => norm(*a).cmp(&norm(*b)),
order => order,
},
order => order,
}
}
#[derive(Copy, Clone, Debug, PartialEq, Eq)]
pub struct Line<T> {
a: T,
b: T,
c: T,
}
impl<T> Line<T> {
pub fn coefficients(&self) -> (&T, &T, &T) {
(&self.a, &self.b, &self.c)
}
}
impl<T> Line<T>
where
T: Copy + Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Neg<Output = T>,
{
pub fn from_points(p: Point<T>, q: Point<T>) -> Option<Self>
where
T: PartialEq,
{
if p == q {
return None;
}
let a = p.y - q.y;
let b = q.x - p.x;
let c = -(a * p.x + b * p.y);
Some(Self { a, b, c })
}
pub fn is_parallel(&self, other: &Self) -> bool
where
T: PartialEq,
{
self.a * other.b == other.a * self.b
}
pub fn is_orthogonal(&self, other: &Self) -> bool
where
T: Add<Output = T> + Default + PartialEq,
{
self.a * other.a + self.b * other.b == T::default()
}
}
#[derive(Copy, Clone, Debug, PartialEq)]
pub enum LineIntersection {
Parallel,
Coincident,
Point(Point<f64>),
}
impl Line<f64> {
pub fn intersection(&self, other: &Self, epsilon: f64) -> LineIntersection {
assert!(epsilon >= 0.0);
let determinant = self.a * other.b - other.a * self.b;
if determinant.abs() <= epsilon {
let ac = self.a * other.c - other.a * self.c;
let bc = self.b * other.c - other.b * self.c;
return if ac.abs() <= epsilon && bc.abs() <= epsilon {
LineIntersection::Coincident
} else {
LineIntersection::Parallel
};
}
LineIntersection::Point(Point::new(
(self.b * other.c - other.b * self.c) / determinant,
(self.c * other.a - other.c * self.a) / determinant,
))
}
}
pub fn convex_hull<T>(points: &[Point<T>]) -> Vec<Point<T>>
where
T: Copy + Ord + Sub<Output = T> + Mul<Output = T> + Default,
{
let mut points = points.to_vec();
points.sort_unstable_by(|a, b| a.x.cmp(&b.x).then_with(|| a.y.cmp(&b.y)));
points.dedup();
if points.len() <= 1 {
return points;
}
let turn = |a: Point<T>, b: Point<T>, c: Point<T>| cross(b - a, c - a);
let mut lower = Vec::new();
for &p in &points {
while lower.len() >= 2
&& turn(lower[lower.len() - 2], lower[lower.len() - 1], p) <= T::default()
{
lower.pop();
}
lower.push(p);
}
let mut upper = Vec::new();
for &p in points.iter().rev() {
while upper.len() >= 2
&& turn(upper[upper.len() - 2], upper[upper.len() - 1], p) <= T::default()
{
upper.pop();
}
upper.push(p);
}
lower.pop();
upper.pop();
lower.extend(upper);
lower
}
}
pub mod global {
#![warn(missing_docs)]
#[macro_export]
macro_rules! global_init {
($name:ident : $ty:ty = $initial:expr $(,)?) => {
std::thread_local! { pub static $name: std::cell::RefCell<$ty> = std::cell::RefCell::new($initial); }
};
}
#[macro_export]
macro_rules! global_with {
($name:ident, |$value:ident| $body:expr $(,)?) => {
$name.with_borrow(|$value| $body)
};
}
#[macro_export]
macro_rules! global_with_mut {
($name:ident, |$value:ident| $body:expr $(,)?) => {
$name.with_borrow_mut(|$value| $body)
};
}
}
pub mod graph {
#![warn(missing_docs)]
use std::collections::{BTreeSet, VecDeque};
use std::error::Error;
use std::fmt;
#[derive(Clone, Debug, Eq, PartialEq)]
pub enum GraphError {
InvalidVertex(usize),
NegativeEdge { from: usize, to: usize, weight: i64 },
DistanceOverflow,
NonSquareMatrix,
}
impl fmt::Display for GraphError {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(f, "{self:?}")
}
}
impl Error for GraphError {}
#[derive(Clone, Debug, Eq, PartialEq)]
pub struct BellmanFordResult {
pub distances: Vec<Option<i64>>,
pub affected_by_negative_cycle: Vec<bool>,
}
fn validate(graph: &[Vec<(usize, i64)>], start: usize) -> Result<(), GraphError> {
if start >= graph.len() {
return Err(GraphError::InvalidVertex(start));
}
for edges in graph {
for &(to, _) in edges {
if to >= graph.len() {
return Err(GraphError::InvalidVertex(to));
}
}
}
Ok(())
}
pub fn bellman_ford(
graph: &[Vec<(usize, i64)>],
start: usize,
) -> Result<BellmanFordResult, GraphError> {
validate(graph, start)?;
let n = graph.len();
let sources = graph
.iter()
.enumerate()
.filter(|(_, edges)| !edges.is_empty())
.collect::<Vec<_>>();
let mut distances: Vec<Option<i64>> = vec![None; n];
distances[start] = Some(0);
for _ in 1..n {
let mut changed = false;
for &(from, edges) in &sources {
let Some(base) = distances[from] else {
continue;
};
for &(to, weight) in edges {
let candidate = base
.checked_add(weight)
.ok_or(GraphError::DistanceOverflow)?;
if distances[to].is_none_or(|old| candidate < old) {
distances[to] = Some(candidate);
changed = true;
}
}
}
if !changed {
break;
}
}
let mut affected = vec![false; n];
let mut queue = VecDeque::new();
for &(from, edges) in &sources {
let Some(base) = distances[from] else {
continue;
};
for &(to, weight) in edges {
let candidate = base
.checked_add(weight)
.ok_or(GraphError::DistanceOverflow)?;
if distances[to].is_some_and(|old| candidate < old) && !affected[to] {
affected[to] = true;
queue.push_back(to);
}
}
}
while let Some(from) = queue.pop_front() {
for &(to, _) in &graph[from] {
if !affected[to] {
affected[to] = true;
queue.push_back(to);
}
}
}
Ok(BellmanFordResult {
distances,
affected_by_negative_cycle: affected,
})
}
pub fn dijkstra(
graph: &[Vec<(usize, i64)>],
start: usize,
) -> Result<Vec<Option<i64>>, GraphError> {
validate(graph, start)?;
for (from, edges) in graph.iter().enumerate() {
for &(to, weight) in edges {
if weight < 0 {
return Err(GraphError::NegativeEdge { from, to, weight });
}
}
}
let mut distances: Vec<Option<i64>> = vec![None; graph.len()];
let mut queue = BTreeSet::new();
distances[start] = Some(0);
queue.insert((0_i64, start));
while let Some((distance, from)) = queue.pop_first() {
for &(to, weight) in &graph[from] {
let candidate = distance
.checked_add(weight)
.ok_or(GraphError::DistanceOverflow)?;
if distances[to].is_none_or(|old| candidate < old) {
if let Some(old) = distances[to] {
queue.remove(&(old, to));
}
distances[to] = Some(candidate);
queue.insert((candidate, to));
}
}
}
Ok(distances)
}
pub fn warshall_floyd(matrix: &mut [Vec<Option<i64>>]) -> Result<(), GraphError> {
let n = matrix.len();
if matrix.iter().any(|row| row.len() != n) {
return Err(GraphError::NonSquareMatrix);
}
for k in 0..n {
for i in 0..n {
for j in 0..n {
if let (Some(a), Some(b)) = (matrix[i][k], matrix[k][j]) {
let candidate = a.checked_add(b).ok_or(GraphError::DistanceOverflow)?;
if matrix[i][j].is_none_or(|old| candidate < old) {
matrix[i][j] = Some(candidate);
}
}
}
}
}
Ok(())
}
}
pub mod graph_search {
#![warn(missing_docs)]
use std::collections::VecDeque;
use std::error::Error;
use std::fmt;
#[derive(Debug, PartialEq, Eq, Clone, Copy)]
pub enum DfsOrder {
Preorder,
Postorder,
}
#[derive(Clone, Debug, Eq, PartialEq)]
pub enum GraphSearchError {
InvalidVertex(usize),
CycleDetected,
}
impl fmt::Display for GraphSearchError {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(f, "{self:?}")
}
}
impl Error for GraphSearchError {}
fn validate(graph: &[Vec<usize>]) -> Result<(), GraphSearchError> {
for edges in graph {
for &to in edges {
if to >= graph.len() {
return Err(GraphSearchError::InvalidVertex(to));
}
}
}
Ok(())
}
pub fn topological_sort(graph: &[Vec<usize>]) -> Result<Vec<usize>, GraphSearchError> {
validate(graph)?;
let mut indegree = vec![0; graph.len()];
for edges in graph {
for &to in edges {
indegree[to] += 1;
}
}
let mut queue = (0..graph.len())
.filter(|&v| indegree[v] == 0)
.collect::<VecDeque<_>>();
let mut order = Vec::with_capacity(graph.len());
while let Some(from) = queue.pop_front() {
order.push(from);
for &to in &graph[from] {
indegree[to] -= 1;
if indegree[to] == 0 {
queue.push_back(to);
}
}
}
if order.len() == graph.len() {
Ok(order)
} else {
Err(GraphSearchError::CycleDetected)
}
}
pub fn graph_dfs(
start: usize,
graph: &[Vec<usize>],
order: DfsOrder,
) -> Result<Vec<usize>, GraphSearchError> {
validate(graph)?;
if start >= graph.len() {
return Err(GraphSearchError::InvalidVertex(start));
}
let mut visited = vec![false; graph.len()];
let mut output = Vec::new();
let mut stack = vec![(start, 0)];
visited[start] = true;
if order == DfsOrder::Preorder {
output.push(start);
}
while let Some((vertex, next_edge)) = stack.last_mut() {
if *next_edge == graph[*vertex].len() {
let (vertex, _) = stack.pop().unwrap();
if order == DfsOrder::Postorder {
output.push(vertex);
}
continue;
}
let next = graph[*vertex][*next_edge];
*next_edge += 1;
if !visited[next] {
visited[next] = true;
if order == DfsOrder::Preorder {
output.push(next);
}
stack.push((next, 0));
}
}
Ok(output)
}
pub fn graph_bfs(start: usize, graph: &[Vec<usize>]) -> Result<Vec<usize>, GraphSearchError> {
validate(graph)?;
if start >= graph.len() {
return Err(GraphSearchError::InvalidVertex(start));
}
let mut visited = vec![false; graph.len()];
let mut queue = VecDeque::from([start]);
let mut output = Vec::new();
visited[start] = true;
while let Some(from) = queue.pop_front() {
output.push(from);
for &to in &graph[from] {
if !visited[to] {
visited[to] = true;
queue.push_back(to);
}
}
}
Ok(output)
}
}
pub mod input {
#![warn(missing_docs)]
use std::error::Error;
use std::fmt;
use std::io::{self, Read};
use std::str::FromStr;
#[derive(Debug)]
pub enum ScanError {
Io(io::Error),
UnexpectedEof,
Parse {
index: usize,
token: String,
message: String,
},
}
impl fmt::Display for ScanError {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
match self {
Self::Io(e) => write!(f, "input I/O error: {e}"),
Self::UnexpectedEof => write!(f, "unexpected end of input"),
Self::Parse {
index,
token,
message,
} => write!(f, "cannot parse token {index} ({token:?}): {message}"),
}
}
}
impl Error for ScanError {
fn source(&self) -> Option<&(dyn Error + 'static)> {
match self {
Self::Io(e) => Some(e),
_ => None,
}
}
}
impl From<io::Error> for ScanError {
fn from(value: io::Error) -> Self {
Self::Io(value)
}
}
#[derive(Clone, Debug)]
pub struct Scanner {
tokens: Vec<String>,
index: usize,
}
impl Scanner {
pub fn new(mut reader: impl Read) -> Result<Self, ScanError> {
let mut source = String::new();
reader.read_to_string(&mut source)?;
Ok(Self {
tokens: source.split_whitespace().map(str::to_owned).collect(),
index: 0,
})
}
pub fn stdin() -> Result<Self, ScanError> {
Self::new(io::stdin())
}
pub fn next_token<T>(&mut self) -> Result<T, ScanError>
where
T: FromStr,
T::Err: fmt::Display,
{
let token = self
.tokens
.get(self.index)
.ok_or(ScanError::UnexpectedEof)?
.clone();
let index = self.index;
self.index += 1;
token.parse().map_err(|e: T::Err| ScanError::Parse {
index,
token,
message: e.to_string(),
})
}
pub fn next_vec<T>(&mut self, len: usize) -> Result<Vec<T>, ScanError>
where
T: FromStr,
T::Err: fmt::Display,
{
(0..len).map(|_| self.next_token()).collect()
}
#[must_use]
pub fn remaining(&self) -> usize {
self.tokens.len() - self.index
}
}
#[macro_export]
macro_rules! inputv{($scanner:expr;)=>{};($scanner:expr;mut $name:ident:$ty:ty $(,$($rest:tt)*)?)=>{let mut $name:$ty=$scanner.next_token().expect("failed to parse input");$crate::inputv!($scanner;$($($rest)*)?);};($scanner:expr;$name:ident:$ty:ty $(,$($rest:tt)*)?)=>{let $name:$ty=$scanner.next_token().expect("failed to parse input");$crate::inputv!($scanner;$($($rest)*)?);};}
pub trait ToCharVec {
fn to_charvec(&self) -> Vec<char>;
}
impl ToCharVec for str {
fn to_charvec(&self) -> Vec<char> {
self.chars().collect()
}
}
impl ToCharVec for String {
fn to_charvec(&self) -> Vec<char> {
self.as_str().to_charvec()
}
}
}
pub mod internal_bit {
#[allow(dead_code)]
pub(crate) fn ceil_pow2(n: u32) -> u32 {
32 - n.saturating_sub(1).leading_zeros()
}
}
pub mod internal_math {
#![allow(dead_code)]
use std::{mem::swap, num::Wrapping as W};
/* const */
pub(crate) fn safe_mod(mut x: i64, m: i64) -> i64 {
x %= m;
if x < 0 {
x += m;
}
x
}
pub(crate) struct Barrett {
pub(crate) _m: u32,
pub(crate) im: u64,
}
impl Barrett {
pub(crate) fn new(m: u32) -> Barrett {
Barrett {
_m: m,
im: (-1i64 as u64 / m as u64).wrapping_add(1),
}
}
pub(crate) fn umod(&self) -> u32 {
self._m
}
#[allow(clippy::many_single_char_names)]
pub(crate) fn mul(&self, a: u32, b: u32) -> u32 {
mul_mod(a, b, self._m, self.im)
}
}
#[allow(clippy::many_single_char_names)]
pub(crate) fn mul_mod(a: u32, b: u32, m: u32, im: u64) -> u32 {
let mut z = a as u64;
z *= b as u64;
let x = (((z as u128) * (im as u128)) >> 64) as u64;
let mut v = z.wrapping_sub(x.wrapping_mul(m as u64)) as u32;
if m <= v {
v = v.wrapping_add(m);
}
v
}
/* const */
#[allow(clippy::many_single_char_names)]
pub(crate) fn pow_mod(x: i64, mut n: i64, m: i32) -> i64 {
if m == 1 {
return 0;
}
let _m = m as u32;
let mut r: u64 = 1;
let mut y: u64 = safe_mod(x, m as i64) as u64;
while n != 0 {
if (n & 1) > 0 {
r = (r * y) % (_m as u64);
}
y = (y * y) % (_m as u64);
n >>= 1;
}
r as i64
}
/* const */
pub(crate) fn is_prime(n: i32) -> bool {
let n = n as i64;
match n {
_ if n <= 1 => return false,
2 | 7 | 61 => return true,
_ if n % 2 == 0 => return false,
_ => {}
}
let mut d = n - 1;
while d % 2 == 0 {
d /= 2;
}
for &a in &[2, 7, 61] {
let mut t = d;
let mut y = pow_mod(a, t, n as i32);
while t != n - 1 && y != 1 && y != n - 1 {
y = y * y % n;
t <<= 1;
}
if y != n - 1 && t % 2 == 0 {
return false;
}
}
true
}
/* const */
#[allow(clippy::many_single_char_names)]
pub(crate) fn inv_gcd(a: i64, b: i64) -> (i64, i64) {
let a = safe_mod(a, b);
if a == 0 {
return (b, 0);
}
let mut s = b;
let mut t = a;
let mut m0 = 0;
let mut m1 = 1;
while t != 0 {
let u = s / t;
s -= t * u;
m0 -= m1 * u; // |m1 * u| <= |m1| * s <= b
swap(&mut s, &mut t);
swap(&mut m0, &mut m1);
}
if m0 < 0 {
m0 += b / s;
}
(s, m0)
}
/* const */
pub(crate) fn primitive_root(m: i32) -> i32 {
match m {
2 => return 1,
167_772_161 => return 3,
469_762_049 => return 3,
754_974_721 => return 11,
998_244_353 => return 3,
_ => {}
}
let mut divs = [0; 20];
divs[0] = 2;
let mut cnt = 1;
let mut x = (m - 1) / 2;
while x % 2 == 0 {
x /= 2;
}
for i in (3..i32::MAX).step_by(2) {
if i as i64 * i as i64 > x as i64 {
break;
}
if x % i == 0 {
divs[cnt] = i;
cnt += 1;
while x % i == 0 {
x /= i;
}
}
}
if x > 1 {
divs[cnt] = x;
cnt += 1;
}
let mut g = 2;
loop {
if (0..cnt).all(|i| pow_mod(g, ((m - 1) / divs[i]) as i64, m) != 1) {
break g as i32;
}
g += 1;
}
}
/* const */
#[allow(clippy::many_single_char_names)]
pub(crate) fn floor_sum_unsigned(
mut n: W<u64>,
mut m: W<u64>,
mut a: W<u64>,
mut b: W<u64>,
) -> W<u64> {
let mut ans = W(0);
loop {
if a >= m {
if n > W(0) {
ans += n * (n - W(1)) / W(2) * (a / m);
}
a %= m;
}
if b >= m {
ans += n * (b / m);
b %= m;
}
let y_max = a * n + b;
if y_max < m {
break;
}
n = y_max / m;
b = y_max % m;
std::mem::swap(&mut m, &mut a);
}
ans
}
}
pub mod internal_queue {
#![allow(dead_code)]
#[derive(Default)]
pub(crate) struct SimpleQueue<T> {
payload: Vec<T>,
pos: usize,
}
impl<T> SimpleQueue<T> {
pub(crate) fn reserve(&mut self, n: usize) {
if n > self.payload.len() {
self.payload.reserve(n - self.payload.len());
}
}
pub(crate) fn size(&self) -> usize {
self.payload.len() - self.pos
}
pub(crate) fn empty(&self) -> bool {
self.pos == self.payload.len()
}
pub(crate) fn push(&mut self, t: T) {
self.payload.push(t);
}
pub(crate) fn front(&self) -> Option<&T> {
if self.pos < self.payload.len() {
Some(&self.payload[self.pos])
} else {
None
}
}
pub(crate) fn clear(&mut self) {
self.payload.clear();
self.pos = 0;
}
pub(crate) fn pop(&mut self) -> Option<&T> {
if self.pos < self.payload.len() {
self.pos += 1;
Some(&self.payload[self.pos - 1])
} else {
None
}
}
}
}
pub mod internal_scc {
pub struct Csr<E> {
start: Vec<usize>,
elist: Vec<E>,
}
impl<E> Csr<E>
where
E: Copy,
{
pub fn new(n: usize, edges: &[(usize, E)], init: E) -> Self {
let mut csr = Csr {
start: vec![0; n + 1],
elist: vec![init; edges.len()],
};
for e in edges.iter() {
csr.start[e.0 + 1] += 1;
}
for i in 1..=n {
csr.start[i] += csr.start[i - 1];
}
let mut counter = csr.start.clone();
for e in edges.iter() {
csr.elist[counter[e.0]] = e.1;
counter[e.0] += 1;
}
csr
}
}
#[derive(Copy, Clone, Debug)]
struct _Edge {
to: usize,
}
#[derive(Clone, Debug)]
pub struct SccGraph {
n: usize,
edges: Vec<(usize, _Edge)>,
}
impl SccGraph {
pub fn new(n: usize) -> Self {
SccGraph { n, edges: vec![] }
}
pub fn num_vertices(&self) -> usize {
self.n
}
pub fn add_edge(&mut self, from: usize, to: usize) {
self.edges.push((from, _Edge { to }));
}
pub fn scc_ids(&self) -> (usize, Vec<usize>) {
struct _Env {
g: Csr<_Edge>,
now_ord: usize,
group_num: usize,
visited: Vec<usize>,
low: Vec<usize>,
ord: Vec<Option<usize>>,
ids: Vec<usize>,
}
let mut env = _Env {
g: Csr::new(self.n, &self.edges, _Edge { to: 0 }),
now_ord: 0,
group_num: 0,
visited: Vec::with_capacity(self.n),
low: vec![0; self.n],
ord: vec![None; self.n],
ids: vec![0; self.n],
};
fn dfs(v: usize, n: usize, env: &mut _Env) {
env.low[v] = env.now_ord;
env.ord[v] = Some(env.now_ord);
env.now_ord += 1;
env.visited.push(v);
for i in env.g.start[v]..env.g.start[v + 1] {
let to = env.g.elist[i].to;
if let Some(x) = env.ord[to] {
env.low[v] = std::cmp::min(env.low[v], x);
} else {
dfs(to, n, env);
env.low[v] = std::cmp::min(env.low[v], env.low[to]);
}
}
if env.low[v] == env.ord[v].unwrap() {
loop {
let u = *env.visited.last().unwrap();
env.visited.pop();
env.ord[u] = Some(n);
env.ids[u] = env.group_num;
if u == v {
break;
}
}
env.group_num += 1;
}
}
for i in 0..self.n {
if env.ord[i].is_none() {
dfs(i, self.n, &mut env);
}
}
for x in env.ids.iter_mut() {
*x = env.group_num - 1 - *x;
}
(env.group_num, env.ids)
}
#[allow(clippy::redundant_closure)]
pub fn scc(&self) -> Vec<Vec<usize>> {
let ids = self.scc_ids();
let group_num = ids.0;
let mut counts = vec![0usize; group_num];
for &x in ids.1.iter() {
counts[x] += 1;
}
let mut groups: Vec<Vec<usize>> = (0..ids.0).map(|_| vec![]).collect();
for i in 0..group_num {
groups[i].reserve(counts[i]);
}
for i in 0..self.n {
groups[ids.1[i]].push(i);
}
groups
}
}
}
pub mod internal_type_traits {
use std::{
fmt,
iter::{Product, Sum},
ops::{
Add, AddAssign, BitAnd, BitAndAssign, BitOr, BitOrAssign, BitXor, BitXorAssign, Div,
DivAssign, Mul, MulAssign, Not, Rem, RemAssign, Shl, ShlAssign, Shr, ShrAssign, Sub,
SubAssign,
},
};
pub trait Integral:
'static
+ Send
+ Sync
+ Copy
+ Ord
+ Not<Output = Self>
+ Add<Output = Self>
+ Sub<Output = Self>
+ Mul<Output = Self>
+ Div<Output = Self>
+ Rem<Output = Self>
+ AddAssign
+ SubAssign
+ MulAssign
+ DivAssign
+ RemAssign
+ Sum
+ Product
+ BitOr<Output = Self>
+ BitAnd<Output = Self>
+ BitXor<Output = Self>
+ BitOrAssign
+ BitAndAssign
+ BitXorAssign
+ Shl<Output = Self>
+ Shr<Output = Self>
+ ShlAssign
+ ShrAssign
+ fmt::Display
+ fmt::Debug
+ fmt::Binary
+ fmt::Octal
+ Zero
+ One
+ BoundedBelow
+ BoundedAbove
{
}
pub trait Zero {
fn zero() -> Self;
}
pub trait One {
fn one() -> Self;
}
pub trait BoundedBelow {
fn min_value() -> Self;
}
pub trait BoundedAbove {
fn max_value() -> Self;
}
macro_rules! impl_integral {
($($ty:ty),*) => {
$(
impl Zero for $ty {
#[inline]
fn zero() -> Self {
0
}
}
impl One for $ty {
#[inline]
fn one() -> Self {
1
}
}
impl BoundedBelow for $ty {
#[inline]
fn min_value() -> Self {
Self::MIN
}
}
impl BoundedAbove for $ty {
#[inline]
fn max_value() -> Self {
Self::MAX
}
}
impl Integral for $ty {}
)*
};
}
impl_integral!(i8, i16, i32, i64, i128, isize, u8, u16, u32, u64, u128, usize);
}
pub mod lazysegtree {
use crate::internal_bit::ceil_pow2;
use crate::Monoid;
pub trait MapMonoid {
type M: Monoid;
type F: Clone;
fn identity_element() -> <Self::M as Monoid>::S {
Self::M::identity()
}
fn binary_operation(
a: &<Self::M as Monoid>::S,
b: &<Self::M as Monoid>::S,
) -> <Self::M as Monoid>::S {
Self::M::binary_operation(a, b)
}
fn identity_map() -> Self::F;
fn mapping(f: &Self::F, x: &<Self::M as Monoid>::S) -> <Self::M as Monoid>::S;
fn composition(f: &Self::F, g: &Self::F) -> Self::F;
}
impl<F: MapMonoid> Default for LazySegtree<F> {
fn default() -> Self {
Self::new(0)
}
}
impl<F: MapMonoid> LazySegtree<F> {
pub fn new(n: usize) -> Self {
vec![F::identity_element(); n].into()
}
}
impl<F: MapMonoid> From<Vec<<F::M as Monoid>::S>> for LazySegtree<F> {
fn from(v: Vec<<F::M as Monoid>::S>) -> Self {
let n = v.len();
let log = ceil_pow2(n as u32) as usize;
let size = 1 << log;
let mut d = vec![F::identity_element(); 2 * size];
let lz = vec![F::identity_map(); size];
d[size..(size + n)].clone_from_slice(&v);
let mut ret = LazySegtree {
n,
size,
log,
d,
lz,
};
for i in (1..size).rev() {
ret.update(i);
}
ret
}
}
impl<F: MapMonoid> LazySegtree<F> {
pub fn set(&mut self, mut p: usize, x: <F::M as Monoid>::S) {
assert!(p < self.n);
p += self.size;
for i in (1..=self.log).rev() {
self.push(p >> i);
}
self.d[p] = x;
for i in 1..=self.log {
self.update(p >> i);
}
}
pub fn get(&mut self, mut p: usize) -> <F::M as Monoid>::S {
assert!(p < self.n);
p += self.size;
for i in (1..=self.log).rev() {
self.push(p >> i);
}
self.d[p].clone()
}
pub fn prod<R>(&mut self, range: R) -> <F::M as Monoid>::S
where
R: RangeBounds<usize>,
{
if range.start_bound() == Bound::Unbounded && range.end_bound() == Bound::Unbounded {
return self.all_prod();
}
let mut r = match range.end_bound() {
Bound::Included(r) => r + 1,
Bound::Excluded(r) => *r,
Bound::Unbounded => self.n,
};
let mut l = match range.start_bound() {
Bound::Included(l) => *l,
Bound::Excluded(l) => l + 1,
Bound::Unbounded => 0,
};
assert!(l <= r && r <= self.n);
if l == r {
return F::identity_element();
}
l += self.size;
r += self.size;
for i in (1..=self.log).rev() {
if ((l >> i) << i) != l {
self.push(l >> i);
}
if ((r >> i) << i) != r {
self.push(r >> i);
}
}
let mut sml = F::identity_element();
let mut smr = F::identity_element();
while l < r {
if l & 1 != 0 {
sml = F::binary_operation(&sml, &self.d[l]);
l += 1;
}
if r & 1 != 0 {
r -= 1;
smr = F::binary_operation(&self.d[r], &smr);
}
l >>= 1;
r >>= 1;
}
F::binary_operation(&sml, &smr)
}
pub fn all_prod(&self) -> <F::M as Monoid>::S {
self.d[1].clone()
}
pub fn apply(&mut self, mut p: usize, f: F::F) {
assert!(p < self.n);
p += self.size;
for i in (1..=self.log).rev() {
self.push(p >> i);
}
self.d[p] = F::mapping(&f, &self.d[p]);
for i in 1..=self.log {
self.update(p >> i);
}
}
pub fn apply_range<R>(&mut self, range: R, f: F::F)
where
R: RangeBounds<usize>,
{
let mut r = match range.end_bound() {
Bound::Included(r) => r + 1,
Bound::Excluded(r) => *r,
Bound::Unbounded => self.n,
};
let mut l = match range.start_bound() {
Bound::Included(l) => *l,
Bound::Excluded(l) => l + 1,
Bound::Unbounded => 0,
};
assert!(l <= r && r <= self.n);
if l == r {
return;
}
l += self.size;
r += self.size;
for i in (1..=self.log).rev() {
if ((l >> i) << i) != l {
self.push(l >> i);
}
if ((r >> i) << i) != r {
self.push((r - 1) >> i);
}
}
{
let l2 = l;
let r2 = r;
while l < r {
if l & 1 != 0 {
self.all_apply(l, f.clone());
l += 1;
}
if r & 1 != 0 {
r -= 1;
self.all_apply(r, f.clone());
}
l >>= 1;
r >>= 1;
}
l = l2;
r = r2;
}
for i in 1..=self.log {
if ((l >> i) << i) != l {
self.update(l >> i);
}
if ((r >> i) << i) != r {
self.update((r - 1) >> i);
}
}
}
pub fn max_right<G>(&mut self, mut l: usize, g: G) -> usize
where
G: Fn(<F::M as Monoid>::S) -> bool,
{
assert!(l <= self.n);
assert!(g(F::identity_element()));
if l == self.n {
return self.n;
}
l += self.size;
for i in (1..=self.log).rev() {
self.push(l >> i);
}
let mut sm = F::identity_element();
while {
while l.is_multiple_of(2) {
l >>= 1;
}
if !g(F::binary_operation(&sm, &self.d[l])) {
while l < self.size {
self.push(l);
l *= 2;
let res = F::binary_operation(&sm, &self.d[l]);
if g(res.clone()) {
sm = res;
l += 1;
}
}
return l - self.size;
}
sm = F::binary_operation(&sm, &self.d[l]);
l += 1;
{
let l = l as isize;
(l & -l) != l
}
} {}
self.n
}
pub fn min_left<G>(&mut self, mut r: usize, g: G) -> usize
where
G: Fn(<F::M as Monoid>::S) -> bool,
{
assert!(r <= self.n);
assert!(g(F::identity_element()));
if r == 0 {
return 0;
}
r += self.size;
for i in (1..=self.log).rev() {
self.push((r - 1) >> i);
}
let mut sm = F::identity_element();
while {
r -= 1;
while r > 1 && !r.is_multiple_of(2) {
r >>= 1;
}
if !g(F::binary_operation(&self.d[r], &sm)) {
while r < self.size {
self.push(r);
r = 2 * r + 1;
let res = F::binary_operation(&self.d[r], &sm);
if g(res.clone()) {
sm = res;
r -= 1;
}
}
return r + 1 - self.size;
}
sm = F::binary_operation(&self.d[r], &sm);
{
let r = r as isize;
(r & -r) != r
}
} {}
0
}
}
#[derive(Clone)]
pub struct LazySegtree<F>
where
F: MapMonoid,
{
n: usize,
size: usize,
log: usize,
d: Vec<<F::M as Monoid>::S>,
lz: Vec<F::F>,
}
impl<F> LazySegtree<F>
where
F: MapMonoid,
{
fn update(&mut self, k: usize) {
self.d[k] = F::binary_operation(&self.d[2 * k], &self.d[2 * k + 1]);
}
fn all_apply(&mut self, k: usize, f: F::F) {
self.d[k] = F::mapping(&f, &self.d[k]);
if k < self.size {
self.lz[k] = F::composition(&f, &self.lz[k]);
}
}
fn push(&mut self, k: usize) {
self.all_apply(2 * k, self.lz[k].clone());
self.all_apply(2 * k + 1, self.lz[k].clone());
self.lz[k] = F::identity_map();
}
}
use std::{
fmt::{Debug, Error, Formatter, Write},
ops::{Bound, RangeBounds},
};
impl<F> Debug for LazySegtree<F>
where
F: MapMonoid,
F::F: Debug,
<F::M as Monoid>::S: Debug,
{
fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), Error> {
for i in 0..self.log {
for j in 0..1 << i {
f.write_fmt(format_args!(
"{:?}[{:?}]\t",
self.d[(1 << i) + j],
self.lz[(1 << i) + j]
))?;
}
f.write_char('\n')?;
}
for i in 0..self.size {
f.write_fmt(format_args!("{:?}\t", self.d[self.size + i]))?;
}
Ok(())
}
}
}
pub mod lca {
#![warn(missing_docs)]
use std::collections::VecDeque;
use std::error::Error;
use std::fmt;
#[derive(Clone, Debug, Eq, PartialEq)]
pub enum LcaError {
InvalidVertex(usize),
NotATree,
}
impl fmt::Display for LcaError {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(f, "{self:?}")
}
}
impl Error for LcaError {}
#[derive(Clone, Debug)]
pub struct Lca {
parent: Vec<Vec<usize>>,
depth: Vec<usize>,
}
impl Lca {
pub fn new(
vertex_count: usize,
edges: &[(usize, usize)],
root: usize,
) -> Result<Self, LcaError> {
if root >= vertex_count {
return Err(LcaError::InvalidVertex(root));
}
if edges.len() != vertex_count.saturating_sub(1) {
return Err(LcaError::NotATree);
}
let mut graph = vec![Vec::new(); vertex_count];
for &(a, b) in edges {
if a >= vertex_count {
return Err(LcaError::InvalidVertex(a));
}
if b >= vertex_count {
return Err(LcaError::InvalidVertex(b));
}
graph[a].push(b);
graph[b].push(a);
}
let levels = (usize::BITS - vertex_count.max(1).leading_zeros()) as usize;
let mut parent = vec![vec![root; vertex_count]; levels];
let mut depth = vec![0; vertex_count];
let mut seen = vec![false; vertex_count];
let mut queue = VecDeque::from([root]);
seen[root] = true;
while let Some(from) = queue.pop_front() {
for &to in &graph[from] {
if seen[to] {
if parent[0][from] != to {
return Err(LcaError::NotATree);
}
continue;
}
seen[to] = true;
parent[0][to] = from;
depth[to] = depth[from] + 1;
queue.push_back(to);
}
}
if seen.iter().any(|&x| !x) {
return Err(LcaError::NotATree);
}
for level in 1..levels {
for v in 0..vertex_count {
parent[level][v] = parent[level - 1][parent[level - 1][v]];
}
}
Ok(Self { parent, depth })
}
pub fn lca(&self, mut a: usize, mut b: usize) -> Result<usize, LcaError> {
self.check(a)?;
self.check(b)?;
if self.depth[a] > self.depth[b] {
std::mem::swap(&mut a, &mut b);
}
b = self.kth_ancestor(b, self.depth[b] - self.depth[a]).unwrap();
if a == b {
return Ok(a);
}
for level in (0..self.parent.len()).rev() {
if self.parent[level][a] != self.parent[level][b] {
a = self.parent[level][a];
b = self.parent[level][b];
}
}
Ok(self.parent[0][a])
}
pub fn depth(&self, vertex: usize) -> Result<usize, LcaError> {
self.check(vertex)?;
Ok(self.depth[vertex])
}
pub fn distance(&self, a: usize, b: usize) -> Result<usize, LcaError> {
let p = self.lca(a, b)?;
Ok(self.depth[a] + self.depth[b] - 2 * self.depth[p])
}
pub fn kth_ancestor(&self, mut vertex: usize, steps: usize) -> Option<usize> {
if vertex >= self.depth.len() || steps > self.depth[vertex] {
return None;
}
for level in 0..self.parent.len() {
if steps >> level & 1 != 0 {
vertex = self.parent[level][vertex];
}
}
Some(vertex)
}
fn check(&self, v: usize) -> Result<(), LcaError> {
if v < self.depth.len() {
Ok(())
} else {
Err(LcaError::InvalidVertex(v))
}
}
}
}
pub mod math {
use crate::internal_math;
use std::mem::swap;
#[allow(clippy::many_single_char_names)]
pub fn pow_mod(x: i64, mut n: i64, m: u32) -> u32 {
assert!(0 <= n && 1 <= m && m <= 2u32.pow(31));
if m == 1 {
return 0;
}
let bt = internal_math::Barrett::new(m);
let mut r = 1;
let mut y = internal_math::safe_mod(x, m as i64) as u32;
while n != 0 {
if n & 1 != 0 {
r = bt.mul(r, y);
}
y = bt.mul(y, y);
n >>= 1;
}
r
}
pub fn inv_mod(x: i64, m: i64) -> i64 {
assert!(1 <= m);
let z = internal_math::inv_gcd(x, m);
assert!(z.0 == 1);
z.1
}
pub fn crt(r: &[i64], m: &[i64]) -> (i64, i64) {
assert_eq!(r.len(), m.len());
let (mut r0, mut m0) = (0, 1);
for (&(mut ri), &(mut mi)) in r.iter().zip(m.iter()) {
assert!(1 <= mi);
ri = internal_math::safe_mod(ri, mi);
if m0 < mi {
swap(&mut r0, &mut ri);
swap(&mut m0, &mut mi);
}
if m0 % mi == 0 {
if r0 % mi != ri {
return (0, 0);
}
continue;
}
let (g, im) = internal_math::inv_gcd(m0, mi);
let u1 = mi / g;
if (ri - r0) % g != 0 {
return (0, 0);
}
let x = (ri - r0) / g % u1 * im % u1;
r0 += x * m0;
m0 *= u1; // -> lcm(m0, mi)
if r0 < 0 {
r0 += m0
};
}
(r0, m0)
}
#[allow(clippy::many_single_char_names)]
pub fn floor_sum(n: i64, m: i64, a: i64, b: i64) -> i64 {
use std::num::Wrapping as W;
assert!((0..1i64 << 32).contains(&n));
assert!((1..1i64 << 32).contains(&m));
let mut ans = W(0_u64);
let (wn, wm, mut wa, mut wb) = (W(n as u64), W(m as u64), W(a as u64), W(b as u64));
if a < 0 {
let a2 = W(internal_math::safe_mod(a, m) as u64);
ans -= wn * (wn - W(1)) / W(2) * ((a2 - wa) / wm);
wa = a2;
}
if b < 0 {
let b2 = W(internal_math::safe_mod(b, m) as u64);
ans -= wn * ((b2 - wb) / wm);
wb = b2;
}
let ret = ans + internal_math::floor_sum_unsigned(wn, wm, wa, wb);
ret.0 as i64
}
}
pub mod matrix {
#![warn(missing_docs)]
use std::error::Error;
use std::fmt;
use std::ops::{Add, AddAssign, Index, IndexMut, Mul};
#[derive(Clone, Debug, Eq, PartialEq)]
pub enum MatrixError {
RaggedRows,
DimensionMismatch,
NotSquare,
}
impl fmt::Display for MatrixError {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(f, "{self:?}")
}
}
impl Error for MatrixError {}
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Matrix<T> {
rows: usize,
cols: usize,
data: Vec<T>,
}
impl<T: Clone> Matrix<T> {
#[must_use]
pub fn new(rows: usize, cols: usize, value: T) -> Self {
Self {
rows,
cols,
data: vec![value; rows.checked_mul(cols).expect("matrix size overflow")],
}
}
pub fn from_rows(rows: Vec<Vec<T>>) -> Result<Self, MatrixError> {
let row_count = rows.len();
let cols = rows.first().map_or(0, Vec::len);
if rows.iter().any(|r| r.len() != cols) {
return Err(MatrixError::RaggedRows);
}
Ok(Self {
rows: row_count,
cols,
data: rows.into_iter().flatten().collect(),
})
}
pub fn transpose(&self) -> Self {
if self.is_empty() {
return Self {
rows: self.cols,
cols: self.rows,
data: Vec::new(),
};
}
let mut data = Vec::with_capacity(self.data.len());
for col in 0..self.cols {
for row in 0..self.rows {
data.push(self[(row, col)].clone());
}
}
Self {
rows: self.cols,
cols: self.rows,
data,
}
}
}
impl<T> Matrix<T> {
#[must_use]
pub fn rows(&self) -> usize {
self.rows
}
#[must_use]
pub fn cols(&self) -> usize {
self.cols
}
#[must_use]
pub fn is_empty(&self) -> bool {
self.rows == 0 || self.cols == 0
}
pub fn swap_rows(&mut self, a: usize, b: usize) {
assert!(a < self.rows && b < self.rows);
for c in 0..self.cols {
self.data.swap(a * self.cols + c, b * self.cols + c);
}
}
pub fn swap_columns(&mut self, a: usize, b: usize) {
assert!(a < self.cols && b < self.cols);
for r in 0..self.rows {
self.data.swap(r * self.cols + a, r * self.cols + b);
}
}
}
impl<T> Index<(usize, usize)> for Matrix<T> {
type Output = T;
fn index(&self, (r, c): (usize, usize)) -> &T {
assert!(r < self.rows && c < self.cols);
&self.data[r * self.cols + c]
}
}
impl<T> IndexMut<(usize, usize)> for Matrix<T> {
fn index_mut(&mut self, (r, c): (usize, usize)) -> &mut T {
assert!(r < self.rows && c < self.cols);
&mut self.data[r * self.cols + c]
}
}
impl<T> Index<usize> for Matrix<T> {
type Output = [T];
fn index(&self, r: usize) -> &[T] {
assert!(r < self.rows);
&self.data[r * self.cols..(r + 1) * self.cols]
}
}
impl<T> IndexMut<usize> for Matrix<T> {
fn index_mut(&mut self, r: usize) -> &mut [T] {
assert!(r < self.rows);
&mut self.data[r * self.cols..(r + 1) * self.cols]
}
}
impl<T: fmt::Display> fmt::Display for Matrix<T> {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
for r in 0..self.rows {
if r > 0 {
writeln!(f)?;
}
for c in 0..self.cols {
if c > 0 {
write!(f, " ")?;
}
write!(f, "{}", self[(r, c)])?;
}
}
Ok(())
}
}
impl<T> Matrix<T>
where
T: Copy + Default + AddAssign + Mul<Output = T>,
{
pub fn checked_mul(&self, rhs: &Self) -> Result<Self, MatrixError> {
if self.cols != rhs.rows {
return Err(MatrixError::DimensionMismatch);
}
let mut out = Self::new(self.rows, rhs.cols, T::default());
if out.is_empty() || self.cols == 0 {
return Ok(out);
}
for i in 0..self.rows {
for k in 0..self.cols {
for j in 0..rhs.cols {
out[(i, j)] += self[(i, k)] * rhs[(k, j)];
}
}
}
Ok(out)
}
pub fn pow(&self, mut exponent: u64, one: T) -> Result<Self, MatrixError> {
if self.rows != self.cols {
return Err(MatrixError::NotSquare);
}
let mut result = Self::new(self.rows, self.cols, T::default());
for i in 0..self.rows {
result[(i, i)] = one;
}
let mut base = self.clone();
while exponent > 0 {
if exponent & 1 != 0 {
result = result.checked_mul(&base)?;
}
exponent >>= 1;
if exponent > 0 {
base = base.checked_mul(&base)?;
}
}
Ok(result)
}
}
impl<T> Mul for &Matrix<T>
where
T: Copy + Default + AddAssign + Mul<Output = T>,
{
type Output = Matrix<T>;
fn mul(self, rhs: Self) -> Self::Output {
self.checked_mul(rhs).expect("matrix dimension mismatch")
}
}
impl<T> Add for &Matrix<T>
where
T: Copy + Add<Output = T>,
{
type Output = Matrix<T>;
fn add(self, rhs: Self) -> Self::Output {
assert_eq!((self.rows, self.cols), (rhs.rows, rhs.cols));
Matrix {
rows: self.rows,
cols: self.cols,
data: self
.data
.iter()
.zip(&rhs.data)
.map(|(&a, &b)| a + b)
.collect(),
}
}
}
}
pub mod maxflow {
#![allow(dead_code, clippy::bool_to_int_with_if)]
use crate::internal_queue::SimpleQueue;
use crate::internal_type_traits::Integral;
use std::cmp::min;
use std::iter;
impl<Cap> MfGraph<Cap>
where
Cap: Integral,
{
pub fn new(n: usize) -> MfGraph<Cap> {
MfGraph {
_n: n,
pos: Vec::new(),
g: iter::repeat_with(Vec::new).take(n).collect(),
}
}
pub fn add_edge(&mut self, from: usize, to: usize, cap: Cap) -> usize {
assert!(from < self._n);
assert!(to < self._n);
assert!(Cap::zero() <= cap);
let m = self.pos.len();
self.pos.push((from, self.g[from].len()));
let rev = self.g[to].len() + usize::from(from == to);
self.g[from].push(_Edge { to, rev, cap });
let rev = self.g[from].len() - 1;
self.g[to].push(_Edge {
to: from,
rev,
cap: Cap::zero(),
});
m
}
}
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Edge<Cap: Integral> {
pub from: usize,
pub to: usize,
pub cap: Cap,
pub flow: Cap,
}
impl<Cap> MfGraph<Cap>
where
Cap: Integral,
{
pub fn get_edge(&self, i: usize) -> Edge<Cap> {
let m = self.pos.len();
assert!(i < m);
let _e = &self.g[self.pos[i].0][self.pos[i].1];
let _re = &self.g[_e.to][_e.rev];
Edge {
from: self.pos[i].0,
to: _e.to,
cap: _e.cap + _re.cap,
flow: _re.cap,
}
}
pub fn edges(&self) -> Vec<Edge<Cap>> {
let m = self.pos.len();
(0..m).map(|i| self.get_edge(i)).collect()
}
pub fn change_edge(&mut self, i: usize, new_cap: Cap, new_flow: Cap) {
let m = self.pos.len();
assert!(i < m);
assert!(Cap::zero() <= new_flow && new_flow <= new_cap);
let (to, rev) = {
let _e = &mut self.g[self.pos[i].0][self.pos[i].1];
_e.cap = new_cap - new_flow;
(_e.to, _e.rev)
};
let _re = &mut self.g[to][rev];
_re.cap = new_flow;
}
pub fn flow(&mut self, s: usize, t: usize) -> Cap {
self.flow_with_capacity(s, t, Cap::max_value())
}
pub fn flow_with_capacity(&mut self, s: usize, t: usize, flow_limit: Cap) -> Cap {
let n_ = self._n;
assert!(s < n_);
assert!(t < n_);
assert_ne!(s, t);
assert!(Cap::zero() <= flow_limit);
let mut calc = FlowCalculator {
graph: self,
s,
t,
flow_limit,
level: vec![0; n_],
iter: vec![0; n_],
que: SimpleQueue::default(),
};
let mut flow = Cap::zero();
while flow < flow_limit {
calc.bfs();
if calc.level[t] == -1 {
break;
}
calc.iter.iter_mut().for_each(|e| *e = 0);
while flow < flow_limit {
let f = calc.dfs(t, flow_limit - flow);
if f == Cap::zero() {
break;
}
flow += f;
}
}
flow
}
pub fn min_cut(&self, s: usize) -> Vec<bool> {
let mut visited = vec![false; self._n];
let mut que = SimpleQueue::default();
que.push(s);
while let Some(&p) = que.pop() {
visited[p] = true;
for e in &self.g[p] {
if e.cap != Cap::zero() && !visited[e.to] {
visited[e.to] = true;
que.push(e.to);
}
}
}
visited
}
}
struct FlowCalculator<'a, Cap> {
graph: &'a mut MfGraph<Cap>,
s: usize,
t: usize,
flow_limit: Cap,
level: Vec<i32>,
iter: Vec<usize>,
que: SimpleQueue<usize>,
}
impl<Cap> FlowCalculator<'_, Cap>
where
Cap: Integral,
{
fn bfs(&mut self) {
self.level.iter_mut().for_each(|e| *e = -1);
self.level[self.s] = 0;
self.que.clear();
self.que.push(self.s);
while let Some(&v) = self.que.pop() {
for e in &self.graph.g[v] {
if e.cap == Cap::zero() || self.level[e.to] >= 0 {
continue;
}
self.level[e.to] = self.level[v] + 1;
if e.to == self.t {
return;
}
self.que.push(e.to);
}
}
}
fn dfs(&mut self, v: usize, up: Cap) -> Cap {
if v == self.s {
return up;
}
let mut res = Cap::zero();
let level_v = self.level[v];
for i in self.iter[v]..self.graph.g[v].len() {
self.iter[v] = i;
let &_Edge {
to: e_to,
rev: e_rev,
..
} = &self.graph.g[v][i];
if level_v <= self.level[e_to] || self.graph.g[e_to][e_rev].cap == Cap::zero() {
continue;
}
let d = self.dfs(e_to, min(up - res, self.graph.g[e_to][e_rev].cap));
if d <= Cap::zero() {
continue;
}
self.graph.g[v][i].cap += d;
self.graph.g[e_to][e_rev].cap -= d;
res += d;
if res == up {
return res;
}
}
self.iter[v] = self.graph.g[v].len();
res
}
}
#[derive(Clone, Debug, Default)]
pub struct MfGraph<Cap> {
_n: usize,
pos: Vec<(usize, usize)>,
g: Vec<Vec<_Edge<Cap>>>,
}
#[derive(Clone, Debug)]
struct _Edge<Cap> {
to: usize,
rev: usize,
cap: Cap,
}
}
pub mod mex {
#![warn(missing_docs)]
use std::collections::BTreeMap;
use std::ops::{Bound, RangeBounds};
const DOMAIN_START: i128 = i64::MIN as i128;
const DOMAIN_END: i128 = i64::MAX as i128 + 1;
fn bounds<R: RangeBounds<i64>>(range: R) -> (i128, i128) {
let start = match range.start_bound() {
Bound::Included(&x) => x as i128,
Bound::Excluded(&x) => x as i128 + 1,
Bound::Unbounded => DOMAIN_START,
};
let end = match range.end_bound() {
Bound::Included(&x) => x as i128 + 1,
Bound::Excluded(&x) => x as i128,
Bound::Unbounded => DOMAIN_END,
};
assert!(start <= end, "range is reversed");
(start, end)
}
#[derive(Clone, Debug, Default)]
pub struct Mex {
intervals: BTreeMap<i128, i128>,
}
impl Mex {
#[must_use]
pub fn new() -> Self {
Self::default()
}
#[must_use]
pub fn contains(&self, value: i64) -> bool {
let x = value as i128;
self.intervals
.range(..=x)
.next_back()
.is_some_and(|(_, &end)| x < end)
}
pub fn insert(&mut self, value: i64) -> bool {
self.insert_range(value..=value)
}
pub fn remove(&mut self, value: i64) -> bool {
self.remove_range(value..=value)
}
pub fn insert_range<R: RangeBounds<i64>>(&mut self, range: R) -> bool {
let (mut start, mut end) = bounds(range);
if start == end {
return false;
}
if self.contains_all_i128(start, end) {
return false;
}
if let Some((&key, &old_end)) = self.intervals.range(..=start).next_back() {
if old_end >= start {
self.intervals.remove(&key);
start = key;
end = end.max(old_end);
}
}
let keys = self
.intervals
.range(start..=end)
.map(|(&interval_start, _)| interval_start)
.collect::<Vec<_>>();
for key in keys {
let old_end = self.intervals.remove(&key).unwrap();
end = end.max(old_end);
}
self.intervals.insert(start, end);
true
}
pub fn remove_range<R: RangeBounds<i64>>(&mut self, range: R) -> bool {
let (start, end) = bounds(range);
if start == end {
return false;
}
let mut keys = Vec::new();
if let Some((&key, &old_end)) = self.intervals.range(..start).next_back() {
if old_end > start {
keys.push(key);
}
}
keys.extend(
self.intervals
.range(start..end)
.map(|(&interval_start, _)| interval_start),
);
if keys.is_empty() {
return false;
}
for key in keys {
let old_end = self.intervals.remove(&key).unwrap();
if key < start {
self.intervals.insert(key, start);
}
if old_end > end {
self.intervals.insert(end, old_end);
}
}
true
}
#[must_use]
pub fn contains_any<R: RangeBounds<i64>>(&self, range: R) -> bool {
let (start, end) = bounds(range);
start < end
&& self
.intervals
.range(..end)
.next_back()
.is_some_and(|(_, &e)| e > start)
}
#[must_use]
pub fn contains_all<R: RangeBounds<i64>>(&self, range: R) -> bool {
let (start, end) = bounds(range);
self.contains_all_i128(start, end)
}
fn contains_all_i128(&self, start: i128, end: i128) -> bool {
start == end
|| self
.intervals
.range(..=start)
.next_back()
.is_some_and(|(_, &e)| e >= end)
}
#[must_use]
pub fn mex<R: RangeBounds<i64>>(&self, range: R) -> Option<i64> {
let (start, end) = bounds(range);
if start == end {
return None;
}
let candidate = self
.intervals
.range(..=start)
.next_back()
.map_or(
start,
|(&s, &e)| if s <= start && start < e { e } else { start },
);
(candidate < end).then_some(candidate as i64)
}
}
}
pub mod mincostflow {
use crate::internal_type_traits::Integral;
#[derive(Clone, Debug)]
pub struct MinCostFlowEdge<T> {
pub from: usize,
pub to: usize,
pub cap: T,
pub flow: T,
pub cost: T,
}
#[derive(Clone, Debug)]
pub struct MinCostFlowGraph<T> {
pos: Vec<(usize, usize)>,
g: Vec<Vec<_Edge<T>>>,
}
#[allow(clippy::redundant_closure)]
impl<T> MinCostFlowGraph<T>
where
T: Integral + std::ops::Neg<Output = T>,
{
pub fn new(n: usize) -> Self {
Self {
pos: vec![],
g: (0..n).map(|_| vec![]).collect(),
}
}
pub fn get_edge(&self, i: usize) -> MinCostFlowEdge<T> {
assert!(i < self.pos.len());
let e = &self.g[self.pos[i].0][self.pos[i].1];
let re = &self.g[e.to][e.rev];
MinCostFlowEdge {
from: self.pos[i].0,
to: e.to,
cap: e.cap + re.cap,
flow: re.cap,
cost: e.cost,
}
}
pub fn edges(&self) -> Vec<MinCostFlowEdge<T>> {
let m = self.pos.len();
let mut result = vec![];
for i in 0..m {
result.push(self.get_edge(i));
}
result
}
pub fn add_edge(&mut self, from: usize, to: usize, cap: T, cost: T) -> usize {
assert!(from < self.g.len());
assert!(to < self.g.len());
assert_ne!(from, to);
assert!(cap >= T::zero());
assert!(cost >= T::zero());
self.pos.push((from, self.g[from].len()));
let rev = self.g[to].len();
self.g[from].push(_Edge { to, rev, cap, cost });
let rev = self.g[from].len() - 1;
self.g[to].push(_Edge {
to: from,
rev,
cap: T::zero(),
cost: -cost,
});
self.pos.len() - 1
}
pub fn flow(&mut self, source: usize, sink: usize, flow_limit: T) -> (T, T) {
self.slope(source, sink, flow_limit).pop().unwrap()
}
pub fn slope(&mut self, source: usize, sink: usize, flow_limit: T) -> Vec<(T, T)> {
let n = self.g.len();
assert!(source < n);
assert!(sink < n);
assert_ne!(source, sink);
let mut dual = vec![T::zero(); n];
let mut prev_v = vec![0; n];
let mut prev_e = vec![0; n];
let mut flow = T::zero();
let mut cost = T::zero();
let mut prev_cost_per_flow: Option<T> = None;
let mut result = vec![(flow, cost)];
while flow < flow_limit {
if !self.refine_dual(source, sink, &mut dual, &mut prev_v, &mut prev_e) {
break;
}
let mut c = flow_limit - flow;
let mut v = sink;
while v != source {
c = std::cmp::min(c, self.g[prev_v[v]][prev_e[v]].cap);
v = prev_v[v];
}
let mut v = sink;
while v != source {
self.g[prev_v[v]][prev_e[v]].cap -= c;
let rev = self.g[prev_v[v]][prev_e[v]].rev;
self.g[v][rev].cap += c;
v = prev_v[v];
}
let d = -dual[source];
flow += c;
cost += d * c;
if prev_cost_per_flow == Some(d) {
assert!(result.pop().is_some());
}
result.push((flow, cost));
prev_cost_per_flow = Some(d);
}
result
}
fn refine_dual(
&self,
source: usize,
sink: usize,
dual: &mut [T],
pv: &mut [usize],
pe: &mut [usize],
) -> bool {
let n = self.g.len();
let mut dist = vec![T::max_value(); n];
let mut vis = vec![false; n];
let mut que = std::collections::BinaryHeap::new();
dist[source] = T::zero();
que.push((std::cmp::Reverse(T::zero()), source));
while let Some((_, v)) = que.pop() {
if vis[v] {
continue;
}
vis[v] = true;
if v == sink {
break;
}
for (i, e) in self.g[v].iter().enumerate() {
if vis[e.to] || e.cap == T::zero() {
continue;
}
let cost = e.cost - dual[e.to] + dual[v];
if dist[e.to] - dist[v] > cost {
dist[e.to] = dist[v] + cost;
pv[e.to] = v;
pe[e.to] = i;
que.push((std::cmp::Reverse(dist[e.to]), e.to));
}
}
}
if !vis[sink] {
return false;
}
for v in 0..n {
if !vis[v] {
continue;
}
dual[v] -= dist[sink] - dist[v];
}
true
}
}
#[derive(Clone, Debug)]
struct _Edge<T> {
to: usize,
rev: usize,
cap: T,
cost: T,
}
}
pub mod modint {
use crate::internal_math;
use std::{
cell::RefCell,
convert::{Infallible, TryInto as _},
fmt,
hash::{Hash, Hasher},
iter::{Product, Sum},
marker::PhantomData,
ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign},
str::FromStr,
sync::atomic::{self, AtomicU32, AtomicU64},
thread::LocalKey,
};
pub type ModInt1000000007 = StaticModInt<Mod1000000007>;
pub type ModInt998244353 = StaticModInt<Mod998244353>;
pub type ModInt = DynamicModInt<DefaultId>;
#[derive(Copy, Clone, Eq, PartialEq)]
#[repr(transparent)]
pub struct StaticModInt<M> {
val: u32,
phantom: PhantomData<fn() -> M>,
}
impl<M: Modulus> StaticModInt<M> {
#[inline(always)]
pub fn modulus() -> u32 {
M::VALUE
}
#[inline]
pub fn new<T: RemEuclidU32>(val: T) -> Self {
Self::raw(val.rem_euclid_u32(M::VALUE))
}
#[inline]
pub fn raw(val: u32) -> Self {
Self {
val,
phantom: PhantomData,
}
}
#[inline]
pub fn val(self) -> u32 {
self.val
}
#[inline]
pub fn pow(self, n: u64) -> Self {
<Self as ModIntBase>::pow(self, n)
}
#[inline]
pub fn inv(self) -> Self {
if M::HINT_VALUE_IS_PRIME {
if self.val() == 0 {
panic!("attempt to divide by zero");
}
debug_assert!(
internal_math::is_prime(M::VALUE.try_into().unwrap()),
"{} is not a prime number",
M::VALUE,
);
self.pow((M::VALUE - 2).into())
} else {
Self::inv_for_non_prime_modulus(self)
}
}
}
impl<M: Modulus> ModIntBase for StaticModInt<M> {
#[inline(always)]
fn modulus() -> u32 {
Self::modulus()
}
#[inline]
fn raw(val: u32) -> Self {
Self::raw(val)
}
#[inline]
fn val(self) -> u32 {
self.val()
}
#[inline]
fn inv(self) -> Self {
self.inv()
}
}
pub trait Modulus: 'static + Copy + Eq {
const VALUE: u32;
const HINT_VALUE_IS_PRIME: bool;
fn butterfly_cache() -> &'static LocalKey<RefCell<Option<ButterflyCache<Self>>>>;
}
#[derive(Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Hash, Debug)]
pub enum Mod1000000007 {}
impl Modulus for Mod1000000007 {
const VALUE: u32 = 1_000_000_007;
const HINT_VALUE_IS_PRIME: bool = true;
fn butterfly_cache() -> &'static LocalKey<RefCell<Option<ButterflyCache<Self>>>> {
thread_local! {
static BUTTERFLY_CACHE: RefCell<Option<ButterflyCache<Mod1000000007>>> = RefCell::default();
}
&BUTTERFLY_CACHE
}
}
#[derive(Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Hash, Debug)]
pub enum Mod998244353 {}
impl Modulus for Mod998244353 {
const VALUE: u32 = 998_244_353;
const HINT_VALUE_IS_PRIME: bool = true;
fn butterfly_cache() -> &'static LocalKey<RefCell<Option<ButterflyCache<Self>>>> {
thread_local! {
static BUTTERFLY_CACHE: RefCell<Option<ButterflyCache<Mod998244353>>> = RefCell::default();
}
&BUTTERFLY_CACHE
}
}
pub struct ButterflyCache<M> {
pub(crate) sum_e: Vec<StaticModInt<M>>,
pub(crate) sum_ie: Vec<StaticModInt<M>>,
}
#[derive(Copy, Clone, Eq, PartialEq)]
#[repr(transparent)]
pub struct DynamicModInt<I> {
val: u32,
phantom: PhantomData<fn() -> I>,
}
impl<I: Id> DynamicModInt<I> {
#[inline]
pub fn modulus() -> u32 {
I::companion_barrett().umod()
}
#[inline]
pub fn set_modulus(modulus: u32) {
if modulus == 0 {
panic!("the modulus must not be 0");
}
I::companion_barrett().update(modulus);
}
#[inline]
pub fn new<T: RemEuclidU32>(val: T) -> Self {
<Self as ModIntBase>::new(val)
}
#[inline]
pub fn raw(val: u32) -> Self {
Self {
val,
phantom: PhantomData,
}
}
#[inline]
pub fn val(self) -> u32 {
self.val
}
#[inline]
pub fn pow(self, n: u64) -> Self {
<Self as ModIntBase>::pow(self, n)
}
#[inline]
pub fn inv(self) -> Self {
Self::inv_for_non_prime_modulus(self)
}
}
impl<I: Id> ModIntBase for DynamicModInt<I> {
#[inline]
fn modulus() -> u32 {
Self::modulus()
}
#[inline]
fn raw(val: u32) -> Self {
Self::raw(val)
}
#[inline]
fn val(self) -> u32 {
self.val()
}
#[inline]
fn inv(self) -> Self {
self.inv()
}
}
pub trait Id: 'static + Copy + Eq {
fn companion_barrett() -> &'static Barrett;
}
#[derive(Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Hash, Debug)]
pub enum DefaultId {}
impl Id for DefaultId {
fn companion_barrett() -> &'static Barrett {
static BARRETT: Barrett = Barrett::default();
&BARRETT
}
}
pub struct Barrett {
m: AtomicU32,
im: AtomicU64,
}
impl Barrett {
#[inline]
pub const fn new(m: u32) -> Self {
Self {
m: AtomicU32::new(m),
im: AtomicU64::new((-1i64 as u64 / m as u64).wrapping_add(1)),
}
}
#[inline]
const fn default() -> Self {
Self::new(998_244_353)
}
#[inline]
fn update(&self, m: u32) {
let im = (-1i64 as u64 / m as u64).wrapping_add(1);
self.m.store(m, atomic::Ordering::SeqCst);
self.im.store(im, atomic::Ordering::SeqCst);
}
#[inline]
fn umod(&self) -> u32 {
self.m.load(atomic::Ordering::SeqCst)
}
#[inline]
fn mul(&self, a: u32, b: u32) -> u32 {
let m = self.m.load(atomic::Ordering::SeqCst);
let im = self.im.load(atomic::Ordering::SeqCst);
internal_math::mul_mod(a, b, m, im)
}
}
impl Default for Barrett {
#[inline]
fn default() -> Self {
Self::default()
}
}
pub trait ModIntBase:
Default
+ FromStr
+ From<i8>
+ From<i16>
+ From<i32>
+ From<i64>
+ From<i128>
+ From<isize>
+ From<u8>
+ From<u16>
+ From<u32>
+ From<u64>
+ From<u128>
+ From<usize>
+ Copy
+ Eq
+ Hash
+ fmt::Display
+ fmt::Debug
+ Neg<Output = Self>
+ Add<Output = Self>
+ Sub<Output = Self>
+ Mul<Output = Self>
+ Div<Output = Self>
+ AddAssign
+ SubAssign
+ MulAssign
+ DivAssign
{
fn modulus() -> u32;
fn raw(val: u32) -> Self;
fn val(self) -> u32;
fn inv(self) -> Self;
#[inline]
fn new<T: RemEuclidU32>(val: T) -> Self {
Self::raw(val.rem_euclid_u32(Self::modulus()))
}
#[inline]
fn pow(self, mut n: u64) -> Self {
let mut x = self;
let mut r = Self::raw(u32::from(Self::modulus() > 1));
while n > 0 {
if n & 1 == 1 {
r *= x;
}
x *= x;
n >>= 1;
}
r
}
}
pub trait RemEuclidU32 {
fn rem_euclid_u32(self, modulus: u32) -> u32;
}
macro_rules! impl_rem_euclid_u32_for_small_signed {
($($ty:tt),*) => {
$(
impl RemEuclidU32 for $ty {
#[inline]
fn rem_euclid_u32(self, modulus: u32) -> u32 {
(self as i64).rem_euclid(i64::from(modulus)) as _
}
}
)*
}
}
impl_rem_euclid_u32_for_small_signed!(i8, i16, i32, i64, isize);
impl RemEuclidU32 for i128 {
#[inline]
fn rem_euclid_u32(self, modulus: u32) -> u32 {
self.rem_euclid(i128::from(modulus)) as _
}
}
macro_rules! impl_rem_euclid_u32_for_small_unsigned {
($($ty:tt),*) => {
$(
impl RemEuclidU32 for $ty {
#[inline]
fn rem_euclid_u32(self, modulus: u32) -> u32 {
self as u32 % modulus
}
}
)*
}
}
macro_rules! impl_rem_euclid_u32_for_large_unsigned {
($($ty:tt),*) => {
$(
impl RemEuclidU32 for $ty {
#[inline]
fn rem_euclid_u32(self, modulus: u32) -> u32 {
(self % (modulus as $ty)) as _
}
}
)*
}
}
impl_rem_euclid_u32_for_small_unsigned!(u8, u16, u32);
impl_rem_euclid_u32_for_large_unsigned!(u64, u128);
#[cfg(target_pointer_width = "32")]
impl_rem_euclid_u32_for_small_unsigned!(usize);
#[cfg(target_pointer_width = "64")]
impl_rem_euclid_u32_for_large_unsigned!(usize);
trait InternalImplementations: ModIntBase {
#[inline]
fn inv_for_non_prime_modulus(this: Self) -> Self {
let (gcd, x) = internal_math::inv_gcd(this.val().into(), Self::modulus().into());
if gcd != 1 {
panic!("the multiplicative inverse does not exist");
}
Self::new(x)
}
#[inline]
fn default_impl() -> Self {
Self::raw(0)
}
#[inline]
fn from_str_impl(s: &str) -> Result<Self, Infallible> {
Ok(s.parse::<i64>()
.map(Self::new)
.unwrap_or_else(|_| todo!("parsing as an arbitrary precision integer?")))
}
#[inline]
fn hash_impl(this: &Self, state: &mut impl Hasher) {
this.val().hash(state)
}
#[inline]
fn display_impl(this: &Self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
fmt::Display::fmt(&this.val(), f)
}
#[inline]
fn debug_impl(this: &Self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
fmt::Debug::fmt(&this.val(), f)
}
#[inline]
fn neg_impl(this: Self) -> Self {
Self::sub_impl(Self::raw(0), this)
}
#[inline]
fn add_impl(lhs: Self, rhs: Self) -> Self {
let modulus = Self::modulus();
let mut val = lhs.val() + rhs.val();
if val >= modulus {
val -= modulus;
}
Self::raw(val)
}
#[inline]
fn sub_impl(lhs: Self, rhs: Self) -> Self {
let modulus = Self::modulus();
let mut val = lhs.val().wrapping_sub(rhs.val());
if val >= modulus {
val = val.wrapping_add(modulus)
}
Self::raw(val)
}
fn mul_impl(lhs: Self, rhs: Self) -> Self;
#[inline]
fn div_impl(lhs: Self, rhs: Self) -> Self {
Self::mul_impl(lhs, rhs.inv())
}
}
impl<M: Modulus> InternalImplementations for StaticModInt<M> {
#[inline]
fn mul_impl(lhs: Self, rhs: Self) -> Self {
Self::raw((u64::from(lhs.val()) * u64::from(rhs.val()) % u64::from(M::VALUE)) as u32)
}
}
impl<I: Id> InternalImplementations for DynamicModInt<I> {
#[inline]
fn mul_impl(lhs: Self, rhs: Self) -> Self {
Self::raw(I::companion_barrett().mul(lhs.val, rhs.val))
}
}
macro_rules! impl_basic_traits {
() => {};
(impl <$generic_param:ident : $generic_param_bound:tt> _ for $self:ty; $($rest:tt)*) => {
impl <$generic_param: $generic_param_bound> Default for $self {
#[inline]
fn default() -> Self {
Self::default_impl()
}
}
impl <$generic_param: $generic_param_bound> FromStr for $self {
type Err = Infallible;
#[inline]
fn from_str(s: &str) -> Result<Self, Infallible> {
Self::from_str_impl(s)
}
}
impl<$generic_param: $generic_param_bound, V: RemEuclidU32> From<V> for $self {
#[inline]
fn from(from: V) -> Self {
Self::new(from)
}
}
#[allow(clippy::derived_hash_with_manual_eq)]
impl<$generic_param: $generic_param_bound> Hash for $self {
#[inline]
fn hash<H: Hasher>(&self, state: &mut H) {
Self::hash_impl(self, state)
}
}
impl<$generic_param: $generic_param_bound> fmt::Display for $self {
#[inline]
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
Self::display_impl(self, f)
}
}
impl<$generic_param: $generic_param_bound> fmt::Debug for $self {
#[inline]
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
Self::debug_impl(self, f)
}
}
impl<$generic_param: $generic_param_bound> Neg for $self {
type Output = $self;
#[inline]
fn neg(self) -> $self {
Self::neg_impl(self)
}
}
impl<$generic_param: $generic_param_bound> Neg for &'_ $self {
type Output = $self;
#[inline]
fn neg(self) -> $self {
<$self>::neg_impl(*self)
}
}
impl_basic_traits!($($rest)*);
};
}
impl_basic_traits! {
impl <M: Modulus> _ for StaticModInt<M> ;
impl <I: Id > _ for DynamicModInt<I>;
}
macro_rules! impl_bin_ops {
() => {};
(for<$($generic_param:ident : $generic_param_bound:tt),*> <$lhs_ty:ty> ~ <$rhs_ty:ty> -> $output:ty { { $lhs_body:expr } ~ { $rhs_body:expr } } $($rest:tt)*) => {
impl <$($generic_param: $generic_param_bound),*> Add<$rhs_ty> for $lhs_ty {
type Output = $output;
#[inline]
fn add(self, rhs: $rhs_ty) -> $output {
<$output>::add_impl(apply($lhs_body, self), apply($rhs_body, rhs))
}
}
impl <$($generic_param: $generic_param_bound),*> Sub<$rhs_ty> for $lhs_ty {
type Output = $output;
#[inline]
fn sub(self, rhs: $rhs_ty) -> $output {
<$output>::sub_impl(apply($lhs_body, self), apply($rhs_body, rhs))
}
}
impl <$($generic_param: $generic_param_bound),*> Mul<$rhs_ty> for $lhs_ty {
type Output = $output;
#[inline]
fn mul(self, rhs: $rhs_ty) -> $output {
<$output>::mul_impl(apply($lhs_body, self), apply($rhs_body, rhs))
}
}
impl <$($generic_param: $generic_param_bound),*> Div<$rhs_ty> for $lhs_ty {
type Output = $output;
#[inline]
fn div(self, rhs: $rhs_ty) -> $output {
<$output>::div_impl(apply($lhs_body, self), apply($rhs_body, rhs))
}
}
impl_bin_ops!($($rest)*);
};
}
macro_rules! impl_assign_ops {
() => {};
(for<$($generic_param:ident : $generic_param_bound:tt),*> <$lhs_ty:ty> ~= <$rhs_ty:ty> { _ ~= { $rhs_body:expr } } $($rest:tt)*) => {
impl <$($generic_param: $generic_param_bound),*> AddAssign<$rhs_ty> for $lhs_ty {
#[inline]
fn add_assign(&mut self, rhs: $rhs_ty) {
*self = *self + apply($rhs_body, rhs);
}
}
impl <$($generic_param: $generic_param_bound),*> SubAssign<$rhs_ty> for $lhs_ty {
#[inline]
fn sub_assign(&mut self, rhs: $rhs_ty) {
*self = *self - apply($rhs_body, rhs);
}
}
impl <$($generic_param: $generic_param_bound),*> MulAssign<$rhs_ty> for $lhs_ty {
#[inline]
fn mul_assign(&mut self, rhs: $rhs_ty) {
*self = *self * apply($rhs_body, rhs);
}
}
impl <$($generic_param: $generic_param_bound),*> DivAssign<$rhs_ty> for $lhs_ty {
#[inline]
fn div_assign(&mut self, rhs: $rhs_ty) {
*self = *self / apply($rhs_body, rhs);
}
}
impl_assign_ops!($($rest)*);
};
}
#[inline]
fn apply<F: FnOnce(X) -> O, X, O>(f: F, x: X) -> O {
f(x)
}
impl_bin_ops! {
for<M: Modulus> <StaticModInt<M> > ~ <StaticModInt<M> > -> StaticModInt<M> { { |x| x } ~ { |x| x } }
for<M: Modulus> <StaticModInt<M> > ~ <&'_ StaticModInt<M> > -> StaticModInt<M> { { |x| x } ~ { |&x| x } }
for<M: Modulus> <&'_ StaticModInt<M> > ~ <StaticModInt<M> > -> StaticModInt<M> { { |&x| x } ~ { |x| x } }
for<M: Modulus> <&'_ StaticModInt<M> > ~ <&'_ StaticModInt<M> > -> StaticModInt<M> { { |&x| x } ~ { |&x| x } }
for<I: Id > <DynamicModInt<I> > ~ <DynamicModInt<I> > -> DynamicModInt<I> { { |x| x } ~ { |x| x } }
for<I: Id > <DynamicModInt<I> > ~ <&'_ DynamicModInt<I>> -> DynamicModInt<I> { { |x| x } ~ { |&x| x } }
for<I: Id > <&'_ DynamicModInt<I>> ~ <DynamicModInt<I> > -> DynamicModInt<I> { { |&x| x } ~ { |x| x } }
for<I: Id > <&'_ DynamicModInt<I>> ~ <&'_ DynamicModInt<I>> -> DynamicModInt<I> { { |&x| x } ~ { |&x| x } }
for<M: Modulus, T: RemEuclidU32> <StaticModInt<M> > ~ <T> -> StaticModInt<M> { { |x| x } ~ { StaticModInt::<M>::new } }
for<I: Id , T: RemEuclidU32> <DynamicModInt<I> > ~ <T> -> DynamicModInt<I> { { |x| x } ~ { DynamicModInt::<I>::new } }
}
impl_assign_ops! {
for<M: Modulus> <StaticModInt<M> > ~= <StaticModInt<M> > { _ ~= { |x| x } }
for<M: Modulus> <StaticModInt<M> > ~= <&'_ StaticModInt<M> > { _ ~= { |&x| x } }
for<I: Id > <DynamicModInt<I>> ~= <DynamicModInt<I> > { _ ~= { |x| x } }
for<I: Id > <DynamicModInt<I>> ~= <&'_ DynamicModInt<I>> { _ ~= { |&x| x } }
for<M: Modulus, T: RemEuclidU32> <StaticModInt<M> > ~= <T> { _ ~= { StaticModInt::<M>::new } }
for<I: Id, T: RemEuclidU32> <DynamicModInt<I>> ~= <T> { _ ~= { DynamicModInt::<I>::new } }
}
macro_rules! impl_folding {
() => {};
(impl<$generic_param:ident : $generic_param_bound:tt> $trait:ident<_> for $self:ty { fn $method:ident(_) -> _ { _($unit:expr, $op:expr) } } $($rest:tt)*) => {
impl<$generic_param: $generic_param_bound> $trait<Self> for $self {
#[inline]
fn $method<S>(iter: S) -> Self
where
S: Iterator<Item = Self>,
{
iter.fold($unit, $op)
}
}
impl<'a, $generic_param: $generic_param_bound> $trait<&'a Self> for $self {
#[inline]
fn $method<S>(iter: S) -> Self
where
S: Iterator<Item = &'a Self>,
{
iter.fold($unit, $op)
}
}
impl_folding!($($rest)*);
};
}
impl_folding! {
impl<M: Modulus> Sum<_> for StaticModInt<M> { fn sum(_) -> _ { _(Self::raw(0), Add::add) } }
impl<M: Modulus> Product<_> for StaticModInt<M> { fn product(_) -> _ { _(Self::raw(u32::from(Self::modulus() > 1)), Mul::mul) } }
impl<I: Id > Sum<_> for DynamicModInt<I> { fn sum(_) -> _ { _(Self::raw(0), Add::add) } }
impl<I: Id > Product<_> for DynamicModInt<I> { fn product(_) -> _ { _(Self::raw(u32::from(Self::modulus() > 1)), Mul::mul) } }
}
}
pub mod ordfloat {
#![warn(missing_docs)]
use std::cmp::Ordering;
use std::hash::{Hash, Hasher};
use std::ops::Deref;
#[derive(Copy, Clone, Debug, Default)]
pub struct OrdFloat(pub f64);
impl OrdFloat {
#[must_use]
pub const fn into_inner(self) -> f64 {
self.0
}
}
impl From<f64> for OrdFloat {
fn from(value: f64) -> Self {
Self(value)
}
}
impl PartialEq for OrdFloat {
fn eq(&self, other: &Self) -> bool {
self.cmp(other).is_eq()
}
}
impl Eq for OrdFloat {}
impl Ord for OrdFloat {
fn cmp(&self, other: &Self) -> Ordering {
self.0.total_cmp(&other.0)
}
}
impl PartialOrd for OrdFloat {
fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
Some(self.cmp(other))
}
}
impl Hash for OrdFloat {
fn hash<H: Hasher>(&self, state: &mut H) {
self.0.to_bits().hash(state);
}
}
impl Deref for OrdFloat {
type Target = f64;
fn deref(&self) -> &Self::Target {
&self.0
}
}
}
pub mod pcg64 {
#![warn(missing_docs)]
use crate::SmallRng;
use std::time::{SystemTime, SystemTimeError, UNIX_EPOCH};
const MULTIPLIER: u128 = 47_026_247_687_942_121_848_144_207_491_837_523_525;
const DEFAULT_STREAM: u128 = 1_442_695_040_888_963_407;
const DEFAULT_SEED: u128 = 88_172_645_463_325_252;
#[derive(Clone, Debug)]
pub struct PCG64 {
seed: u128,
state: u128,
increment: u128,
}
impl Default for PCG64 {
fn default() -> Self {
Self::seed_from_u128(DEFAULT_SEED)
}
}
impl PCG64 {
#[must_use]
pub fn seed_from_u128(seed: u128) -> Self {
Self::with_stream(seed, DEFAULT_STREAM)
}
#[must_use]
pub fn with_stream(seed: u128, stream: u128) -> Self {
let mut rng = Self {
seed,
state: 0,
increment: (stream << 1) | 1,
};
let _ = rng.next_u64();
rng.state = rng.state.wrapping_add(seed);
let _ = rng.next_u64();
rng
}
pub fn from_system_time() -> Result<Self, SystemTimeError> {
Ok(Self::seed_from_u128(
SystemTime::now().duration_since(UNIX_EPOCH)?.as_nanos(),
))
}
#[must_use]
pub fn seed(&self) -> u128 {
self.seed
}
}
impl SmallRng for PCG64 {
fn next_u64(&mut self) -> u64 {
self.state = self
.state
.wrapping_mul(MULTIPLIER)
.wrapping_add(self.increment);
(((self.state >> 64) ^ self.state) as u64).rotate_right((self.state >> 122) as u32)
}
}
}
pub mod primenumber {
#![warn(missing_docs)]
use std::convert::TryFrom;
use std::iter::Iterator;
type ValueType = u64;
pub trait GetDivisor {
fn get_divisor(&self) -> Divisor;
}
macro_rules! GetDivisor_macro{
($($t:ty),*) => {
$(
impl GetDivisor for $t {
fn get_divisor(&self) -> Divisor {
Divisor::calc(
ValueType::try_from(*self)
.expect("the value must be a non-negative u64"),
)
}
})*
};
}
GetDivisor_macro!(u32, u64, u128, usize, i32, i64, i128, isize);
pub trait GetPrimeFactorization {
fn prime_factorization(&self) -> PrimeFactorization<'_>;
}
macro_rules! PrimeFactorization_macro{
($($t:ty),*) => {
$(
impl GetPrimeFactorization for $t {
fn prime_factorization(&self) -> PrimeFactorization<'_> {
PrimeFactorization::calc(
ValueType::try_from(*self)
.expect("the value must be a non-negative u64"),
)
}
})*
};
}
PrimeFactorization_macro!(u32, u64, u128, usize, i32, i64, i128, isize);
pub struct Divisor {
n: ValueType,
cur: ValueType,
flag: bool,
}
impl Divisor {
pub fn calc(n: ValueType) -> Divisor {
Divisor {
n,
cur: 1,
flag: false,
}
}
}
impl Iterator for Divisor {
type Item = ValueType;
fn next(&mut self) -> Option<Self::Item> {
if self.cur > self.n / self.cur {
None
} else if self.flag {
if self.cur == self.n / self.cur {
return None;
}
self.flag = false;
self.cur += 1;
Some(self.n / (self.cur - 1))
} else {
while !self.n.is_multiple_of(self.cur) {
self.cur += 1;
if self.cur > self.n / self.cur {
return None;
}
}
self.flag = true;
Some(self.cur)
}
}
}
pub struct PrimeFactorization<'a> {
n: ValueType,
cur: ValueType,
p_list: Option<&'a [ValueType]>,
idx: usize,
}
impl<'a> PrimeFactorization<'a> {
pub fn calc(n: ValueType) -> PrimeFactorization<'a> {
PrimeFactorization {
n,
cur: 1,
p_list: None,
idx: 0,
}
}
pub fn calc_fast(n: ValueType, p_list: &'a [ValueType]) -> PrimeFactorization<'a> {
PrimeFactorization {
n,
cur: 1,
p_list: Some(p_list),
idx: 0,
}
}
}
impl Iterator for PrimeFactorization<'_> {
type Item = ValueType;
fn next(&mut self) -> Option<Self::Item> {
loop {
if self.cur == 0 || self.cur > self.n {
return None;
}
if let Some(p_list) = self.p_list {
if self.idx >= p_list.len() {
self.p_list = None;
self.cur = if self.cur < 2 {
2
} else if self.cur == 2 {
3
} else if self.cur.is_multiple_of(2) {
self.cur.checked_add(1)?
} else {
self.cur.checked_add(2)?
};
} else {
self.cur = p_list[self.idx];
self.idx += 1;
}
} else {
self.cur = if self.cur < 2 {
2
} else if self.cur == 2 {
3
} else {
self.cur.checked_add(2)?
};
}
if self.cur == 0 {
return None;
}
if self.cur > self.n / self.cur {
if self.n != 1 {
self.cur = 0;
return Some(self.n);
}
return None;
}
if self.n.is_multiple_of(self.cur) {
self.n /= self.cur;
if self.p_list.is_some() {
self.idx -= 1;
return Some(self.cur);
}
let factor = self.cur;
self.cur = if factor <= 3 { 1 } else { factor - 2 };
return Some(factor);
}
}
}
}
pub fn get_primelist(u: ValueType) -> Vec<ValueType> {
let u = usize::try_from(u).expect("the upper bound must fit in usize");
let len = u.checked_add(1).expect("the upper bound is too large");
let mut composite = vec![false; len];
let mut primes = Vec::new();
for value in 2..=u {
if !composite[value] {
primes.push(value);
}
for &prime in &primes {
if prime > u / value {
break;
}
composite[value * prime] = true;
if value.is_multiple_of(prime) {
break;
}
}
}
primes.into_iter().map(|value| value as u64).collect()
}
pub fn get_mobius(n: ValueType) -> Vec<isize> {
let n_usize = usize::try_from(n).expect("the upper bound must fit in usize");
let len = n_usize
.checked_add(1)
.expect("the upper bound is too large");
let mut mobius = vec![0; len];
if n_usize == 0 {
return mobius;
}
mobius[1] = 1;
let mut primes = Vec::new();
let mut composite = vec![false; len];
for value in 2..=n_usize {
if !composite[value] {
primes.push(value);
mobius[value] = -1;
}
for &prime in &primes {
if prime > n_usize / value {
break;
}
let product = value * prime;
composite[product] = true;
if value.is_multiple_of(prime) {
mobius[product] = 0;
break;
}
mobius[product] = -mobius[value];
}
}
mobius
}
fn modpow_128bit(mut s: u128, mut n: u128, p: u128) -> u128 {
if p == 0 {
return 1;
}
let mut t = s;
s = 1;
while n > 0 {
if n & 1 != 0 {
s *= t;
s %= p;
}
n >>= 1;
t *= t;
t %= p;
}
s
}
pub fn miller_rabin(n: u64) -> bool {
if n == 2 {
return true;
}
if n == 1 || n.is_multiple_of(2) {
return false;
}
let (mut s, mut t) = (0, n - 1);
while t % 2 == 0 {
s += 1;
t >>= 1;
}
let bases: &[u64] = if n < 4_759_123_141 {
&[2, 7, 61]
} else if n < 341_550_071_728_321 {
&[2, 3, 5, 7, 11, 13, 17]
} else if n < 3_825_123_056_546_413_051 {
&[2, 3, 5, 7, 11, 13, 17, 19, 23]
} else {
&[2, 325, 9_375, 28_178, 450_775, 9_780_504, 1_795_265_022]
};
'witness: for &base in bases {
if base >= n {
continue;
}
let mut value = modpow_128bit(base as u128, t as u128, n as u128) as u64;
if value == 1 || value == n - 1 {
continue;
}
for _ in 1..s {
value = (u128::from(value) * u128::from(value) % u128::from(n)) as u64;
if value == n - 1 {
continue 'witness;
}
}
return false;
}
true
}
fn gcd_u64(a: u64, b: u64) -> u64
where
{
if b == 0 {
return a;
}
gcd_u64(b, a % b)
}
pub struct PollardRho {
arr: Vec<u64>,
}
impl PollardRho {
pub fn calc(n: u64) -> PollardRho {
PollardRho { arr: vec![n] }
}
}
fn brent_factor(n: u64) -> u64 {
if n.is_multiple_of(2) {
return 2;
}
const BATCH: u64 = 128;
for constant in 1_u64.. {
let step = |x: u64| {
((u128::from(x) * u128::from(x) + u128::from(constant)) % u128::from(n)) as u64
};
let (mut y, mut cycle, mut divisor) = (2_u64, 1_u64, 1_u64);
let (mut x, mut saved_y) = (0_u64, 0_u64);
while divisor == 1 {
x = y;
for _ in 0..cycle {
y = step(y);
}
let mut offset = 0;
while offset < cycle && divisor == 1 {
saved_y = y;
let mut product = 1_u64;
for _ in 0..BATCH.min(cycle - offset) {
y = step(y);
product = ((u128::from(product) * u128::from(x.abs_diff(y)))
% u128::from(n)) as u64;
}
divisor = gcd_u64(product, n);
offset += BATCH;
}
cycle = cycle.saturating_mul(2);
}
if divisor == n {
loop {
saved_y = step(saved_y);
divisor = gcd_u64(x.abs_diff(saved_y), n);
if divisor > 1 {
break;
}
}
}
if divisor != n {
return divisor;
}
}
unreachable!()
}
impl Iterator for PollardRho {
type Item = ValueType;
fn next(&mut self) -> Option<Self::Item> {
if self.arr.is_empty() || self.arr[0] == 0 {
return None;
}
let n = self.arr.pop().unwrap();
if n == 1 {
return None;
}
if miller_rabin(n) {
let r = n;
return Some(r);
}
if n.is_multiple_of(2) {
self.arr.push(n / 2);
return Some(2);
}
let divisor = brent_factor(n);
self.arr.push(n / divisor);
self.arr.push(divisor);
self.next()
}
}
}
pub mod raq_bit {
#![warn(missing_docs)]
use crate::FenwickTree;
use std::convert::TryFrom;
use std::ops::{Bound, RangeBounds};
#[derive(Clone)]
pub struct RaqBit {
n: usize,
bit0: FenwickTree<i128>,
bit1: FenwickTree<i128>,
}
impl RaqBit {
#[must_use]
pub fn new(n: usize) -> Self {
Self {
n,
bit0: FenwickTree::new(n, 0),
bit1: FenwickTree::new(n, 0),
}
}
#[must_use]
pub fn len(&self) -> usize {
self.n
}
#[must_use]
pub fn is_empty(&self) -> bool {
self.n == 0
}
pub fn prefix_sum(&self, end: usize) -> i64 {
i64::try_from(self.prefix_sum_i128(end)).expect("range sum overflow")
}
fn prefix_sum_i128(&self, end: usize) -> i128 {
assert!(end <= self.n);
let index = i128::try_from(end).expect("index exceeds i128");
self.bit0
.accum(end)
.checked_add(
self.bit1
.accum(end)
.checked_mul(index)
.expect("range sum overflow"),
)
.expect("range sum overflow")
}
pub fn add<R: RangeBounds<usize>>(&mut self, range: R, value: i64) {
let (l, r) = self.bounds(range);
if l == r {
return;
}
let li = i128::try_from(l).expect("index exceeds i128");
let value = i128::from(value);
self.bit0.add(
l,
value
.checked_mul(li)
.and_then(|x| x.checked_neg())
.expect("range update overflow"),
);
self.bit1.add(l, value);
if r < self.n {
let ri = i128::try_from(r).expect("index exceeds i128");
self.bit0
.add(r, value.checked_mul(ri).expect("range update overflow"));
self.bit1
.add(r, value.checked_neg().expect("range update overflow"));
}
}
pub fn sum<R: RangeBounds<usize>>(&self, range: R) -> i64 {
let (l, r) = self.bounds(range);
let sum = self
.prefix_sum_i128(r)
.checked_sub(self.prefix_sum_i128(l))
.expect("range sum overflow");
i64::try_from(sum).expect("range sum overflow")
}
pub fn get(&self, index: usize) -> i64 {
assert!(index < self.n);
self.sum(index..=index)
}
fn bounds<R: RangeBounds<usize>>(&self, range: R) -> (usize, usize) {
let l = match range.start_bound() {
Bound::Included(&x) => x,
Bound::Excluded(&x) => x.checked_add(1).expect("range bound overflow"),
Bound::Unbounded => 0,
};
let r = match range.end_bound() {
Bound::Included(&x) => x.checked_add(1).expect("range bound overflow"),
Bound::Excluded(&x) => x,
Bound::Unbounded => self.n,
};
assert!(l <= r && r <= self.n, "range out of bounds");
(l, r)
}
}
}
pub mod rng_trait {
#![warn(missing_docs)]
use std::ops::{Bound, RangeBounds};
pub trait SmallRng {
fn next_u64(&mut self) -> u64;
fn next_u8(&mut self) -> u8 {
self.next_u64() as u8
}
fn next_i8(&mut self) -> i8 {
self.next_u64() as i8
}
fn next_u16(&mut self) -> u16 {
self.next_u64() as u16
}
fn next_i16(&mut self) -> i16 {
self.next_u64() as i16
}
fn next_u32(&mut self) -> u32 {
self.next_u64() as u32
}
fn next_i32(&mut self) -> i32 {
self.next_u64() as i32
}
fn next_i64(&mut self) -> i64 {
self.next_u64() as i64
}
fn next_u128(&mut self) -> u128 {
(u128::from(self.next_u64()) << 64) | u128::from(self.next_u64())
}
fn next_i128(&mut self) -> i128 {
self.next_u128() as i128
}
fn next_usize(&mut self) -> usize {
self.next_u64() as usize
}
fn next_isize(&mut self) -> isize {
self.next_u64() as isize
}
fn next_bool(&mut self) -> bool {
self.next_u64() & 1 != 0
}
fn range_u64<R: RangeBounds<u64>>(&mut self, range: R) -> u64 {
let start = match range.start_bound() {
Bound::Included(&x) => u128::from(x),
Bound::Excluded(&x) => u128::from(x) + 1,
Bound::Unbounded => 0,
};
let end = match range.end_bound() {
Bound::Included(&x) => u128::from(x) + 1,
Bound::Excluded(&x) => u128::from(x),
Bound::Unbounded => 1_u128 << 64,
};
assert!(start < end, "random range must be non-empty and ordered");
(start + sample_below(self, end - start)) as u64
}
fn range_i64<R: RangeBounds<i64>>(&mut self, range: R) -> i64 {
const BIAS: i128 = 1_i128 << 63;
let start = match range.start_bound() {
Bound::Included(&x) => i128::from(x) + BIAS,
Bound::Excluded(&x) => i128::from(x) + BIAS + 1,
Bound::Unbounded => 0,
};
let end = match range.end_bound() {
Bound::Included(&x) => i128::from(x) + BIAS + 1,
Bound::Excluded(&x) => i128::from(x) + BIAS,
Bound::Unbounded => 1_i128 << 64,
};
assert!(start < end, "random range must be non-empty and ordered");
(start + sample_below(self, (end - start) as u128) as i128 - BIAS) as i64
}
fn range_usize<R: RangeBounds<usize>>(&mut self, range: R) -> usize {
let start = match range.start_bound() {
Bound::Included(&x) => u128::from(x as u64),
Bound::Excluded(&x) => u128::from(x as u64) + 1,
Bound::Unbounded => 0,
};
let end = match range.end_bound() {
Bound::Included(&x) => u128::from(x as u64) + 1,
Bound::Excluded(&x) => u128::from(x as u64),
Bound::Unbounded => 1_u128 << usize::BITS,
};
assert!(start < end, "random range must be non-empty and ordered");
(start + sample_below(self, end - start)) as usize
}
fn range_isize<R: RangeBounds<isize>>(&mut self, range: R) -> isize {
let bias = 1_i128 << (isize::BITS - 1);
let start = match range.start_bound() {
Bound::Included(&x) => x as i128 + bias,
Bound::Excluded(&x) => x as i128 + bias + 1,
Bound::Unbounded => 0,
};
let end = match range.end_bound() {
Bound::Included(&x) => x as i128 + bias + 1,
Bound::Excluded(&x) => x as i128 + bias,
Bound::Unbounded => 1_i128 << isize::BITS,
};
assert!(start < end, "random range must be non-empty and ordered");
(start + sample_below(self, (end - start) as u128) as i128 - bias) as isize
}
}
fn sample_below<R: SmallRng + ?Sized>(rng: &mut R, width: u128) -> u128 {
debug_assert!((1..=1_u128 << 64).contains(&width));
if width == 1_u128 << 64 {
return u128::from(rng.next_u64());
}
let width = width as u64;
let threshold = width.wrapping_neg() % width;
loop {
let value = rng.next_u64();
if value >= threshold {
return u128::from(value % width);
}
}
}
pub trait Shuffle {
fn shuffle<R: SmallRng>(&mut self, rng: &mut R);
}
impl<T> Shuffle for [T] {
fn shuffle<R: SmallRng>(&mut self, rng: &mut R) {
for i in (1..self.len()).rev() {
let j = rng.range_usize(0..=i);
self.swap(i, j);
}
}
}
}
pub mod rollinghash {
#![warn(missing_docs)]
use std::error::Error;
use std::fmt;
use std::ops::{Bound, RangeBounds};
pub const ROLLING_HASH_MODULUS: u64 = (1_u64 << 61) - 1;
#[derive(Copy, Clone, Debug, Eq, PartialEq)]
pub enum RollingHashError {
InvalidBase(u64),
RangeOutOfBounds { endpoint: usize, len: usize },
ReversedRange { start: usize, end: usize },
LengthOutOfBounds { length: usize, max_length: usize },
}
impl fmt::Display for RollingHashError {
fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
match self {
Self::InvalidBase(base) => write!(
formatter,
"rolling-hash base {base} is outside 2..{ROLLING_HASH_MODULUS}"
),
Self::RangeOutOfBounds { endpoint, len } => {
write!(
formatter,
"range endpoint {endpoint} exceeds byte length {len}"
)
}
Self::ReversedRange { start, end } => {
write!(formatter, "rolling-hash range {start}..{end} is reversed")
}
Self::LengthOutOfBounds { length, max_length } => write!(
formatter,
"requested length {length} exceeds cached length {max_length}"
),
}
}
}
impl Error for RollingHashError {}
#[derive(Clone, Debug, Eq, PartialEq)]
pub struct RollingHash {
prefix: Vec<u64>,
powers: Vec<u64>,
base: u64,
}
impl RollingHash {
pub fn new(input: impl AsRef<[u8]>, base: u64) -> Result<Self, RollingHashError> {
if !(2..ROLLING_HASH_MODULUS).contains(&base) {
return Err(RollingHashError::InvalidBase(base));
}
let input = input.as_ref();
let mut prefix = Vec::with_capacity(input.len() + 1);
let mut powers = Vec::with_capacity(input.len() + 1);
prefix.push(0);
powers.push(1);
for &byte in input {
prefix.push(add_mod(
mul_mod(*prefix.last().unwrap(), base),
u64::from(byte),
));
powers.push(mul_mod(*powers.last().unwrap(), base));
}
Ok(Self {
prefix,
powers,
base,
})
}
#[must_use]
pub fn len(&self) -> usize {
self.prefix.len() - 1
}
#[must_use]
pub fn is_empty(&self) -> bool {
self.len() == 0
}
#[must_use]
pub const fn base(&self) -> u64 {
self.base
}
pub fn hash<R: RangeBounds<usize>>(&self, range: R) -> Result<u64, RollingHashError> {
let (start, end) = self.normalize_range(range)?;
Ok(sub_mod(
self.prefix[end],
mul_mod(self.prefix[start], self.powers[end - start]),
))
}
pub fn concatenate(
&self,
left_hash: u64,
right_hash: u64,
right_len: usize,
) -> Result<u64, RollingHashError> {
let Some(&power) = self.powers.get(right_len) else {
return Err(RollingHashError::LengthOutOfBounds {
length: right_len,
max_length: self.len(),
});
};
Ok(add_mod(mul_mod(left_hash, power), right_hash))
}
fn normalize_range<R: RangeBounds<usize>>(
&self,
range: R,
) -> Result<(usize, usize), RollingHashError> {
let start = match range.start_bound() {
Bound::Included(&value) => value,
Bound::Excluded(&value) => {
value
.checked_add(1)
.ok_or(RollingHashError::RangeOutOfBounds {
endpoint: value,
len: self.len(),
})?
}
Bound::Unbounded => 0,
};
let end = match range.end_bound() {
Bound::Included(&value) => {
value
.checked_add(1)
.ok_or(RollingHashError::RangeOutOfBounds {
endpoint: value,
len: self.len(),
})?
}
Bound::Excluded(&value) => value,
Bound::Unbounded => self.len(),
};
if start > end {
return Err(RollingHashError::ReversedRange { start, end });
}
if end > self.len() {
return Err(RollingHashError::RangeOutOfBounds {
endpoint: end,
len: self.len(),
});
}
Ok((start, end))
}
}
fn mul_mod(left: u64, right: u64) -> u64 {
reduce(u128::from(left) * u128::from(right))
}
fn add_mod(left: u64, right: u64) -> u64 {
let sum = left + right;
if sum >= ROLLING_HASH_MODULUS {
sum - ROLLING_HASH_MODULUS
} else {
sum
}
}
fn sub_mod(left: u64, right: u64) -> u64 {
if left >= right {
left - right
} else {
left + ROLLING_HASH_MODULUS - right
}
}
fn reduce(product: u128) -> u64 {
let modulus = u128::from(ROLLING_HASH_MODULUS);
let folded = (product >> 61) + (product & modulus);
let folded = (folded >> 61) + (folded & modulus);
let value = folded as u64;
if value >= ROLLING_HASH_MODULUS {
value - ROLLING_HASH_MODULUS
} else {
value
}
}
}
pub mod scc {
use crate::internal_scc;
#[derive(Clone, Debug)]
pub struct SccGraph {
internal: internal_scc::SccGraph,
}
impl SccGraph {
pub fn new(n: usize) -> Self {
SccGraph {
internal: internal_scc::SccGraph::new(n),
}
}
pub fn add_edge(&mut self, from: usize, to: usize) {
let n = self.internal.num_vertices();
assert!(from < n);
assert!(to < n);
self.internal.add_edge(from, to);
}
pub fn scc(&self) -> Vec<Vec<usize>> {
self.internal.scc()
}
}
}
pub mod segtree {
use crate::internal_bit::ceil_pow2;
use crate::internal_type_traits::{BoundedAbove, BoundedBelow, One, Zero};
use std::cmp::{max, min};
use std::convert::Infallible;
use std::iter::FromIterator;
use std::marker::PhantomData;
use std::ops::{Add, BitAnd, BitOr, BitXor, Bound, Mul, Not, RangeBounds};
pub trait Monoid {
type S: Clone;
fn identity() -> Self::S;
fn binary_operation(a: &Self::S, b: &Self::S) -> Self::S;
}
pub struct Max<S>(Infallible, PhantomData<fn() -> S>);
impl<S> Monoid for Max<S>
where
S: Copy + Ord + BoundedBelow,
{
type S = S;
fn identity() -> Self::S {
S::min_value()
}
fn binary_operation(a: &Self::S, b: &Self::S) -> Self::S {
max(*a, *b)
}
}
pub struct Min<S>(Infallible, PhantomData<fn() -> S>);
impl<S> Monoid for Min<S>
where
S: Copy + Ord + BoundedAbove,
{
type S = S;
fn identity() -> Self::S {
S::max_value()
}
fn binary_operation(a: &Self::S, b: &Self::S) -> Self::S {
min(*a, *b)
}
}
pub struct Additive<S>(Infallible, PhantomData<fn() -> S>);
impl<S> Monoid for Additive<S>
where
S: Copy + Add<Output = S> + Zero,
{
type S = S;
fn identity() -> Self::S {
S::zero()
}
fn binary_operation(a: &Self::S, b: &Self::S) -> Self::S {
*a + *b
}
}
pub struct Multiplicative<S>(Infallible, PhantomData<fn() -> S>);
impl<S> Monoid for Multiplicative<S>
where
S: Copy + Mul<Output = S> + One,
{
type S = S;
fn identity() -> Self::S {
S::one()
}
fn binary_operation(a: &Self::S, b: &Self::S) -> Self::S {
*a * *b
}
}
pub struct BitwiseOr<S>(Infallible, PhantomData<fn() -> S>);
impl<S> Monoid for BitwiseOr<S>
where
S: Copy + BitOr<Output = S> + Zero,
{
type S = S;
fn identity() -> Self::S {
S::zero()
}
fn binary_operation(a: &Self::S, b: &Self::S) -> Self::S {
*a | *b
}
}
pub struct BitwiseAnd<S>(Infallible, PhantomData<fn() -> S>);
impl<S> Monoid for BitwiseAnd<S>
where
S: Copy + BitAnd<Output = S> + Not<Output = S> + Zero,
{
type S = S;
fn identity() -> Self::S {
!S::zero()
}
fn binary_operation(a: &Self::S, b: &Self::S) -> Self::S {
*a & *b
}
}
pub struct BitwiseXor<S>(Infallible, PhantomData<fn() -> S>);
impl<S> Monoid for BitwiseXor<S>
where
S: Copy + BitXor<Output = S> + Zero,
{
type S = S;
fn identity() -> Self::S {
S::zero()
}
fn binary_operation(a: &Self::S, b: &Self::S) -> Self::S {
*a ^ *b
}
}
impl<M: Monoid> Default for Segtree<M> {
fn default() -> Self {
Segtree::new(0)
}
}
impl<M: Monoid> Segtree<M> {
pub fn new(n: usize) -> Segtree<M> {
vec![M::identity(); n].into()
}
}
impl<M: Monoid> From<Vec<M::S>> for Segtree<M> {
fn from(v: Vec<M::S>) -> Self {
let n = v.len();
let log = ceil_pow2(n as u32) as usize;
let size = 1 << log;
let mut d = vec![M::identity(); 2 * size];
d[size..][..n].clone_from_slice(&v);
let mut ret = Segtree { n, size, log, d };
for i in (1..size).rev() {
ret.update(i);
}
ret
}
}
impl<M: Monoid> FromIterator<M::S> for Segtree<M> {
fn from_iter<T: IntoIterator<Item = M::S>>(iter: T) -> Self {
let v = iter.into_iter().collect::<Vec<_>>();
v.into()
}
}
impl<M: Monoid> Segtree<M> {
pub fn set(&mut self, mut p: usize, x: M::S) {
assert!(p < self.n);
p += self.size;
self.d[p] = x;
for i in 1..=self.log {
self.update(p >> i);
}
}
pub fn get(&self, p: usize) -> M::S {
assert!(p < self.n);
self.d[p + self.size].clone()
}
pub fn get_slice(&self) -> &[M::S] {
&self.d[self.size..][..self.n]
}
pub fn prod<R>(&self, range: R) -> M::S
where
R: RangeBounds<usize>,
{
if range.start_bound() == Bound::Unbounded && range.end_bound() == Bound::Unbounded {
return self.all_prod();
}
let mut r = match range.end_bound() {
Bound::Included(r) => r + 1,
Bound::Excluded(r) => *r,
Bound::Unbounded => self.n,
};
let mut l = match range.start_bound() {
Bound::Included(l) => *l,
Bound::Excluded(l) => l + 1,
Bound::Unbounded => 0,
};
assert!(l <= r && r <= self.n);
let mut sml = M::identity();
let mut smr = M::identity();
l += self.size;
r += self.size;
while l < r {
if l & 1 != 0 {
sml = M::binary_operation(&sml, &self.d[l]);
l += 1;
}
if r & 1 != 0 {
r -= 1;
smr = M::binary_operation(&self.d[r], &smr);
}
l >>= 1;
r >>= 1;
}
M::binary_operation(&sml, &smr)
}
pub fn all_prod(&self) -> M::S {
self.d[1].clone()
}
pub fn max_right<F>(&self, mut l: usize, f: F) -> usize
where
F: Fn(&M::S) -> bool,
{
assert!(l <= self.n);
assert!(f(&M::identity()));
if l == self.n {
return self.n;
}
l += self.size;
let mut sm = M::identity();
while {
while l.is_multiple_of(2) {
l >>= 1;
}
if !f(&M::binary_operation(&sm, &self.d[l])) {
while l < self.size {
l *= 2;
let res = M::binary_operation(&sm, &self.d[l]);
if f(&res) {
sm = res;
l += 1;
}
}
return l - self.size;
}
sm = M::binary_operation(&sm, &self.d[l]);
l += 1;
{
let l = l as isize;
(l & -l) != l
}
} {}
self.n
}
pub fn min_left<F>(&self, mut r: usize, f: F) -> usize
where
F: Fn(&M::S) -> bool,
{
assert!(r <= self.n);
assert!(f(&M::identity()));
if r == 0 {
return 0;
}
r += self.size;
let mut sm = M::identity();
while {
r -= 1;
while r > 1 && r % 2 == 1 {
r >>= 1;
}
if !f(&M::binary_operation(&self.d[r], &sm)) {
while r < self.size {
r = 2 * r + 1;
let res = M::binary_operation(&self.d[r], &sm);
if f(&res) {
sm = res;
r -= 1;
}
}
return r + 1 - self.size;
}
sm = M::binary_operation(&self.d[r], &sm);
{
let r = r as isize;
(r & -r) != r
}
} {}
0
}
fn update(&mut self, k: usize) {
self.d[k] = M::binary_operation(&self.d[2 * k], &self.d[2 * k + 1]);
}
}
#[derive(Clone)]
pub struct Segtree<M>
where
M: Monoid,
{
n: usize,
size: usize,
log: usize,
d: Vec<M::S>,
}
}
pub mod segtree_template {
use super::lazysegtree::*;
use super::segtree::*;
pub struct MaxAdd;
impl MapMonoid for MaxAdd {
type M = Max<i64>;
type F = i64;
fn identity_map() -> Self::F {
0
}
fn mapping(&f: &i64, &x: &i64) -> i64 {
f + x
}
fn composition(&f: &i64, &g: &i64) -> i64 {
f + g
}
}
pub struct MinAdd;
impl MapMonoid for MinAdd {
type M = Max<i64>;
type F = i64;
fn identity_map() -> Self::F {
0
}
fn mapping(&f: &i64, &x: &i64) -> i64 {
f + x
}
fn composition(&f: &i64, &g: &i64) -> i64 {
f + g
}
}
const INF: i64 = i64::MAX;
pub struct MaxForLazy;
impl Monoid for MaxForLazy {
type S = i64;
fn identity() -> Self::S {
INF
}
fn binary_operation(a: &Self::S, b: &Self::S) -> Self::S {
if *a == INF {
*b
} else if *b == INF {
*a
} else {
std::cmp::max(*a, *b)
}
}
}
pub struct MinForLazy;
impl Monoid for MinForLazy {
type S = i64;
fn identity() -> Self::S {
INF
}
fn binary_operation(a: &Self::S, b: &Self::S) -> Self::S {
if *a == INF {
*b
} else if *b == INF {
*a
} else {
std::cmp::min(*a, *b)
}
}
}
pub struct MaxUpdate;
impl MapMonoid for MaxUpdate {
type M = MaxForLazy;
type F = i64;
fn identity_map() -> Self::F {
INF
}
fn mapping(&f: &i64, &x: &i64) -> i64 {
if f == INF {
x
} else {
f
}
}
fn composition(&f: &i64, &g: &i64) -> i64 {
if f == INF {
g
} else {
f
}
}
}
pub struct MinUpdate;
impl MapMonoid for MinUpdate {
type M = MinForLazy;
type F = i64;
fn identity_map() -> Self::F {
INF
}
fn mapping(&f: &i64, &x: &i64) -> i64 {
if f == INF {
x
} else {
f
}
}
fn composition(&f: &i64, &g: &i64) -> i64 {
if f == INF {
g
} else {
f
}
}
}
}
pub mod string {
#![allow(
clippy::many_single_char_names,
clippy::implicit_saturating_sub,
clippy::bool_to_int_with_if
)]
fn sa_naive<T: Ord>(s: &[T]) -> Vec<usize> {
let n = s.len();
let mut sa: Vec<usize> = (0..n).collect();
sa.sort_by(|&(mut l), &(mut r)| {
if l == r {
return std::cmp::Ordering::Equal;
}
while l < n && r < n {
if s[l] != s[r] {
return s[l].cmp(&s[r]);
}
l += 1;
r += 1;
}
if l == n {
std::cmp::Ordering::Less
} else {
std::cmp::Ordering::Greater
}
});
sa
}
fn sa_doubling(s: &[i32]) -> Vec<usize> {
let n = s.len();
let mut sa: Vec<usize> = (0..n).collect();
let mut rnk: Vec<i32> = s.to_vec();
let mut tmp = vec![0; n];
let mut k = 1;
while k < n {
let cmp = |&x: &usize, &y: &usize| {
if rnk[x] != rnk[y] {
return rnk[x].cmp(&rnk[y]);
}
let rx = if x + k < n { rnk[x + k] } else { -1 };
let ry = if y + k < n { rnk[y + k] } else { -1 };
rx.cmp(&ry)
};
sa.sort_by(cmp);
tmp[sa[0]] = 0;
for i in 1..n {
tmp[sa[i]] =
tmp[sa[i - 1]] + i32::from(cmp(&sa[i - 1], &sa[i]) == std::cmp::Ordering::Less);
}
std::mem::swap(&mut tmp, &mut rnk);
k *= 2;
}
sa
}
trait Threshold {
fn threshold_naive() -> usize;
fn threshold_doubling() -> usize;
}
enum DefaultThreshold {}
impl Threshold for DefaultThreshold {
fn threshold_naive() -> usize {
10
}
fn threshold_doubling() -> usize {
40
}
}
#[allow(clippy::cognitive_complexity)]
fn sa_is<T: Threshold>(s: &[usize], upper: usize) -> Vec<usize> {
let n = s.len();
match n {
0 => return vec![],
1 => return vec![0],
2 => return if s[0] < s[1] { vec![0, 1] } else { vec![1, 0] },
_ => (),
}
if n < T::threshold_naive() {
return sa_naive(s);
}
if n < T::threshold_doubling() {
let s: Vec<i32> = s.iter().map(|&x| x as i32).collect();
return sa_doubling(&s);
}
let mut sa = vec![0; n];
let mut ls = vec![false; n];
for i in (0..n - 1).rev() {
ls[i] = if s[i] == s[i + 1] {
ls[i + 1]
} else {
s[i] < s[i + 1]
};
}
let mut sum_l = vec![0; upper + 1];
let mut sum_s = vec![0; upper + 1];
for i in 0..n {
if !ls[i] {
sum_s[s[i]] += 1;
} else {
sum_l[s[i] + 1] += 1;
}
}
for i in 0..=upper {
sum_s[i] += sum_l[i];
if i < upper {
sum_l[i + 1] += sum_s[i];
}
}
let induce = |sa: &mut [usize], lms: &[usize]| {
for elem in sa.iter_mut() {
*elem = 0;
}
let mut buf = sum_s.clone();
for &d in lms {
if d == n {
continue;
}
let old = buf[s[d]];
buf[s[d]] += 1;
sa[old] = d + 1;
}
buf.copy_from_slice(&sum_l);
let old = buf[s[n - 1]];
buf[s[n - 1]] += 1;
sa[old] = n;
for i in 0..n {
let v = sa[i];
if v >= 2 && !ls[v - 2] {
let old = buf[s[v - 2]];
buf[s[v - 2]] += 1;
sa[old] = v - 1;
}
}
buf.copy_from_slice(&sum_l);
for i in (0..n).rev() {
let v = sa[i];
if v >= 2 && ls[v - 2] {
buf[s[v - 2] + 1] -= 1;
sa[buf[s[v - 2] + 1]] = v - 1;
}
}
};
let mut lms_map = vec![0; n + 1];
let mut m = 0;
for i in 1..n {
if !ls[i - 1] && ls[i] {
lms_map[i] = m + 1;
m += 1;
}
}
let mut lms = Vec::with_capacity(m);
for i in 1..n {
if !ls[i - 1] && ls[i] {
lms.push(i);
}
}
assert_eq!(lms.len(), m);
induce(&mut sa, &lms);
if m > 0 {
let mut sorted_lms = Vec::with_capacity(m);
for &v in &sa {
if lms_map[v - 1] != 0 {
sorted_lms.push(v - 1);
}
}
let mut rec_s = vec![0; m];
let mut rec_upper = 0;
rec_s[lms_map[sorted_lms[0]] - 1] = 0;
for i in 1..m {
let mut l = sorted_lms[i - 1];
let mut r = sorted_lms[i];
let end_l = if lms_map[l] < m { lms[lms_map[l]] } else { n };
let end_r = if lms_map[r] < m { lms[lms_map[r]] } else { n };
let same = if end_l - l != end_r - r {
false
} else {
while l < end_l {
if s[l] != s[r] {
break;
}
l += 1;
r += 1;
}
l != n && s[l] == s[r]
};
if !same {
rec_upper += 1;
}
rec_s[lms_map[sorted_lms[i]] - 1] = rec_upper;
}
let rec_sa = sa_is::<T>(&rec_s, rec_upper);
for i in 0..m {
sorted_lms[i] = lms[rec_sa[i]];
}
induce(&mut sa, &mut sorted_lms);
}
for elem in sa.iter_mut() {
*elem -= 1;
}
sa
}
fn sa_is_i32<T: Threshold>(s: &[i32], upper: i32) -> Vec<usize> {
let s: Vec<usize> = s.iter().map(|&x| x as usize).collect();
sa_is::<T>(&s, upper as usize)
}
pub fn suffix_array_manual(s: &[i32], upper: i32) -> Vec<usize> {
assert!(upper >= 0);
for &elem in s {
assert!(0 <= elem && elem <= upper);
}
sa_is_i32::<DefaultThreshold>(s, upper)
}
pub fn suffix_array_arbitrary<T: Ord>(s: &[T]) -> Vec<usize> {
let n = s.len();
let mut idx: Vec<usize> = (0..n).collect();
idx.sort_by_key(|&i| &s[i]);
let mut s2 = vec![0; n];
let mut now = 0;
for i in 0..n {
if i > 0 && s[idx[i - 1]] != s[idx[i]] {
now += 1;
}
s2[idx[i]] = now;
}
sa_is_i32::<DefaultThreshold>(&s2, now)
}
pub fn suffix_array(s: &str) -> Vec<usize> {
let s2: Vec<usize> = s.bytes().map(|x| x as usize).collect();
sa_is::<DefaultThreshold>(&s2, 255)
}
pub fn lcp_array_arbitrary<T: Ord>(s: &[T], sa: &[usize]) -> Vec<usize> {
let n = s.len();
assert!(n >= 1);
let mut rnk = vec![0; n];
for i in 0..n {
rnk[sa[i]] = i;
}
let mut lcp = vec![0; n - 1];
let mut h: usize = 0;
for i in 0..n - 1 {
h = h.saturating_sub(1);
if rnk[i] == 0 {
continue;
}
let j = sa[rnk[i] - 1];
while j + h < n && i + h < n {
if s[j + h] != s[i + h] {
break;
}
h += 1;
}
lcp[rnk[i] - 1] = h;
}
lcp
}
pub fn lcp_array(s: &str, sa: &[usize]) -> Vec<usize> {
let s: &[u8] = s.as_bytes();
lcp_array_arbitrary(s, sa)
}
pub fn z_algorithm_arbitrary<T: Ord>(s: &[T]) -> Vec<usize> {
let n = s.len();
if n == 0 {
return vec![];
}
let mut z = vec![0; n];
z[0] = 0;
let mut j = 0;
for i in 1..n {
let mut k = if j + z[j] <= i {
0
} else {
std::cmp::min(j + z[j] - i, z[i - j])
};
while i + k < n && s[k] == s[i + k] {
k += 1;
}
z[i] = k;
if j + z[j] < i + z[i] {
j = i;
}
}
z[0] = n;
z
}
pub fn z_algorithm(s: &str) -> Vec<usize> {
let s: &[u8] = s.as_bytes();
z_algorithm_arbitrary(s)
}
}
pub mod timer {
#![warn(missing_docs)]
use std::time::{Duration, Instant};
#[derive(Clone, Debug)]
pub struct Timer {
limit: Duration,
started: Instant,
}
impl Timer {
#[must_use]
pub fn new(limit: Duration) -> Self {
Self {
limit,
started: Instant::now(),
}
}
#[must_use]
pub fn from_millis(milliseconds: u64) -> Self {
Self::new(Duration::from_millis(milliseconds))
}
#[must_use]
pub fn elapsed(&self) -> Duration {
self.started.elapsed()
}
#[must_use]
pub fn remaining(&self) -> Duration {
self.limit.saturating_sub(self.elapsed())
}
#[must_use]
pub fn is_within_limit(&self) -> bool {
self.elapsed() < self.limit
}
pub fn reset(&mut self) {
self.started = Instant::now();
}
#[must_use]
pub fn formatted_elapsed(&self) -> String {
format_duration(self.elapsed())
}
}
#[must_use]
pub fn format_duration(duration: Duration) -> String {
let nanos = duration.as_nanos();
if nanos < 1_000 {
format!("{nanos} ns")
} else if nanos < 1_000_000 {
format!("{} μs", duration.as_micros())
} else if nanos < 1_000_000_000 {
format!("{} ms", duration.as_millis())
} else {
format!("{} s", duration.as_secs())
}
}
}
pub mod twosat {
use crate::internal_scc;
#[derive(Clone, Debug)]
pub struct TwoSat {
n: usize,
scc: internal_scc::SccGraph,
answer: Vec<bool>,
}
impl TwoSat {
pub fn new(n: usize) -> Self {
TwoSat {
n,
answer: vec![false; n],
scc: internal_scc::SccGraph::new(2 * n),
}
}
pub fn add_clause(&mut self, i: usize, f: bool, j: usize, g: bool) {
assert!(i < self.n && j < self.n);
self.scc.add_edge(2 * i + !f as usize, 2 * j + g as usize);
self.scc.add_edge(2 * j + !g as usize, 2 * i + f as usize);
}
pub fn satisfiable(&mut self) -> bool {
let id = self.scc.scc_ids().1;
for i in 0..self.n {
if id[2 * i] == id[2 * i + 1] {
return false;
}
self.answer[i] = id[2 * i] < id[2 * i + 1];
}
true
}
pub fn answer(&self) -> &[bool] {
&self.answer
}
}
}
pub mod vectools {
#![warn(missing_docs)]
use std::cmp::Ordering;
pub trait BSearch<T> {
fn lower_bound(&self, key: &T) -> usize;
fn upper_bound(&self, key: &T) -> usize;
fn equal_range(&self, key: &T) -> std::ops::Range<usize>;
}
impl<T: Ord> BSearch<T> for [T] {
fn lower_bound(&self, key: &T) -> usize {
self.partition_point(|x| x < key)
}
fn upper_bound(&self, key: &T) -> usize {
self.partition_point(|x| x <= key)
}
fn equal_range(&self, key: &T) -> std::ops::Range<usize> {
self.lower_bound(key)..self.upper_bound(key)
}
}
pub fn next_permutation<T: Ord>(values: &mut [T]) -> bool {
next_permutation_by(values, Ord::cmp)
}
pub fn next_permutation_by<T, F>(values: &mut [T], mut compare: F) -> bool
where
F: FnMut(&T, &T) -> Ordering,
{
let Some(pivot) = (1..values.len()).rfind(|&i| compare(&values[i - 1], &values[i]).is_lt())
else {
return false;
};
let pivot = pivot - 1;
let swap = (pivot + 1..values.len())
.rfind(|&i| compare(&values[pivot], &values[i]).is_lt())
.unwrap();
values.swap(pivot, swap);
values[pivot + 1..].reverse();
true
}
pub fn prev_permutation<T: Ord>(values: &mut [T]) -> bool {
prev_permutation_by(values, Ord::cmp)
}
pub fn prev_permutation_by<T, F>(values: &mut [T], mut compare: F) -> bool
where
F: FnMut(&T, &T) -> Ordering,
{
let Some(pivot) = (1..values.len()).rfind(|&i| compare(&values[i - 1], &values[i]).is_gt())
else {
return false;
};
let pivot = pivot - 1;
let swap = (pivot + 1..values.len())
.rfind(|&i| compare(&values[pivot], &values[i]).is_gt())
.unwrap();
values.swap(pivot, swap);
values[pivot + 1..].reverse();
true
}
#[must_use]
pub fn coordinate_compress<T: Ord + Clone>(values: &[T]) -> (Vec<usize>, Vec<T>) {
let mut unique = values.to_vec();
unique.sort_unstable();
unique.dedup();
let ranks = values
.iter()
.map(|x| unique.binary_search(x).unwrap())
.collect();
(ranks, unique)
}
#[must_use]
pub fn uppercase_to_usize(c: char) -> Option<usize> {
c.is_ascii_uppercase().then(|| (c as u8 - b'A') as usize)
}
#[must_use]
pub fn lowercase_to_usize(c: char) -> Option<usize> {
c.is_ascii_lowercase().then(|| (c as u8 - b'a') as usize)
}
#[must_use]
pub fn usize_to_uppercase(index: usize) -> Option<char> {
(index < 26).then(|| (b'A' + index as u8) as char)
}
#[must_use]
pub fn usize_to_lowercase(index: usize) -> Option<char> {
(index < 26).then(|| (b'a' + index as u8) as char)
}
}
pub mod xorshift {
#![warn(missing_docs)]
use crate::SmallRng;
use std::time::{SystemTime, SystemTimeError, UNIX_EPOCH};
const DEFAULT_SEED: u64 = 88_172_645_463_325_252;
const NON_ZERO_FALLBACK: u64 = 0x9e37_79b9_7f4a_7c15;
#[derive(Clone, Debug)]
pub struct XorShift {
initial_seed: u64,
state: u64,
}
impl Default for XorShift {
fn default() -> Self {
Self::seed_from_u64(DEFAULT_SEED)
}
}
impl XorShift {
#[must_use]
pub fn seed_from_u64(seed: u64) -> Self {
Self {
initial_seed: seed,
state: if seed == 0 { NON_ZERO_FALLBACK } else { seed },
}
}
pub fn from_system_time() -> Result<Self, SystemTimeError> {
let duration = SystemTime::now().duration_since(UNIX_EPOCH)?;
Ok(Self::seed_from_u64(
duration.as_secs() ^ u64::from(duration.subsec_nanos()).rotate_left(32),
))
}
#[must_use]
pub fn seed(&self) -> u64 {
self.initial_seed
}
}
impl SmallRng for XorShift {
fn next_u64(&mut self) -> u64 {
let mut x = self.state;
x ^= x >> 12;
x ^= x << 25;
x ^= x >> 27;
self.state = x;
x.wrapping_mul(0x2545_f491_4f6c_dd1d)
}
}
}