// codesnip-guard: PathWeightedMatching #[cfg_attr(any(),rustfmt::skip)]pub use self::path_weighted_matching::{PathWeightedMatching,path_weighted_matching_queries};#[cfg_attr(any(),rustfmt::skip)]mod path_weighted_matching{use super::{Associative,Magma,RangeMinimumQuery,SegmentTree,Unital,WaveletMatrix,Zero};use std::{cmp::Reverse,ops::{Add,Sub}};#[derive(Debug,Clone)]pub struct PathWeightedMatching{weights:Vec}implPathWeightedMatching{pub fn from_edge_weights(weights:impl IntoIterator)->Self where T:Copy+Zero+Sub{let mut vertex=T::zero();let vertices=std::iter::once(vertex).chain(weights.into_iter().map(|weight|{vertex=weight-vertex;vertex}));Self::from_vertex_weights(vertices)}#[doc=" Edge weights are the sums of adjacent vertex weights."]pub fn from_vertex_weights(weights:impl IntoIterator)->Self{Self{weights:weights.into_iter().collect()}}#[doc=" Returns the maximum weight of a matching with exactly `k` edges in O(n) time."]pub fn maximum_weight(self,k:usize)->T where T:Copy+Ord+Zero+Add{assert!(k<=self.weights.len()/2);let mut gains=path_matching_candidates(self.weights,|x|x,|a,b|a+b);if k!=0&&kVecwhere T:Copy+Ord+Zero+Add{let mut gains=path_matching_candidates(self.weights,|x|x,|a,b|a+b);gains.sort_unstable_by_key(|x|Reverse(*x));gains.insert(0,T::zero());for i in 1..gains.len(){gains[i]=gains[i-1]+gains[i];}gains}#[doc=" Returns the maximum weight and sorted edge indices for exactly `k` edges in O(n) time."]pub fn maximum_weight_matching(self,k:usize)->(T,Vec)where T:Copy+Ord+Zero+Add{let n=self.weights.len();assert!(k<=n/2);let mut pairs=path_matching_candidates((0..n).collect(),|i|self.weights[i],|u,v|(self.weights[u]+self.weights[v],u,v));if k!=0&&k(weights:&[T],queries:&[(usize,usize,usize)])->Vecwhere T:Copy+Ord+Zero+Add+Sub{let n=weights.len();let mut prefix=Vec::with_capacity(n+1);prefix.push(T::zero());for&weight in weights{prefix.push(*prefix.last().unwrap()+weight);}let mut answers=vec![T::zero();queries.len()];let mut pending=false;for(i,&(l,r,k))in queries.iter().enumerate(){assert!(l<=r&&r<=n&&k<=(r-l)/2);if 2*k==r-l{answers[i]=prefix[r]-prefix[l];}else if k!=0{pending=true;}}if!pending{return answers;}let mut pairs=path_matching_candidates((0..n).collect(),|i|weights[i],|u,v|(weights[u]+weights[v],u,v));pairs.sort_unstable_by_key(|&(w,u,v)|(Reverse(w),u,v));let mut times=vec![pairs.len();n];for(t,&(_,u,v))in pairs.iter().enumerate(){times[u]=t;times[v]=t;}let times=WaveletMatrix::new(times);let mut heads=vec![!0;pairs.len()];let mut next=vec![!0;queries.len()];for(i,&(l,r,k))in queries.iter().enumerate(){if k!=0&&2*k!=r-l{let t=times.quantile(l..r,2*k-1);next[i]=heads[t];heads[t]=i;}}let minima:[_;2]=std::array::from_fn(|parity|{RangeMinimumQuery::new((parity..n).step_by(2).map(|i|(weights[i],Reverse(i))).collect())});let mut state=SegmentTree::>::from_vec(weights.iter().enumerate().map(|(i,&w)|RangeState::unmatched(w,i)).collect());let mut rollback=Vec::with_capacity(6);for(t,&(_,u,v))in pairs.iter().enumerate(){let mut query=heads[t];while query!=!0{let(l,r,k)=queries[query];let current=state.fold(l..r);let used=r-l-current.free;let left=(l-state.fold(..l).free)%2;let right=(left+used)%2;let mut answer=current.sum;if left!=0{let first=current.first.unwrap().1;let(_,Reverse(i))=minima[l%2].fold(l/2,first.div_ceil(2));rollback.push((i,state.get(i)));state.set(i,RangeState::unmatched(weights[i],i));answer=answer-weights[i];}if right!=0{let last=current.last.unwrap().1;let(_,Reverse(i))=minima[(r-1)%2].fold(last.div_ceil(2),r.div_ceil(2));rollback.push((i,state.get(i)));state.set(i,RangeState::unmatched(weights[i],i));answer=answer-weights[i];}for _ in 0..(2*k-used+left+right)/2{let(gain,Reverse(a),Reverse(b))=state.fold(l..r).best.unwrap();answer=answer+gain;for i in[a,b]{rollback.push((i,state.get(i)));state.set(i,RangeState::matched(weights[i]));}}answers[query]=answer;while let Some((i,value))=rollback.pop(){state.set(i,value);}query=next[query];}state.set(u,RangeState::matched(weights[u]));state.set(v,RangeState::matched(weights[v]));}answers}fn path_matching_candidates(mut stack:Vec,weight:impl Fn(V)->T,pair:impl Fn(V,V)->R)->Vecwhere V:Copy,T:Ord{let mut pairs=Vec::with_capacity(stack.len()/2+1);let mut n=0;for i in 0..stack.len(){stack[n]=stack[i];n+=1;while n>=3{let local_maximum=weight(stack[n-3])>=weight(stack[n-1])&&(n==3||weight(stack[n-2])>weight(stack[n-4]));if!local_maximum{break;}pairs.push(pair(stack[n-3],stack[n-2]));stack[n-3]=stack[n-1];n-=2;}}while n>=2{pairs.push(pair(stack[n-2],stack[n-1]));n-=2;}pairs}#[derive(Clone)]struct RangeState{free:usize,sum:T,first:Option<(T,usize)>,last:Option<(T,usize)>,best:Option<(T,Reverse,Reverse)>}implRangeState{fn unmatched(weight:T,index:usize)->Self where T:Copy+Zero{Self{free:1,sum:T::zero(),first:Some((weight,index)),last:Some((weight,index)),best:None}}fn matched(weight:T)->Self{Self{free:0,sum:weight,first:None,last:None,best:None}}}implMagma for RangeStatewhere T:Copy+Ord+Add{type T=Self;fn operate(x:&Self,y:&Self)->Self{let crossing=x.last.zip(y.first).map(|((a,i),(b,j))|(a+b,Reverse(i),Reverse(j)));Self{free:x.free+y.free,sum:x.sum+y.sum,first:x.first.or(y.first),last:y.last.or(x.last),best:x.best.max(y.best).max(crossing)}}}implUnital for RangeStatewhere T:Copy+Ord+Zero+Add{fn unit()->Self{Self::matched(T::zero())}}implAssociative for RangeStatewhere T:Copy+Ord+Add{}} // codesnip-guard: BinaryIndexedTree #[cfg_attr(any(),rustfmt::skip)]pub use self::binary_indexed_tree::BinaryIndexedTree;#[cfg_attr(any(),rustfmt::skip)]mod binary_indexed_tree{use super::{AbelianGroup,Group,Monoid};use std::fmt::{self,Debug,Formatter};pub struct BinaryIndexedTreewhere M:Monoid{n:usize,bit:Vec}implClone for BinaryIndexedTreewhere M:Monoid{fn clone(&self)->Self{Self{n:self.n,bit:self.bit.clone()}}}implDebug for BinaryIndexedTreewhere M:Monoid{fn fmt(&self,f:&mut Formatter<'_>)->fmt::Result{f.debug_struct("BinaryIndexedTree").field("n",&self.n).field("bit",&self.bit).finish()}}implBinaryIndexedTreewhere M:Monoid{#[inline]pub fn new(n:usize)->Self{let bit=vec![M::unit();n+1];Self{n,bit}}#[inline]pub fn from_slice(slice:&[M::T])->Self{let n=slice.len();let mut bit=vec![M::unit();n+1];for(i,x)in slice.iter().enumerate(){let k=i+1;M::operate_assign(&mut bit[k],x);let j=k+(k&(!k+1));if j<=n{bit[j]=M::operate(&bit[j],&bit[k]);}}Self{n,bit}}#[inline]#[doc=" fold [0, k)"]pub fn accumulate0(&self,mut k:usize)->M::T{debug_assert!(k<=self.n);let mut res=M::unit();while k>0{res=M::operate(&self.bit[k],&res);k-=k&(!k+1);}res}#[inline]#[doc=" fold [0, k]"]pub fn accumulate(&self,k:usize)->M::T{self.accumulate0(k+1)}#[inline]pub fn update(&mut self,k:usize,x:M::T){debug_assert!(k(&self,mut pred:P)->usize where P:FnMut(&M::T)->bool{let n=self.n;let mut acc=M::unit();let mut pos=0;let mut k=n.next_power_of_two();while k>0{if k+pos<=n{let nacc=M::operate(&acc,&self.bit[k+pos]);if pred(&nacc){pos+=k;acc=nacc;}}k>>=1;}pos}}implBinaryIndexedTree{#[inline]pub fn fold(&self,l:usize,r:usize)->G::T{debug_assert!(l<=self.n&&r<=self.n);G::operate(&G::inverse(&self.accumulate0(l)),&self.accumulate0(r))}#[inline]pub fn fold_abelian(&self,mut l:usize,mut r:usize)->G::T where G:AbelianGroup{debug_assert!(l<=self.n&&r<=self.n);if l==r{return G::unit();}let common=l&!(usize::MAX>>(l^r).leading_zeros());let mut left=G::unit();let mut right=G::unit();while l!=common{G::operate_assign(&mut left,&self.bit[l]);l&=l-1;}while r!=common{G::operate_assign(&mut right,&self.bit[r]);r&=r-1;}G::rinv_operate(&right,&left)}#[inline]pub fn get(&self,k:usize)->G::T{self.fold(k,k+1)}#[inline]pub fn set(&mut self,k:usize,x:G::T){self.update(k,G::operate(&G::inverse(&self.get(k)),&x));}}} // codesnip-guard: BitVector #[cfg_attr(any(),rustfmt::skip)]pub use self::bit_vector::{BitVector,BitVectorBlock,RankSelectDictionaries};#[cfg_attr(any(),rustfmt::skip)]mod bit_vector{use std::iter::FromIterator;#[doc=" rank_i(select_i(k)) = k"]#[doc=" rank_i(select_i(k) + 1) = k + 1"]pub trait RankSelectDictionaries{fn bit_length(&self)->usize;#[doc=" get k-th bit"]fn access(&self,k:usize)->bool;#[doc=" Returns the k-th bit and the number of ones before it."]fn access_rank1(&self,k:usize)->(bool,usize){(self.access(k),self.rank1(k))}#[doc=" the number of 1 in [0, k)"]fn rank1(&self,k:usize)->usize{(0..k).filter(|&i|self.access(i)).count()}#[doc=" the number of 0 in [0, k)"]fn rank0(&self,k:usize)->usize{k-self.rank1(k)}#[doc=" index of k-th 1"]fn select1(&self,k:usize)->Option{let n=self.bit_length();if self.rank1(n)<=k{return None;}let(mut l,mut r)=(0,n);while r-l>1{let m=l.midpoint(r);if self.rank1(m)<=k{l=m;}else{r=m;}}Some(l)}#[doc=" index of k-th 0"]fn select0(&self,k:usize)->Option{let n=self.bit_length();if self.rank0(n)<=k{return None;}let(mut l,mut r)=(0,n);while r-l>1{let m=l.midpoint(r);if self.rank0(m)<=k{l=m;}else{r=m;}}Some(l)}}macro_rules!impl_rank_select_for_bits{($($t:ty)*)=>{$(impl RankSelectDictionaries for$t{fn bit_length(&self)->usize{const WORD_SIZE:usize=(0 as$t).count_zeros()as usize;WORD_SIZE}fn access(&self,k:usize)->bool{const WORD_SIZE:usize=(0 as$t).count_zeros()as usize;if kusize{const WORD_SIZE:usize=(0 as$t).count_zeros()as usize;if kusize{let count=word.count_ones()as usize;debug_assert!(rank>(u64::BITS as usize-width);let count=(word&mask).count_ones()as usize;if rank>=width;rank-=count;offset+=width;}width/=2;}offset}#[cfg(target_arch="x86_64")]#[allow(unsafe_op_in_unsafe_fn)]mod simd{#[cfg(target_arch="x86_64")]use std::arch::x86_64::_pdep_u64;#[target_feature(enable="bmi2")]#[inline]#[cfg(target_arch="x86_64")]pub unsafe fn select_word(word:u64,rank:usize)->usize{_pdep_u64(1<,len:usize,sum:usize,select_samples:[Vec;2]}impl BitVector{const WORD_SIZE:usize=u64::BITS as usize;#[doc=" Builds a bit vector from low-bit-first words. `words.len()` must equal"]#[doc=" `len.div_ceil(64)`. Unused high bits in the final word are ignored."]pub fn from_words(words:&[u64],len:usize)->Self{assert_eq!(words.len(),len.div_ceil(Self::WORD_SIZE));let mut sum=0;let mut blocks=Vec::with_capacity(len/Self::WORD_SIZE+1);for(i,&bits)in words.iter().enumerate(){let count=(len-i*Self::WORD_SIZE).min(Self::WORD_SIZE);let bits=if count==Self::WORD_SIZE{bits}else{bits&((1u64<,len:usize,sum:usize)->Self{let mut select_samples=[Vec::new(),Vec::new()];for(i,block)in blocks.iter().enumerate().take(len.div_ceil(Self::WORD_SIZE)){let start=[i*Self::WORD_SIZE-block.rank,block.rank];let end1=blocks.get(i+1).map_or(sum,|next|next.rank);let end=[((i+1)*Self::WORD_SIZE).min(len)-end1,end1];for bit in 0..2{if start[bit].div_ceil(256)!=end[bit].div_ceil(256){select_samples[bit].push(i);}}}Self{blocks,len,sum,select_samples}}pub fn with_capacity(bits:usize)->Self{let mut blocks=Vec::with_capacity(bits.div_ceil(Self::WORD_SIZE)+1);blocks.push(BitVectorBlock{bits:0,rank:0});Self{blocks,len:0,sum:0,select_samples:[Vec::new(),Vec::new()]}}pub fn push(&mut self,bit:bool){let word=self.len/Self::WORD_SIZE;let rank=if bit{self.sum}else{self.len-self.sum};if rank.is_multiple_of(256){self.select_samples[bit as usize].push(word);}self.blocks[word].bits|=(bit as u64)<<(self.len%Self::WORD_SIZE);self.sum+=bit as usize;self.len+=1;if self.len.is_multiple_of(Self::WORD_SIZE){self.blocks.push(BitVectorBlock{bits:0,rank:self.sum});}}#[doc=" Words paired with the number of ones preceding each word. The last block"]#[doc=" is partial, or an empty sentinel when the bit length is a multiple of 64."]pub fn blocks(&self)->&[BitVectorBlock]{&self.blocks}#[doc=" Returns the position of the zero-based occurrence `rank` in `bits`."]#[doc=" `rank` must be less than the number of set bits."]#[inline]pub fn select_word(bits:u64,rank:usize)->usize{#[cfg(target_arch="x86_64")]if is_x86_feature_detected!("bmi2"){return unsafe{simd::select_word(bits,rank)};}select_word_scalar(bits,rank)}}impl RankSelectDictionaries for BitVector{fn bit_length(&self)->usize{self.len}#[inline]fn access(&self,k:usize)->bool{self.blocks[k/Self::WORD_SIZE].bits&(1u64<<(k%Self::WORD_SIZE))!=0}#[inline]fn access_rank1(&self,k:usize)->(bool,usize){let block=&self.blocks[k/Self::WORD_SIZE];let offset=k%Self::WORD_SIZE;(block.bits&(1u64<usize{self.access_rank1(k).1}fn select1(&self,k:usize)->Option{if k>=self.sum{return None;}let sample=k/256;let start=self.select_samples[1][sample];let end=self.select_samples[1].get(sample+1).map_or(self.blocks.len(),|&word|word+1);let word=start+self.blocks[start..end].partition_point(|block|block.rank<=k)-1;let rank=k-self.blocks[word].rank;Some(word*Self::WORD_SIZE+Self::select_word(self.blocks[word].bits,rank))}fn select0(&self,k:usize)->Option{if k>=self.len-self.sum{return None;}let sample=k/256;let mut word=self.select_samples[0][sample];let end=self.select_samples[0].get(sample+1).map_or(self.blocks.len(),|&word|word+1);let mut size=end-word;while size>1{let half=size/2;let middle=word+half;word=if middle*Self::WORD_SIZE-self.blocks[middle].rank<=k{middle}else{word};size-=half;}let rank=k-(word*Self::WORD_SIZE-self.blocks[word].rank);Some(word*Self::WORD_SIZE+Self::select_word(!self.blocks[word].bits,rank))}}impl FromIteratorfor BitVector{fn from_iter>(iter:I)->Self{let iter=iter.into_iter();let mut blocks=Vec::with_capacity(iter.size_hint().0/Self::WORD_SIZE+1);let mut len=0usize;let mut sum=0;let mut iter=iter.fuse();while let Some(first)=iter.next(){let mut bits=first as u64;let mut count=1;for(i,bit)in iter.by_ref().take(Self::WORD_SIZE-1).enumerate(){bits|=(bit as u64)<<(i+1);count+=1;}blocks.push(BitVectorBlock{bits,rank:sum});len+=count;sum+=bits.count_ones()as usize;}if len.is_multiple_of(Self::WORD_SIZE){blocks.push(BitVectorBlock{bits:0,rank:sum});}Self::from_blocks(blocks,len,sum)}}} // codesnip-guard: RangeMinimumQuery #[cfg_attr(any(),rustfmt::skip)]pub use self::range_minimum_query::RangeMinimumQuery;#[cfg_attr(any(),rustfmt::skip)]mod range_minimum_query{const BLOCK_SIZE:usize=64;#[derive(Clone,Debug)]pub struct RangeMinimumQuery{data:Vec,suffix:Vec,prefix:Vec,table:Vec,blocks:usize}implRangeMinimumQuerywhere T:Ord+Copy{pub fn new(data:Vec)->Self{let n=data.len();let blocks=n.div_ceil(BLOCK_SIZE);if blocks==0{return Self{data,suffix:vec![],prefix:vec![],table:vec![],blocks};}let levels=usize::BITS as usize-blocks.leading_zeros()as usize;let mut table=vec![data[0];levels*blocks];let mut prefix=Vec::with_capacity(n);let mut minimum=data[0];for(i,&value)in data.iter().enumerate(){minimum=if i%BLOCK_SIZE==0{value}else{minimum.min(value)};prefix.push(minimum);if i%BLOCK_SIZE==BLOCK_SIZE-1||i+1==n{table[i/BLOCK_SIZE]=minimum;}}let mut suffix=data.clone();for block in suffix.chunks_mut(BLOCK_SIZE){minimum=*block.last().unwrap();for value in block.iter_mut().rev(){minimum=minimum.min(*value);*value=minimum;}}for level in 1..levels{let current=level*blocks;let previous=current-blocks;let half=1<<(level-1);for i in 0..blocks-(1<T{let r=r-1;let left_block=l/BLOCK_SIZE;let right_block=r/BLOCK_SIZE;if left_block+1where M:Monoid{n:usize,seg:Vec}implClone for SegmentTreewhere M:Monoid{fn clone(&self)->Self{Self{n:self.n,seg:self.seg.clone()}}}implDebug for SegmentTreewhere M:Monoid{fn fmt(&self,f:&mut Formatter<'_>)->fmt::Result{f.debug_struct("SegmentTree").field("n",&self.n).field("seg",&self.seg).finish()}}implSegmentTreewhere M:Monoid{pub fn new(n:usize)->Self{let seg=vec![M::unit();2*n];Self{n,seg}}pub fn from_vec(v:Vec)->Self{let n=v.len();let mut seg=vec![M::unit();2*n];for(i,x)in v.into_iter().enumerate(){seg[n+i]=x;}for i in(1..n).rev(){seg[i]=M::operate(&seg[2*i],&seg[2*i+1]);}Self{n,seg}}pub fn set(&mut self,k:usize,x:M::T){assert!(k0{self.seg[k]=M::operate(&self.seg[2*k],&self.seg[2*k+1]);k/=2;}}pub fn clear(&mut self,k:usize){self.set(k,M::unit());}pub fn update(&mut self,k:usize,x:M::T){assert!(k0{self.seg[k]=M::operate(&self.seg[2*k],&self.seg[2*k+1]);k/=2;}}pub fn get(&self,k:usize)->M::T{assert!(k(&self,range:R)->M::T where R:RangeBounds{let range=range.to_range_bounded(0,self.n).expect("invalid range");let mut l=range.start+self.n;let mut r=range.end+self.n;let mut vl=M::unit();let mut vr=M::unit();while l(&self,mut pos:usize,mut acc:M::T,mut pred:P)->(usize,M::T)where P:FnMut(&M::T)->bool{while pos(&self,mut pos:usize,mut acc:M::T,mut pred:P)->(usize,M::T)where P:FnMut(&M::T)->bool{while pos(&self,left:usize,mut pred:P)->usize where P:FnMut(&M::T)->bool{let mut l=left+self.n;let r=2*self.n;let mut k=0usize;let mut acc=M::unit();while l>k{if l&1!=0{let nacc=M::operate(&acc,&self.seg[l]);if!pred(&nacc){return self.partition_point_perfect(l,acc,pred).0;}acc=nacc;l+=1;}l>>=1;k+=1;}for k in(0..k).rev(){let r=r>>k;if r&1!=0{let nacc=M::operate(&acc,&self.seg[r-1]);if!pred(&nacc){return self.partition_point_perfect(r-1,acc,pred).0;}acc=nacc;}}self.n}pub fn rpartition_point_acc

(&self,right:usize,mut pred:P)->usize where P:FnMut(&M::T)->bool{let mut l=self.n;let mut r=right+self.n;let mut c=0usize;let mut k=0usize;let mut acc=M::unit();while l>>k>=1;k+=1;}for k in(0..k).rev(){if c&1!=0{l-=1<>k;let nacc=M::operate(&self.seg[l],&acc);if!pred(&nacc){return self.rpartition_point_perfect(l,acc,pred).0+1;}acc=nacc;}c>>=1;}0}pub fn as_slice(&self)->&[M::T]{&self.seg[self.n..]}}implSegmentTreewhere M:AbelianMonoid{pub fn fold_all(&self)->M::T{self.seg[1].clone()}}} // codesnip-guard: WaveletMatrix #[cfg_attr(any(),rustfmt::skip)]pub use self::wavelet_matrix::{WaveletMatrix,WaveletMatrixPointAdd};#[cfg_attr(any(),rustfmt::skip)]mod wavelet_matrix{use super::{AbelianGroup,BinaryIndexedTree,BitVector,Compressor,OrderedCompressor,RankSelectDictionaries,VecCompress};use std::{mem::{self,MaybeUninit},ops::Range};const ACCESS:u8=0;const RANK:u8=1;const QUANTILE:u8=2;const RANK_LESSTHAN:u8=3;#[derive(Debug,Clone)]#[repr(C)]struct WaveletMatrixQuadBlock{lo:u64,hi:u64,rank:[u32;4]}#[derive(Debug,Clone)]struct WaveletMatrixQuadVector{blocks:Vec,starts:[usize;4],select_samples:[Vec;4]}impl WaveletMatrixQuadVector{fn from_words(low:&[u64],high:Option<&[u64]>,len:usize)->Self{let mut blocks=Vec::with_capacity(len/64+1);let mut rank=[0;4];for(word,&lo)in low.iter().enumerate(){let count=(len-word*64).min(64);let hi=high.map_or(0,|high|high[word]);blocks.push(WaveletMatrixQuadBlock{lo,hi,rank});let low=lo.count_ones();let high=hi.count_ones();let both=(lo&hi).count_ones();rank[0]+=count as u32-low-(high-both);rank[1]+=low-both;rank[2]+=high-both;rank[3]+=both;}if len.is_multiple_of(64){blocks.push(WaveletMatrixQuadBlock{lo:0,hi:0,rank});}let starts=[0,(rank[0]as usize+rank[2]as usize),rank[0]as usize,rank[0]as usize+rank[2]as usize+rank[1]as usize];let mut select_samples:[Vec;4]=std::array::from_fn(|_|Vec::new());for(i,block)in blocks.iter().enumerate().take(len.div_ceil(64)){for digit in 0..4{let start=block.rank[digit]as usize;let end=blocks.get(i+1).map_or(rank[digit],|next|next.rank[digit])as usize;if start.div_ceil(128)!=end.div_ceil(128){select_samples[digit].push(i);}}}Self{blocks,starts,select_samples}}#[inline]fn ranks(&self,position:usize)->[usize;4]{let block=&self.blocks[position/64];let mask=!(u64::MAX<<(position%64));let low=(block.lo&mask).count_ones()as usize;let high=(block.hi&mask).count_ones()as usize;let both=(block.lo&block.hi&mask).count_ones()as usize;[block.rank[0]as usize+position%64-low-(high-both),block.rank[1]as usize+low-both,block.rank[2]as usize+high-both,block.rank[3]as usize+both]}#[inline]fn rank(&self,digit:usize,position:usize)->usize{let block=&self.blocks[position/64];let low=if digit&1!=0{block.lo}else{!block.lo};let high=if digit&2!=0{block.hi}else{!block.hi};block.rank[digit]as usize+(low&high&!(u64::MAX<<(position%64))).count_ones()as usize}fn select(&self,digit:usize,k:usize)->usize{let sample=k/128;let start=self.select_samples[digit][sample];let end=self.select_samples[digit].get(sample+1).map_or(self.blocks.len(),|&word|word+1);let word=start+self.blocks[start..end].partition_point(|block|block.rank[digit]as usize<=k)-1;let block=&self.blocks[word];let lo=if digit&1!=0{block.lo}else{!block.lo};let hi=if digit&2!=0{block.hi}else{!block.hi};word*64+BitVector::select_word(lo&hi,k-block.rank[digit]as usize)}#[inline]fn access_rank(&self,position:usize)->(usize,usize){let block=&self.blocks[position/64];let offset=position%64;let digit=((block.lo>>offset)&1)as usize|(((block.hi>>offset)&1)as usize)<<1;(digit,self.rank(digit,position))}}#[cfg(target_arch="x86_64")]mod simd{#![allow(unsafe_op_in_unsafe_fn)]#[cfg(target_arch="x86_64")]use super::BitVector;#[cfg(target_arch="x86_64")]use std::arch::x86_64::*;#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn pack_words(indices:&[u32],bit:usize)->Vec{let mut result=Vec::with_capacity(indices.len().div_ceil(64));let shift=_mm_cvtsi32_si128((31-bit)as i32);for chunk in indices.chunks(64){let mut word=0u64;let mut i=0;while i+8<=chunk.len(){let v=_mm256_loadu_si256(chunk.as_ptr().add(i).cast());word|=(_mm256_movemask_ps(_mm256_castsi256_ps(_mm256_sll_epi32(v,shift)))as u64)<>bit)&1)as u64)<<(i+j);}result.push(word);}result}#[target_feature(enable="avx2")]#[cfg(target_arch="x86_64")]pub unsafe fn partition_avx2(indices:&[u32],words:&[u64],mut one:usize,next:&mut[u32]){const PERM:[[i32;8];256]={let mut table=[[0;8];256];let mut mask=0;while mask<256{let mut pos=0;let mut bit=0;while bit<2{let mut i=0;while i<8{if(mask>>i)&1==bit{table[mask][pos]=i;pos+=1;}i+=1;}bit+=1;}mask+=1;}table};let order=_mm256_setr_epi32(0,1,2,3,4,5,6,7);let mut zero=0;for(chunk,&word)in indices.chunks(64).zip(words){if word==0||word==u64::MAX{let at=if word==0{&mut zero}else{&mut one};next[*at..*at+chunk.len()].copy_from_slice(chunk);*at+=chunk.len();continue;}let mut i=0;while i+8<=chunk.len(){let value=_mm256_loadu_si256(chunk.as_ptr().add(i).cast());let mask=((word>>i)&255)as usize;let ones=mask.count_ones()as usize;let zeros=8-ones;let packed=_mm256_permutevar8x32_epi32(value,_mm256_loadu_si256(PERM[mask].as_ptr().cast()));let rotated=_mm256_permutevar8x32_epi32(packed,_mm256_add_epi32(order,_mm256_set1_epi32(zeros as i32)));_mm256_maskstore_epi32(next.as_mut_ptr().add(zero).cast(),_mm256_cmpgt_epi32(_mm256_set1_epi32(zeros as i32),order),packed);_mm256_maskstore_epi32(next.as_mut_ptr().add(one).cast(),_mm256_cmpgt_epi32(_mm256_set1_epi32(ones as i32),order),rotated);zero+=zeros;one+=ones;i+=8;}for(j,&value)in chunk[i..].iter().enumerate(){let bit=(word>>(i+j))&1!=0;next[if bit{one}else{zero}]=value;zero+=!bit as usize;one+=bit as usize;}}}#[target_feature(enable="avx2")]#[inline]#[cfg(target_arch="x86_64")]unsafe fn popcount_avx2(value:__m256i)->__m256i{let lookup=_mm256_setr_epi8(0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4);let mask=_mm256_set1_epi8(15);let low=_mm256_shuffle_epi8(lookup,_mm256_and_si256(value,mask));let high=_mm256_shuffle_epi8(lookup,_mm256_and_si256(_mm256_srli_epi16::<4>(value),mask));_mm256_sad_epu8(_mm256_add_epi8(low,high),_mm256_setzero_si256())}#[target_feature(enable="avx512f")]#[inline]#[cfg(target_arch="x86_64")]unsafe fn greater_avx512(a:__m512i,b:__m512i)->__m512i{_mm512_maskz_set1_epi64(_mm512_cmpgt_epi64_mask(a,b),-1)}#[target_feature(enable="avx512f")]#[inline]#[cfg(target_arch="x86_64")]unsafe fn gather_avx512(base:*const i64,index:__m512i)->__m512i{_mm512_i64gather_epi64::(index,base)}#[cfg(target_arch="x86_64")]macro_rules!rank_lessthan{($name:ident,$feature:literal,$lanes:literal,$load:ident,$store:ident,$set:ident,$gather:ident,$add:ident,$sub:ident,$and:ident,$andnot:ident,$or:ident,$left:ident,$right:ident,$gt:ident,$popcount:ident)=>{#[target_feature(enable=$feature)]pub unsafe fn$name(vectors:&[BitVector],zeros:&[usize],states:&mut[[usize;4]]){for chunk in states.chunks_exact_mut($lanes){let mut starts=[0u64;$lanes];let mut ends=[0u64;$lanes];let mut keys=[0u64;$lanes];let mut values=[0u64;$lanes];for(i,state)in chunk.iter().enumerate(){starts[i]=state[0]as u64;ends[i]=state[1]as u64;keys[i]=state[2]as u64;}let mut start=$load(starts.as_ptr().cast());let mut end=$load(ends.as_ptr().cast());let key=$load(keys.as_ptr().cast());let mut value=$set(0);let one=$set(1);for(level,vector)in vectors.iter().enumerate(){let d=vectors.len()-1-level;let base=vector.blocks().as_ptr().cast::();let rank1=|position|{let offset=$and(position,$set(63));let index=$left($right(position,$set(6)),one);let bits=$gather::<8>(base,index);let prefix=$gather::<8>(base,$add(index,one));$add(prefix,$popcount($and(bits,$sub($left(one,offset),one))))};let rank=rank1(start);let end_rank=rank1(end);let start0=$sub(start,rank);let end0=$sub(end,end_rank);let count=$sub(end0,start0);let mask=$gt($and($right(key,$set(d as i64)),one),$set(0));let zero=$set(zeros[level]as i64);start=$or($andnot(mask,start0),$and(mask,$add(zero,rank)));end=$or($andnot(mask,end0),$and(mask,$add(zero,end_rank)));value=$add(value,$and(mask,count));}$store(values.as_mut_ptr().cast(),value);for(i,state)in chunk.iter_mut().enumerate(){state[3]=values[i]as usize;}}}};}#[cfg(target_arch="x86_64")]rank_lessthan!(rank_lessthan_avx2,"avx2",4,_mm256_loadu_si256,_mm256_storeu_si256,_mm256_set1_epi64x,_mm256_i64gather_epi64,_mm256_add_epi64,_mm256_sub_epi64,_mm256_and_si256,_mm256_andnot_si256,_mm256_or_si256,_mm256_sllv_epi64,_mm256_srlv_epi64,_mm256_cmpgt_epi64,popcount_avx2);#[cfg(target_arch="x86_64")]rank_lessthan!(rank_lessthan_avx512,"avx512f,avx512vpopcntdq",8,_mm512_loadu_si512,_mm512_storeu_si512,_mm512_set1_epi64,gather_avx512,_mm512_add_epi64,_mm512_sub_epi64,_mm512_and_si512,_mm512_andnot_si512,_mm512_or_si512,_mm512_sllv_epi64,_mm512_srlv_epi64,greater_avx512,_mm512_popcnt_epi64);#[target_feature(enable="avx512f,avx512vpopcntdq")]#[cfg(target_arch="x86_64")]pub unsafe fn quad_avx512(layers:&[super::WaveletMatrixQuadVector],states:&mut[[usize;4]]){for chunk in states.as_chunks_mut::<8>().0{let mut starts=[0u64;8];let mut ends=[0u64;8];let mut keys=[0u64;8];let mut result=[0u64;8];for(i,s)in chunk.iter().enumerate(){starts[i]=s[0]as u64;ends[i]=s[1]as u64;keys[i]=s[2]as u64;}let mut start=_mm512_loadu_si512(starts.as_ptr().cast());let mut end=_mm512_loadu_si512(ends.as_ptr().cast());let mut key=_mm512_loadu_si512(keys.as_ptr().cast());let mut code=_mm512_set1_epi64(0);let one=_mm512_set1_epi64(1);macro_rules!blend{($m:expr,$x:expr,$y:expr)=>{_mm512_or_si512(_mm512_andnot_si512($m,$x),_mm512_and_si512($m,$y))};}for layer in layers{let base=layer.blocks.as_ptr().cast::();let ranks=|pos|{let offset=_mm512_and_si512(pos,_mm512_set1_epi64(63));let mask=_mm512_sub_epi64(_mm512_sllv_epi64(one,offset),one);let i=_mm512_sllv_epi64(_mm512_srlv_epi64(pos,_mm512_set1_epi64(6)),_mm512_set1_epi64(2));let lo=_mm512_and_si512(gather_avx512::<8>(base,i),mask);let hi=_mm512_and_si512(gather_avx512::<8>(base,_mm512_add_epi64(i,one)),mask);let a=_mm512_popcnt_epi64(lo);let b=_mm512_popcnt_epi64(hi);let c=_mm512_popcnt_epi64(_mm512_and_si512(lo,hi));let p01=gather_avx512::<8>(base,_mm512_add_epi64(i,_mm512_set1_epi64(2)));let p23=gather_avx512::<8>(base,_mm512_add_epi64(i,_mm512_set1_epi64(3)));let r0=_mm512_add_epi64(_mm512_and_si512(p01,_mm512_set1_epi64(u32::MAX as i64)),_mm512_add_epi64(_mm512_sub_epi64(_mm512_sub_epi64(offset,a),b),c));let r1=_mm512_add_epi64(_mm512_srlv_epi64(p01,_mm512_set1_epi64(32)),_mm512_sub_epi64(a,c));let r2=_mm512_add_epi64(_mm512_and_si512(p23,_mm512_set1_epi64(u32::MAX as i64)),_mm512_sub_epi64(b,c));let r3=_mm512_sub_epi64(_mm512_sub_epi64(_mm512_sub_epi64(pos,r0),r1),r2);[r0,r1,r2,r3]};let l=ranks(start);let r=ranks(end);let low_count=_mm512_sub_epi64(_mm512_add_epi64(r[0],r[1]),_mm512_add_epi64(l[0],l[1]));let high=greater_avx512(key,_mm512_sub_epi64(low_count,one));key=_mm512_sub_epi64(key,_mm512_and_si512(high,low_count));let l0=blend!(high,l[0],l[2]);let r0=blend!(high,r[0],r[2]);let l1=blend!(high,l[1],l[3]);let r1=blend!(high,r[1],r[3]);let count0=_mm512_sub_epi64(r0,l0);let low=greater_avx512(key,_mm512_sub_epi64(count0,one));key=_mm512_sub_epi64(key,_mm512_and_si512(low,count0));let base0=blend!(high,_mm512_set1_epi64(layer.starts[0]as i64),_mm512_set1_epi64(layer.starts[2]as i64));let base1=blend!(high,_mm512_set1_epi64(layer.starts[1]as i64),_mm512_set1_epi64(layer.starts[3]as i64));let offset=blend!(low,base0,base1);start=_mm512_add_epi64(offset,blend!(low,l0,l1));end=_mm512_add_epi64(offset,blend!(low,r0,r1));code=_mm512_or_si512(_mm512_sllv_epi64(code,_mm512_set1_epi64(2)),_mm512_or_si512(_mm512_and_si512(high,_mm512_set1_epi64(2)),_mm512_and_si512(low,one)));}_mm512_storeu_si512(result.as_mut_ptr().cast(),code);for(i,s)in chunk.iter_mut().enumerate(){s[3]=result[i]as usize;}}}}#[derive(Debug,Clone)]pub struct WaveletMatrix{len:usize,bit_length:usize,zeros:Vec,bit_vectors:Vec,quad_vectors:Vec,compress:VecCompress,#[cfg(target_arch="x86_64")]backend:super::SimdBackend}implWaveletMatrixwhere T:Ord+Clone{pub fn new(v:Vec)->Self{if v.len()<=u32::MAX as usize{#[cfg(target_arch="x86_64")]let backend=super::simd_backend();Self::from_values(v,|i|i as u32,|i|i as usize,|indices,d|{#[cfg(target_arch="x86_64")]if backend!=super::SimdBackend::Scalar&&is_x86_feature_detected!("avx2"){return unsafe{simd::pack_words(indices,d)};}Self::pack_words(indices,|i|(i>>d)&1!=0)},|indices,words,zeros,next|{#[cfg(target_arch="x86_64")]if backend!=super::SimdBackend::Scalar&&is_x86_feature_detected!("avx2"){unsafe{simd::partition_avx2(indices,words,zeros,next)};return;}Self::partition(indices,words,zeros,next);})}else{Self::from_values(v,|i|i,|i|i,|indices,d|Self::pack_words(indices,|i|(i>>d)&1!=0),Self::partition)}}fn pack_words(indices:&[I],bit:impl Fn(I)->bool)->Vec{indices.chunks(64).map(|chunk|{chunk.iter().enumerate().fold(0,|word,(i,&index)|word|((bit(index)as u64)<(indices:&[I],words:&[u64],mut one:usize,next:&mut[I]){let mut zero=0;for(chunk,&word)in indices.chunks(64).zip(words){if word==0{next[zero..zero+chunk.len()].copy_from_slice(chunk);zero+=chunk.len();}else if word==u64::MAX{next[one..one+chunk.len()].copy_from_slice(chunk);one+=chunk.len();}else{for(i,&index)in chunk.iter().enumerate(){let bit=(word>>i)&1!=0;next[if bit{one}else{zero}]=index;zero+=!bit as usize;one+=bit as usize;}}}}fn from_values(v:Vec,code:impl Fn(usize)->I,index:impl Fn(I)->usize,pack:impl Fn(&[I],usize)->Vec,partition:impl Fn(&[I],&[u64],usize,&mut[I]))->Self{let len=v.len();let mut sorted:Vec<_>=v.into_iter().enumerate().map(|(i,value)|(value,code(i))).collect();sorted.sort_unstable_by(|a,b|a.0 .cmp(&b.0));let mut values=Vec::with_capacity(len);let mut indices=vec![code(0);len];for(value,i)in sorted{if values.last().is_none_or(|last|last!=&value){values.push(value);}indices[index(i)]=code(values.len()-1);}let compress=VecCompress::from_sorted_unique(values);let bit_length=usize::BITS as usize-compress.size().leading_zeros()as usize;let mut bit_vectors=Vec::with_capacity(bit_length);let mut zeros=Vec::with_capacity(bit_length);let quad_bits=usize::BITS as usize-compress.size().saturating_sub(1).leading_zeros()as usize;let mut quad_vectors=Vec::with_capacity(quad_bits.div_ceil(2));let mut next=indices.clone();for d in(0..bit_length).rev(){let words=pack(&indices,d);if len<=u32::MAX as usize&&d{super::SimdBackend::Avx2}backend=>backend,}}}pub fn new_with_init(v:Vec,mut f:F)->Self where F:FnMut(usize,usize,T){let this=Self::new(v.clone());if!this.quad_vectors.is_empty(){let bits=usize::BITS as usize-this.compress.size().saturating_sub(1).leading_zeros()as usize;for(mut k,value)in v.into_iter().enumerate(){if this.bit_length>bits{f(this.bit_length-1,k,value.clone());}for(level,vector)in this.quad_vectors.iter().enumerate(){let d=(this.quad_vectors.len()-level-1)*2;let(digit,rank)=vector.access_rank(k);if d+1usize{self.bit_length-1-d}fn rank1(&self,level:usize,k:usize)->usize{self.bit_vectors[level].rank1(k)}fn reorder(&self,level:usize,current:Vec)->Vec{assert_eq!(current.len(),self.len);let mut next=Vec::with_capacity(self.len);next.resize_with(self.len,MaybeUninit::uninit);let mut zero=0;let mut one=self.zeros[level];let mut current=current.into_iter();for block in self.bit_vectors[level].blocks(){let count=current.len().min(64);if block.bits==0||block.bits==u64::MAX{let offset=if block.bits==0{&mut zero}else{&mut one};for(slot,value)in next[*offset..*offset+count].iter_mut().zip(current.by_ref().take(count)){slot.write(value);}*offset+=count;}else{for(i,value)in current.by_ref().take(count).enumerate(){let bit=(block.bits>>i)&1!=0;next[if bit{one}else{zero}].write(value);zero+=!bit as usize;one+=bit as usize;}}}unsafe{let mut next=mem::ManuallyDrop::new(next);Vec::from_raw_parts(next.as_mut_ptr().cast(),next.len(),next.capacity())}}fn range_by_index(&self,idx:usize,mut range:Range)->Range{if!self.quad_vectors.is_empty(){for(level,vector)in self.quad_vectors.iter().enumerate(){if range.is_empty(){break;}let digit=(idx>>((self.quad_vectors.len()-level-1)*2))&3;range=vector.starts[digit]+vector.rank(digit,range.start)..vector.starts[digit]+vector.rank(digit,range.end);}return range;}for d in(0..self.bit_length-self.compress.size().is_power_of_two()as usize).rev(){if range.is_empty(){break;}let level=self.level(d);let start1=self.rank1(level,range.start);let end1=self.rank1(level,range.end);if((idx>>d)&1)!=0{range.start=self.zeros[level]+start1;range.end=self.zeros[level]+end1;}else{range.start-=start1;range.end-=end1;}}range}fn batch(&self,states:&mut[[usize;4];16],count:usize){#[cfg(target_arch="x86_64")]if!cfg!(target_feature="popcnt")&&is_x86_feature_detected!("popcnt"){unsafe{self.batch_popcnt::(states,count);}return;}self.batch_inner::(states,count);}#[cfg(target_arch="x86_64")]#[target_feature(enable="popcnt")]unsafe fn batch_popcnt(&self,states:&mut[[usize;4];16],count:usize){self.batch_inner::(states,count);}#[inline(always)]fn batch_inner(&self,states:&mut[[usize;4];16],count:usize){if OP==RANK_LESSTHAN&&self.compress.size()<=1{for state in&mut states[..count]{state[3]=if state[2]==0{0}else{state[1]-state[0]};}return;}if!self.quad_vectors.is_empty()&&matches!(OP,ACCESS|RANK|QUANTILE){#[cfg(target_arch="x86_64")]if OP==QUANTILE&&count>=8&&self.backend==super::SimdBackend::Avx512{unsafe{simd::quad_avx512(&self.quad_vectors,&mut states[..count.next_multiple_of(8)]);}return;}for(level,vector)in self.quad_vectors.iter().enumerate(){let d=(self.quad_vectors.len()-level-1)*2;#[cfg(target_arch="x86_64")]let next=self.quad_vectors.get(level+1).filter(|_|self.len>=1048576||(OP==ACCESS&&self.len>=262144));for state in&mut states[..count]{if OP==ACCESS{let(digit,rank)=vector.access_rank(state[0]);state[0]=vector.starts[digit]+rank;state[3]=state[3]*4+digit;}else if OP==RANK{let digit=(state[2]>>d)&3;state[0]=vector.starts[digit]+vector.rank(digit,state[0]);state[1]=vector.starts[digit]+vector.rank(digit,state[1]);}else{let start=vector.ranks(state[0]);let end=vector.ranks(state[1]);let prefix=[0,end[0]-start[0],end[0]+end[1]-start[0]-start[1],state[1]-state[0]-(end[3]-start[3])];let digit=(state[2]>=prefix[1])as usize+(state[2]>=prefix[2])as usize+(state[2]>=prefix[3])as usize;state[2]-=prefix[digit];state[0]=vector.starts[digit]+start[digit];state[1]=vector.starts[digit]+end[digit];state[3]=state[3]*4+digit;}#[cfg(target_arch="x86_64")]if let Some(next)=next{unsafe{std::arch::x86_64::_mm_prefetch::<{std::arch::x86_64::_MM_HINT_T0}>(next.blocks.as_ptr().add(state[0]/64).cast());if OP!=ACCESS{std::arch::x86_64::_mm_prefetch::<{std::arch::x86_64::_MM_HINT_T0}>(next.blocks.as_ptr().add(state[1]/64).cast());}}}}}return;}let first=(matches!(OP,ACCESS|RANK|QUANTILE)&&self.compress.size().is_power_of_two())as usize;#[cfg(target_arch="x86_64")]if OP==RANK_LESSTHAN&&count>=8&&self.len<=u32::MAX as usize{unsafe{match self.backend{super::SimdBackend::Avx512=>{simd::rank_lessthan_avx512(&self.bit_vectors[first..],&self.zeros[first..],&mut states[..count.next_multiple_of(8)]);return;}super::SimdBackend::Avx2 if self.len<1048576=>{simd::rank_lessthan_avx2(&self.bit_vectors[first..],&self.zeros[first..],&mut states[..count.next_multiple_of(4)]);return;}_=>{}}}}for d in(0..self.bit_length-first).rev(){let level=self.level(d);for state in&mut states[..count]{let(bit,start1)=self.bit_vectors[level].access_rank1(state[0]);let start0=state[0]-start1;let end1=if OP==ACCESS{0}else{self.rank1(level,state[1])};let end0=if OP==ACCESS{0}else{state[1]-end1};let count0=if OP==ACCESS{0}else{end0-start0};let bit=match OP{ACCESS=>bit,RANK|RANK_LESSTHAN=>(state[2]>>d)&1!=0,_=>state[2]>=count0,};state[0]=if bit{self.zeros[level]+start1}else{start0};if OP!=ACCESS{state[1]=if bit{self.zeros[level]+end1}else{end0};}if OP==ACCESS||OP==QUANTILE{state[3]|=(bit as usize)<T{if!self.quad_vectors.is_empty(){let mut index=0;for vector in&self.quad_vectors{let(digit,rank)=vector.access_rank(k);index=index*4+digit;k=vector.starts[digit]+rank;}return self.compress.values()[index].clone();}let mut idx=0;for d in(0..self.bit_length-self.compress.size().is_power_of_two()as usize).rev(){let level=self.level(d);let(bit,rank1)=self.bit_vectors[level].access_rank1(k);idx|=(bit as usize)<)->Vec{let indices:Vec<_>=indices.into_iter().collect();let mut result=Vec::with_capacity(indices.len());for indices in indices.chunks(16){let mut states=[[0;4];16];for(state,&index)in states.iter_mut().zip(indices){assert!(index(&mut states,indices.len());result.extend(states[..indices.len()].iter().map(|state|self.compress.values()[state[3]].clone()));}result}#[doc=" the number of val in range"]pub fn rank(&self,val:T,range:Range)->usize{match self.compress.index_exact(&val){Some(idx)=>self.range_by_index(idx,range).len(),None=>0,}}#[doc=" Returns the number of exact matches for each `(value, range)` query."]pub fn rank_batch(&self,queries:impl IntoIterator)>)->Vec{let queries:Vec<_>=queries.into_iter().collect();let mut result=Vec::with_capacity(queries.len());for queries in queries.chunks(16){let mut states=[[0;4];16];for(state,(value,range))in states.iter_mut().zip(queries){assert!(range.start<=range.end&&range.end<=self.len);if let Some(index)=self.compress.index_exact(value){*state=[range.start,range.end,index,0];}}self.batch::(&mut states,queries.len());result.extend(states[..queries.len()].iter().map(|state|state[1]-state[0]));}result}#[doc=" index of k-th val"]pub fn select(&self,val:T,k:usize)->Option{let idx=self.compress.index_exact(&val)?;let range=self.range_by_index(idx,0..self.len);if range.len()<=k{return None;}let mut i=range.start+k;if!self.quad_vectors.is_empty(){for(level,vector)in self.quad_vectors.iter().enumerate().rev(){let digit=(idx>>((self.quad_vectors.len()-level-1)*2))&3;i=vector.select(digit,i-vector.starts[digit]);}return Some(i);}for level in(self.compress.size().is_power_of_two()as usize..self.bit_length).rev(){if i>=self.zeros[level]{i=self.bit_vectors[level].select1(i-self.zeros[level]).unwrap();}else{i=self.bit_vectors[level].select0(i).unwrap();}}Some(i)}#[doc=" get k-th smallest value in range"]pub fn quantile(&self,mut range:Range,mut k:usize)->T{if!self.quad_vectors.is_empty(){return self.quad_quantile(range,k,0,0);}let mut idx=0;for d in(0..self.bit_length-self.compress.size().is_power_of_two()as usize).rev(){let level=self.level(d);let start1=self.rank1(level,range.start);let end1=self.rank1(level,range.end);let start0=range.start-start1;let end0=range.end-end1;let z=end0-start0;let bit=z<=k;k-=if bit{z}else{0};idx|=(bit as usize)<,mut k:usize,level:usize,mut index:usize)->T{for vector in&self.quad_vectors[level..]{let start=vector.ranks(range.start);let end=vector.ranks(range.end);let mut digit=0;while digit<3&&k>=end[digit]-start[digit]{k-=end[digit]-start[digit];digit+=1;}index=index*4+digit;range=vector.starts[digit]+start[digit]..vector.starts[digit]+end[digit];}self.compress.values()[index].clone()}pub fn quantile_batch(&self,queries:impl IntoIterator,usize)>)->Vec{let queries:Vec<_>=queries.into_iter().collect();let mut result=Vec::with_capacity(queries.len());for queries in queries.chunks(16){let mut states=[[0;4];16];for(state,(range,k))in states.iter_mut().zip(queries){assert!(range.start<=range.end&&range.end<=self.len&&*k(&mut states,queries.len());result.extend(states[..queries.len()].iter().map(|state|self.compress.values()[state[3]].clone()));}result}#[doc=" Returns the requested order statistics of one range. `ranks` must be sorted"]#[doc=" in nondecreasing order and every rank must be less than the range length."]pub fn quantiles_sorted(&self,range:Range,ranks:&[usize])->Vec{assert!(range.start<=range.end&&range.end<=self.len);assert!(ranks.windows(2).all(|pair|pair[0]<=pair[1]));assert!(ranks.last().is_none_or(|&k|k,mut k:usize)->T{if!self.quad_vectors.is_empty(){let mut outer=0..self.len;let mut index=0;for vector in&self.quad_vectors{let start=vector.ranks(range.start);let end=vector.ranks(range.end);let outer_start=vector.ranks(outer.start);let outer_end=vector.ranks(outer.end);let mut digit=0;while digit<3{let count=outer_end[digit]-outer_start[digit]-(end[digit]-start[digit]);if k)->usize{let index=self.compress.index_lower_bound(&val);if index==self.compress.size(){return range.len();}if self.len>=262144&&!self.quad_vectors.is_empty(){self.quad_rank_lessthan_index(index,range,0)}else{let bits=usize::BITS as usize-self.compress.size().saturating_sub(1).leading_zeros()as usize;self.rank_lessthan_index(index,range,bits)}}#[doc=" Returns the number of values below each query's threshold."]pub fn rank_lessthan_batch(&self,queries:impl IntoIterator)>)->Vec{let queries:Vec<_>=queries.into_iter().collect();let mut result=Vec::with_capacity(queries.len());for queries in queries.chunks(16){let mut states=[[0;4];16];for(state,(value,range))in states.iter_mut().zip(queries){assert!(range.start<=range.end&&range.end<=self.len);*state=[range.start,range.end,self.compress.index_lower_bound(value),0];}self.batch::(&mut states,queries.len());result.extend(states[..queries.len()].iter().map(|state|state[3]));}result}fn quad_rank_lessthan_index(&self,index:usize,mut range:Range,start_level:usize)->usize{let mut result=0;for(level,vector)in self.quad_vectors.iter().enumerate().skip(start_level){let d=(self.quad_vectors.len()-level-1)*2;if d+2<=index.trailing_zeros()as usize{break;}let digit=(index>>d)&3;let ranks=|position:usize|{let block=&vector.blocks[position/64];let mask=!(u64::MAX<<(position%64));let low=if digit&1==0{!block.lo}else{block.lo};let high=if digit&2==0{!block.hi}else{block.hi};let exact=block.rank[digit]as usize+(low&high&mask).count_ones()as usize;let less=match digit{0=>0,1=>{block.rank[0]as usize+(!(block.lo|block.hi)&mask).count_ones()as usize}2=>{block.rank[0]as usize+block.rank[1]as usize+(!block.hi&mask).count_ones()as usize}_=>position-exact,};(exact,less)};let(start,start_less)=ranks(range.start);let(end,end_less)=ranks(range.end);result+=end_less-start_less;range=vector.starts[digit]+start..vector.starts[digit]+end;}result}fn rank_lessthan_index(&self,idx:usize,mut range:Range,bits:usize)->usize{let mut res=0;for d in(idx.trailing_zeros()as usize..bits).rev(){let level=self.level(d);let start1=self.rank1(level,range.start);let end1=self.rank1(level,range.end);let bit=(idx>>d)&1!=0;let start0=range.start-start1;let end0=range.end-end1;res+=if bit{end0-start0}else{0};range.start=if bit{self.zeros[level]+start1}else{start0};range.end=if bit{self.zeros[level]+end1}else{end0};}res}#[doc=" the number of valrange in range"]pub fn rank_range(&self,valrange:Range,mut range:Range)->usize{let lower=self.compress.index_lower_bound(&valrange.start);let upper=self.compress.index_lower_bound(&valrange.end);if lower>=upper{return 0;}if!self.quad_vectors.is_empty(){if upper==self.compress.size(){return range.len()-self.quad_rank_lessthan_index(lower,range,0);}for(level,vector)in self.quad_vectors.iter().enumerate(){let d=(self.quad_vectors.len()-level-1)*2;let low=(lower>>d)&3;let high=(upper>>d)&3;if low==high{range=vector.starts[low]+vector.rank(low,range.start)..vector.starts[low]+vector.rank(low,range.end);continue;}let start=vector.ranks(range.start);let end=vector.ranks(range.end);let count:usize=(low..high).map(|digit|end[digit]-start[digit]).sum();return count-self.quad_rank_lessthan_index(lower,vector.starts[low]+start[low]..vector.starts[low]+end[low],level+1)+self.quad_rank_lessthan_index(upper,vector.starts[high]+start[high]..vector.starts[high]+end[high],level+1);}return 0;}for d in(0..self.bit_length).rev(){let level=self.level(d);let start1=self.rank1(level,range.start);let end1=self.rank1(level,range.end);let zero_range=range.start-start1..range.end-end1;let one_range=self.zeros[level]+start1..self.zeros[level]+end1;if((lower^upper)>>d)&1!=0{return zero_range.len()-self.rank_lessthan_index(lower,zero_range,d)+self.rank_lessthan_index(upper,one_range,d);}range=if(lower>>d)&1!=0{one_range}else{zero_range};}0}#[doc=" Returns each query's count in a half-open value range."]pub fn rank_range_batch(&self,queries:impl IntoIterator,Range)>)->Vec{let queries:Vec<_>=queries.into_iter().collect();let mut result=Vec::with_capacity(queries.len());for queries in queries.chunks(8){let mut states=[[0;4];16];for(i,(values,range))in queries.iter().enumerate(){assert!(range.start<=range.end&&range.end<=self.len);let lower=self.compress.index_lower_bound(&values.start);let upper=self.compress.index_lower_bound(&values.end);if lower(&mut states,queries.len()*2);result.extend(states[..queries.len()*2].as_chunks::<2>().0 .iter().map(|pair|pair[1][3]-pair[0][3]));}result}pub fn query_less_than(&self,val:T,mut range:Range,mut f:F)where F:FnMut(usize,Range){let idx=self.compress.index_lower_bound(&val);if!self.quad_vectors.is_empty(){if idx==self.compress.size()&&idx.is_power_of_two(){f(self.bit_length-1,range);return;}for(level,vector)in self.quad_vectors.iter().enumerate().take(self.quad_vectors.len().saturating_sub(idx.trailing_zeros()as usize/2)){let d=(self.quad_vectors.len()-level-1)*2;let digit=(idx>>d)&3;let start=vector.ranks(range.start);let end=vector.ranks(range.end);if digit&2!=0{f(d+1,start[0]+start[1]..end[0]+end[1]);}if digit&1!=0{let zero=digit&2;f(d,vector.starts[zero]+start[zero]..vector.starts[zero]+end[zero]);}range=vector.starts[digit]+start[digit]..vector.starts[digit]+end[digit];}return;}for d in(idx.trailing_zeros()as usize..self.bit_length).rev(){let level=self.level(d);let start1=self.rank1(level,range.start);let end1=self.rank1(level,range.end);let start0=range.start-start1;let end0=range.end-end1;if((idx>>d)&1)!=0{f(d,start0..end0);range.start=self.zeros[level]+start1;range.end=self.zeros[level]+end1;}else{range.start=start0;range.end=end0;}}}pub fn build_fold(&self,weights:&[M::T])->WaveletMatrixFold<'_,T,M>where M:AbelianGroup{assert_eq!(weights.len(),self.len);let mut offsets=Vec::with_capacity(self.bit_length);let mut prefix=Vec::with_capacity(self.zeros.iter().map(|&zero|zero+1).sum());let mut current:Vec=weights.to_vec();for level in 0..self.bit_length{current=self.reorder(level,current);offsets.push(prefix.len());let mut acc=M::unit();prefix.push(acc.clone());for w in¤t[..self.zeros[level]]{acc=M::operate(&acc,w);prefix.push(acc.clone());}}WaveletMatrixFold{wavelet_matrix:self,prefix,offsets}}pub fn build_point_add(&self,weights:&[M::T])->WaveletMatrixPointAdd<'_,T,M>where M:AbelianGroup{assert_eq!(weights.len(),self.len);let mut current=weights.to_vec();let mut bits=Vec::with_capacity(self.bit_length);for level in 0..self.bit_length{current=self.reorder(level,current);bits.push(BinaryIndexedTree::from_slice(¤t[..self.zeros[level]]));}WaveletMatrixPointAdd{wavelet_matrix:self,bits}}}pub struct WaveletMatrixPointAdd<'a,T,M>where T:Ord+Clone,M:AbelianGroup{wavelet_matrix:&'a WaveletMatrix,bits:Vec>}impl<'a,T,M>WaveletMatrixPointAdd<'a,T,M>where T:Ord+Clone,M:AbelianGroup{pub fn update(&mut self,mut index:usize,value:M::T){debug_assert!(index)->M::T{let mut result=M::unit();self.wavelet_matrix.query_less_than(value,range,|d,range|{M::operate_assign(&mut result,&self.bits[self.wavelet_matrix.level(d)].fold_abelian(range.start,range.end));});result}pub fn fold_range(&self,values:Range,range:Range)->M::T{let lower=self.wavelet_matrix.compress.index_lower_bound(&values.start);let upper=self.wavelet_matrix.compress.index_lower_bound(&values.end);if lower>=upper{return M::unit();}let mut range=range;for d in(0..self.wavelet_matrix.bit_length).rev(){let level=self.wavelet_matrix.level(d);let start1=self.wavelet_matrix.bit_vectors[level].rank1(range.start);let end1=self.wavelet_matrix.bit_vectors[level].rank1(range.end);let start0=range.start-start1;let end0=range.end-end1;if((lower>>d)&1)==((upper>>d)&1){if((lower>>d)&1)==0{range=start0..end0;}else{range=self.wavelet_matrix.zeros[level]+start1..self.wavelet_matrix.zeros[level]+end1;}continue;}let zero_range=start0..end0;let one_range=self.wavelet_matrix.zeros[level]+start1..self.wavelet_matrix.zeros[level]+end1;let lower_sum=self.fold_lessthan_index(lower,zero_range.clone(),d);let upper_sum=self.fold_lessthan_index(upper,one_range,d);let zero_sum=self.bits[level].fold_abelian(zero_range.start,zero_range.end);let mut result=M::rinv_operate(&zero_sum,&lower_sum);M::operate_assign(&mut result,&upper_sum);return result;}M::unit()}fn fold_lessthan_index(&self,idx:usize,mut range:Range,bits:usize)->M::T{let mut result=M::unit();for d in(idx.trailing_zeros()as usize..bits).rev(){let level=self.wavelet_matrix.level(d);let start1=self.wavelet_matrix.bit_vectors[level].rank1(range.start);let end1=self.wavelet_matrix.bit_vectors[level].rank1(range.end);let start0=range.start-start1;let end0=range.end-end1;if((idx>>d)&1)!=0{M::operate_assign(&mut result,&self.bits[level].fold_abelian(start0,end0));range.start=self.wavelet_matrix.zeros[level]+start1;range.end=self.wavelet_matrix.zeros[level]+end1;}else{range.start=start0;range.end=end0;}}result}}#[derive(Debug,Clone)]pub struct WaveletMatrixFold<'a,T,M>where T:Ord+Clone,M:AbelianGroup{wavelet_matrix:&'a WaveletMatrix,prefix:Vec,offsets:Vec}impl<'a,T,M>WaveletMatrixFold<'a,T,M>where T:Ord+Clone,M:AbelianGroup{pub fn fold_lessthan(&self,val:T,range:Range)->M::T{self.fold_lessthan_with_count(val,range).1}pub fn fold_lessthan_with_count(&self,val:T,range:Range)->(usize,M::T){debug_assert!(range.end<=self.wavelet_matrix.len);let[result]=self.fold_lessthan_indices_with_count([self.wavelet_matrix.compress.index_lower_bound(&val)],[range],self.wavelet_matrix.bit_length);result}pub fn fold_range(&self,valrange:Range,range:Range)->M::T{self.fold_range_with_count(valrange,range).1}pub fn fold_range_with_count(&self,valrange:Range,mut range:Range)->(usize,M::T){debug_assert!(range.end<=self.wavelet_matrix.len);let lower=self.wavelet_matrix.compress.index_lower_bound(&valrange.start);let upper=self.wavelet_matrix.compress.index_lower_bound(&valrange.end);if lower>=upper{return(0,M::unit());}for d in(0..self.wavelet_matrix.bit_length).rev(){let level=self.wavelet_matrix.level(d);let start1=self.wavelet_matrix.bit_vectors[level].rank1(range.start);let end1=self.wavelet_matrix.bit_vectors[level].rank1(range.end);let start0=range.start-start1;let end0=range.end-end1;if((lower>>d)&1)==((upper>>d)&1){if((lower>>d)&1)==0{range=start0..end0;}else{range=self.wavelet_matrix.zeros[level]+start1..self.wavelet_matrix.zeros[level]+end1;}continue;}let zero_range=start0..end0;let one_range=self.wavelet_matrix.zeros[level]+start1..self.wavelet_matrix.zeros[level]+end1;let[(lower_count,lower_sum),(upper_count,upper_sum)]=self.fold_lessthan_indices_with_count([lower,upper],[zero_range.clone(),one_range],d);let zero_sum=self.range_sum(level,zero_range.clone());return(zero_range.len()-lower_count+upper_count,M::operate(&M::rinv_operate(&zero_sum,&lower_sum),&upper_sum));}(0,M::unit())}#[inline]fn range_sum(&self,level:usize,range:Range)->M::T{let offset=self.offsets[level];M::rinv_operate(&self.prefix[offset+range.end],&self.prefix[offset+range.start])}fn fold_lessthan_indices_with_count(&self,indices:[usize;N],mut ranges:[Range;N],bits:usize)->[(usize,M::T);N]{let mut results=std::array::from_fn(|_|(0,M::unit()));let last=indices.iter().map(|index|index.trailing_zeros()as usize).min().unwrap_or(bits);for d in(last..bits).rev(){let level=self.wavelet_matrix.level(d);for((&index,range),(count,sum))in indices.iter().zip(&mut ranges).zip(&mut results){let start1=self.wavelet_matrix.bit_vectors[level].rank1(range.start);let end1=self.wavelet_matrix.bit_vectors[level].rank1(range.end);let start0=range.start-start1;let end0=range.end-end1;if((index>>d)&1)!=0{*count+=end0-start0;M::operate_assign(sum,&self.range_sum(level,start0..end0));range.start=self.wavelet_matrix.zeros[level]+start1;range.end=self.wavelet_matrix.zeros[level]+end1;}else{range.start=start0;range.end=end0;}}}results}#[doc=" Folds the weights below each query's threshold, traversing the queries together."]pub fn fold_lessthan_batch(&self,queries:impl IntoIterator)>)->Vec{let queries:Vec<_>=queries.into_iter().collect();let mut result=Vec::with_capacity(queries.len());for queries in queries.chunks(16){let mut indices=[0;16];let mut ranges=std::array::from_fn(|_|0..0);for(i,(value,range))in queries.iter().enumerate(){assert!(range.start<=range.end&&range.end<=self.wavelet_matrix.len);indices[i]=self.wavelet_matrix.compress.index_lower_bound(value);ranges[i]=range.clone();}result.extend(self.fold_lessthan_indices_with_count(indices,ranges,self.wavelet_matrix.bit_length).into_iter().take(queries.len()).map(|(_,sum)|sum));}result}}} // codesnip-guard: algebra #[cfg_attr(any(),rustfmt::skip)]pub use self::magma::*;#[cfg_attr(any(),rustfmt::skip)]mod magma{#![doc=" algebraic traits"]#[doc=" binary operaion: $T \\circ T \\to T$"]pub trait Magma{#[doc=" type of operands: $T$"]type T:Clone;#[doc=" binary operaion: $\\circ$"]fn operate(x:&Self::T,y:&Self::T)->Self::T;fn reverse_operate(x:&Self::T,y:&Self::T)->Self::T{Self::operate(y,x)}fn operate_assign(x:&mut Self::T,y:&Self::T){*x=Self::operate(x,y);}}#[doc=" $\\forall a,\\forall b,\\forall c \\in T, (a \\circ b) \\circ c = a \\circ (b \\circ c)$"]pub trait Associative:Magma{}#[doc=" associative binary operation"]pub trait SemiGroup:Magma+Associative{}implSemiGroup for S where S:Magma+Associative{}#[doc=" $\\exists e \\in T, \\forall a \\in T, e \\circ a = a \\circ e = a$"]pub trait Unital:Magma{#[doc=" identity element: $e$"]fn unit()->Self::T;fn is_unit(x:&Self::T)->bool where Self::T:PartialEq{x==&Self::unit()}fn set_unit(x:&mut Self::T){*x=Self::unit();}}pub trait ExpBits{type Iter:Iterator;fn bits(self)->Self::Iter;}pub trait SignedExpBits{type T:ExpBits;fn neg_and_bits(self)->(bool,Self::T);}pub struct Bits{n:T}macro_rules!impl_exp_bits_for_uint{($($t:ty)*)=>{$(impl Iterator for Bits<$t>{type Item=bool;fn next(&mut self)->Option{if self.n==0{None}else{let bit=(self.n&1)==1;self.n>>=1;Some(bit)}}}impl ExpBits for$t{type Iter=Bits<$t>;fn bits(self)->Self::Iter{Bits{n:self}}}impl SignedExpBits for$t{type T=$t;fn neg_and_bits(self)->(bool,Self::T){(false,self)}})*};}impl_exp_bits_for_uint!(u8 u16 u32 u64 u128 usize);macro_rules!impl_signed_exp_bits_for_sint{($($s:ty,$u:ty;)*)=>{$(impl SignedExpBits for$s{type T=$u;fn neg_and_bits(self)->(bool,Self::T){(self<0,self.unsigned_abs())}})*};}impl_signed_exp_bits_for_sint!(i8,u8;i16,u16;i32,u32;i64,u64;i128,u128;isize,usize;);#[doc=" associative binary operation and an identity element"]pub trait Monoid:SemiGroup+Unital{#[doc=" binary exponentiation: $x^n = x\\circ\\ddots\\circ x$"]fn pow(mut x:Self::T,exp:E)->Self::T where E:ExpBits{let mut bits=exp.bits();let Some(first)=bits.next()else{return Self::unit();};let mut res=if first{x.clone()}else{Self::unit()};for bit in bits{x=Self::operate(&x,&x);if bit{res=Self::operate(&res,&x);}}res}fn fold(iter:I)->Self::T where I:IntoIterator{let mut iter=iter.into_iter();if let Some(item)=iter.next(){iter.fold(item,|acc,x|Self::operate(&acc,&x))}else{Self::unit()}}}implMonoid for M where M:SemiGroup+Unital{}#[doc=" $\\exists e \\in T, \\forall a \\in T, \\exists b,c \\in T, b \\circ a = a \\circ c = e$"]pub trait Invertible:Magma+Unital{#[doc=" $a$ where $a \\circ x = e$"]fn inverse(x:&Self::T)->Self::T;fn rinv_operate(x:&Self::T,y:&Self::T)->Self::T{Self::operate(x,&Self::inverse(y))}fn rinv_operate_assign(x:&mut Self::T,y:&Self::T){*x=Self::rinv_operate(x,y);}}#[doc=" associative binary operation and an identity element and inverse elements"]pub trait Group:Monoid+Invertible{fn signed_pow(x:Self::T,exp:E)->Self::T where E:SignedExpBits{let(neg,exp)=E::neg_and_bits(exp);let res=Self::pow(x,exp);if neg{Self::inverse(&res)}else{res}}}implGroup for G where G:Monoid+Invertible{}#[doc=" $\\forall a,\\forall b \\in T, a \\circ b = b \\circ a$"]pub trait Commutative:Magma{}#[doc=" commutative monoid"]pub trait AbelianMonoid:Monoid+Commutative{}implAbelianMonoid for M where M:Monoid+Commutative{}#[doc=" commutative group"]pub trait AbelianGroup:Group+Commutative{}implAbelianGroup for G where G:Group+Commutative{}#[doc=" $\\forall a \\in T, a \\circ a = a$"]pub trait Idempotent:Magma{}#[doc=" idempotent monoid"]pub trait IdempotentMonoid:Monoid+Idempotent{}implIdempotentMonoid for M where M:Monoid+Idempotent{}#[macro_export]macro_rules!monoid_fold{($m:ty)=>{<$m as Unital>::unit()};($m:ty,)=>{<$m as Unital>::unit()};($m:ty,$f:expr)=>{$f};($m:ty,$f:expr,$($ff:expr),*)=>{<$m as Magma>::operate(&($f),&monoid_fold!($m,$($ff),*))};}#[macro_export]macro_rules!define_monoid{($Name:ident,$t:ty,|$x:ident,$y:ident|$op:expr,$unit:expr)=>{struct$Name;impl Magma for$Name{type T=$t;fn operate($x:&Self::T,$y:&Self::T)->Self::T{$op}}impl Unital for$Name{fn unit()->Self::T{$unit}}impl Associative for$Name{}};}} // codesnip-guard: avx_helper #[cfg_attr(any(),rustfmt::skip)]pub use self::avx_helper::{SimdBackend,avx512_enabled,avx512_supported,disable_avx512,enable_avx512,simd_backend};#[cfg_attr(any(),rustfmt::skip)]mod avx_helper{use std::sync::atomic::{AtomicBool,Ordering};#[derive(Copy,Clone,Debug,Eq,PartialEq)]pub enum SimdBackend{Scalar,Avx2,Avx512}static AVX512_ENABLED:AtomicBool=AtomicBool::new(true);#[inline]pub fn disable_avx512(){AVX512_ENABLED.store(false,Ordering::Relaxed);}#[inline]pub fn enable_avx512(){AVX512_ENABLED.store(true,Ordering::Relaxed);}#[inline]pub fn avx512_enabled()->bool{AVX512_ENABLED.load(Ordering::Relaxed)}#[inline]pub fn avx512_supported()->bool{#[cfg(any(target_arch="x86",target_arch="x86_64"))]return is_x86_feature_detected!("avx512f")&&is_x86_feature_detected!("avx512dq")&&is_x86_feature_detected!("avx512cd")&&is_x86_feature_detected!("avx512bw")&&is_x86_feature_detected!("avx512vl");#[cfg(not(any(target_arch="x86",target_arch="x86_64")))]false}#[inline]pub fn simd_backend()->SimdBackend{#[cfg(any(target_arch="x86",target_arch="x86_64"))]{if avx512_enabled()&&avx512_supported(){return SimdBackend::Avx512;}if is_x86_feature_detected!("avx2"){return SimdBackend::Avx2;}}SimdBackend::Scalar}#[macro_export]macro_rules!avx_helper{(@dispatch$backend:path,$kind:ident;$avx512:expr,$avx2:expr,$scalar:expr)=>{{#[cfg(target_arch="x86_64")]{match$backend(){$kind::Avx512=>$avx512,$kind::Avx2=>$avx2,$kind::Scalar=>$scalar,}}#[cfg(not(target_arch="x86_64"))]$scalar}};(@dispatch_avx2_fma$avx2:expr,$scalar:expr)=>{{#[cfg(target_arch="x86_64")]{if is_x86_feature_detected!("avx2")&&is_x86_feature_detected!("fma"){$avx2}else{$scalar}}#[cfg(not(target_arch="x86_64"))]$scalar}};(@avx512$(#[$meta:meta])*$vis:vis fn$name:ident$(<$($T:ident),+>)?($($i:ident:$t:ty),*)->$ret:ty where[$($clauses:tt)*]$body:block)=>{$(#[$meta])*$vis fn$name$(<$($T)*>)?($($i:$t),*)->$ret where$($clauses)*{if$crate::avx512_supported(){$crate::avx_helper!(@def_avx512 fn avx512$(<$($T)*>)?($($i:$t),*)->$ret where[$($clauses)*]$body);unsafe{avx512$(::<$($T),*>)?($($i),*)}}else if is_x86_feature_detected!("avx2"){$crate::avx_helper!(@def_avx2 fn avx2$(<$($T)*>)?($($i:$t),*)->$ret where[$($clauses)*]$body);unsafe{avx2$(::<$($T),*>)?($($i),*)}}else{$body}}};(@avx2$(#[$meta:meta])*$vis:vis fn$name:ident$(<$($T:ident),+>)?($($i:ident:$t:ty),*)->$ret:ty where[$($clauses:tt)*]$body:block)=>{$(#[$meta])*$vis fn$name$(<$($T)*>)?($($i:$t),*)->$ret where$($clauses)*{if is_x86_feature_detected!("avx2"){$crate::avx_helper!(@def_avx2 fn avx2$(<$($T)*>)?($($i:$t),*)->$ret where[$($clauses)*]$body);unsafe{avx2$(::<$($T),*>)?($($i),*)}}else{$body}}};(@def_avx512 fn$name:ident$(<$($T:ident),+>)?($($args:tt)*)->$ret:ty where[$($clauses:tt)*]$body:block)=>{#[target_feature(enable="avx512f,avx512dq,avx512cd,avx512bw,avx512vl")]unsafe fn$name$(<$($T)*>)?($($args)*)->$ret where$($clauses)*$body};(@def_avx2 fn$name:ident$(<$($T:ident),+>)?($($args:tt)*)->$ret:ty where[$($clauses:tt)*]$body:block)=>{#[target_feature(enable="avx2")]unsafe fn$name$(<$($T)*>)?($($args)*)->$ret where$($clauses)*$body};(@$tag:ident$(#[$meta:meta])*$vis:vis fn$name:ident$(<$($T:ident),+>)?($($args:tt)*)->$ret:ty$body:block)=>{$crate::avx_helper!(@$tag$(#[$meta])*$vis fn$name$(<$($T)*>)?($($args)*)->$ret where[]$body);};(@$tag:ident$(#[$meta:meta])*$vis:vis fn$name:ident$(<$($T:ident),+>)?($($args:tt)*)$($t:tt)*)=>{$crate::avx_helper!(@$tag$(#[$meta])*$vis fn$name$(<$($T)*>)?($($args)*)->()$($t)*);};($($t:tt)*)=>{::std::compile_error!($($t)*);}}} // codesnip-guard: bounded #[cfg_attr(any(),rustfmt::skip)]pub use self::bounded::Bounded;#[cfg_attr(any(),rustfmt::skip)]mod bounded{#[doc=" Trait for max/min bounds"]pub trait Bounded:Sized+PartialOrd{fn maximum()->Self;fn minimum()->Self;fn is_maximum(&self)->bool{self==&Self::maximum()}fn is_minimum(&self)->bool{self==&Self::minimum()}fn set_maximum(&mut self){*self=Self::maximum()}fn set_minimum(&mut self){*self=Self::minimum()}}macro_rules!impl_bounded_num{($($t:ident)*)=>{$(impl Bounded for$t{fn maximum()->Self{$t::MAX}fn minimum()->Self{$t::MIN}})*};}impl_bounded_num!(u8 u16 u32 u64 u128 usize i8 i16 i32 i64 i128 isize f32 f64);macro_rules!impl_bounded_tuple{(@impl$($T:ident)*)=>{impl<$($T:Bounded),*>Bounded for($($T,)*){fn maximum()->Self{($(<$T as Bounded>::maximum(),)*)}fn minimum()->Self{($(<$T as Bounded>::minimum(),)*)}}};(@inner$($T:ident)*,)=>{impl_bounded_tuple!(@impl$($T)*);};(@inner$($T:ident)*,$U:ident$($Rest:ident)*)=>{impl_bounded_tuple!(@impl$($T)*);impl_bounded_tuple!(@inner$($T)*$U,$($Rest)*);};($T:ident$($Rest:ident)*)=>{impl_bounded_tuple!(@inner$T,$($Rest)*);};}impl_bounded_tuple!(A B C D E F G H I J);impl Bounded for(){fn maximum()->Self{}fn minimum()->Self{}}impl Bounded for bool{fn maximum()->Self{true}fn minimum()->Self{false}}implBounded for Optionwhere T:Bounded{fn maximum()->Self{Some(::maximum())}fn minimum()->Self{None}}implBounded for std::cmp::Reversewhere T:Bounded{fn maximum()->Self{std::cmp::Reverse(::minimum())}fn minimum()->Self{std::cmp::Reverse(::maximum())}}} // codesnip-guard: compress #[cfg_attr(any(),rustfmt::skip)]pub use self::compress::{Compressor,HashCompress,OrderedCompressor,VecCompress};#[cfg_attr(any(),rustfmt::skip)]mod compress{use std::{collections::HashMap,fmt::{self,Debug},hash::Hash,iter::FromIterator};pub trait Compressorwhere Self:FromIterator,T:Ord{fn index_exact(&self,index:&T)->Option;fn size(&self)->usize;}pub trait OrderedCompressor:Compressorwhere T:Ord{fn index_lower_bound(&self,index:&T)->usize;}#[derive(Debug,Clone)]pub struct VecCompress{data:Vec}implVecCompress{pub fn from_sorted_unique(data:Vec)->Self{Self{data}}pub fn values(&self)->&[T]{&self.data}}implFromIteratorfor VecCompresswhere T:Ord{fn from_iter(iter:I)->Self where I:IntoIterator{let mut data:Vec<_>=iter.into_iter().collect();data.sort_unstable();data.dedup();Self{data}}}implCompressorfor VecCompresswhere T:Ord{fn index_exact(&self,index:&T)->Option{self.data.binary_search(index).ok()}fn size(&self)->usize{self.data.len()}}implOrderedCompressorfor VecCompresswhere T:Ord{fn index_lower_bound(&self,index:&T)->usize{self.data.partition_point(|x|x{data:HashMap}implDebug for HashCompresswhere T:Debug+Eq+Hash{fn fmt(&self,f:&mut fmt::Formatter<'_>)->fmt::Result{f.debug_struct("HashCompress").field("data",&self.data).finish()}}implFromIteratorfor HashCompresswhere T:Ord+Hash{fn from_iter(iter:I)->Self where I:IntoIterator{let mut data:Vec<_>=iter.into_iter().collect();data.sort_unstable();data.dedup();let data=data.into_iter().enumerate().map(|(i,t)|(t,i)).collect();Self{data}}}implCompressorfor HashCompresswhere T:Ord+Hash{fn index_exact(&self,index:&T)->Option{self.data.get(index).copied()}fn size(&self)->usize{self.data.len()}}} // codesnip-guard: discrete_steps #[cfg_attr(any(),rustfmt::skip)]pub use self::discrete_steps::{DiscreteSteps,RangeBoundsExt};#[cfg_attr(any(),rustfmt::skip)]mod discrete_steps{use super::Bounded;use std::{convert::TryFrom,ops::{Bound,Range,RangeBounds,RangeInclusive}};pub trait DiscreteSteps:Clone{fn delta()->Delta;fn steps_between(start:&Self,end:&Self)->Option;fn forward_checked(start:Self,delta:Delta)->Option;fn backward_checked(start:Self,delta:Delta)->Option;fn forward(start:Self,delta:Delta)->Self{Self::forward_checked(start,delta).expect("overflow in `DiscreteSteps::forward`")}fn backward(start:Self,delta:Delta)->Self{Self::backward_checked(start,delta).expect("overflow in `DiscreteSteps::backward`")}fn forward_delta_checked(start:Self)->Option{Self::forward_checked(start,Self::delta())}fn backward_delta_checked(start:Self)->Option{Self::backward_checked(start,Self::delta())}fn forward_delta(start:Self)->Self{Self::forward(start,Self::delta())}fn backward_delta(start:Self)->Self{Self::backward(start,Self::delta())}}macro_rules!impl_discrete_steps_integer{(@common$u_source:ident)=>{fn delta()->$u_source{1}fn forward(start:Self,delta:$u_source)->Self{assert!(Self::forward_checked(start,delta).is_some(),"attempt to add with overflow");start.wrapping_add(delta as Self)}fn backward(start:Self,delta:$u_source)->Self{assert!(Self::backward_checked(start,delta).is_some(),"attempt to subtract with overflow");start.wrapping_sub(delta as Self)}};($u_source:ident$i_source:ident;$($u_narrower:ident$i_narrower:ident),*;$($u_wider:ident$i_wider:ident),*)=>{$(impl DiscreteSteps<$u_source>for$u_narrower{impl_discrete_steps_integer!(@common$u_source);fn steps_between(start:&Self,end:&Self)->Option<$u_source>{if*start<=*end{$u_source::try_from(*end-*start).ok()}else{None}}fn forward_checked(start:Self,delta:$u_source)->Option{Self::try_from(delta).ok().and_then(|delta|start.checked_add(delta))}fn backward_checked(start:Self,delta:$u_source)->Option{Self::try_from(delta).ok().and_then(|delta|start.checked_sub(delta))}}impl DiscreteSteps<$u_source>for$i_narrower{impl_discrete_steps_integer!(@common$u_source);fn steps_between(start:&Self,end:&Self)->Option<$u_source>{if*start<=*end{$u_source::try_from(end.abs_diff(*start)).ok()}else{None}}fn forward_checked(start:Self,delta:$u_source)->Option{$u_narrower::try_from(delta).ok().and_then(|delta|{let wrapped=start.wrapping_add(delta as Self);if wrapped>=start{Some(wrapped)}else{None}})}fn backward_checked(start:Self,delta:$u_source)->Option{$u_narrower::try_from(delta).ok().and_then(|delta|{let wrapped=start.wrapping_sub(delta as Self);if wrapped<=start{Some(wrapped)}else{None}})}})*$(impl DiscreteSteps<$u_source>for$u_wider{impl_discrete_steps_integer!(@common$u_source);fn steps_between(start:&Self,end:&Self)->Option<$u_source>{if*start<=*end{$u_source::try_from(*end-*start).ok()}else{None}}fn forward_checked(start:Self,delta:$u_source)->Option{start.checked_add(delta as Self)}fn backward_checked(start:Self,delta:$u_source)->Option{start.checked_sub(delta as Self)}}impl DiscreteSteps<$u_source>for$i_wider{impl_discrete_steps_integer!(@common$u_source);fn steps_between(start:&Self,end:&Self)->Option<$u_source>{if*start<=*end{end.checked_sub(*start).and_then(|result|$u_source::try_from(result).ok())}else{None}}fn forward_checked(start:Self,delta:$u_source)->Option{start.checked_add(delta as Self)}fn backward_checked(start:Self,delta:$u_source)->Option{start.checked_sub(delta as Self)}})*};}impl_discrete_steps_integer!(u16 i16;u8 i8,u16 i16,usize isize;u32 i32,u64 i64,u128 i128);impl_discrete_steps_integer!(u32 i32;u8 i8,u16 i16,u32 i32,usize isize;u64 i64,u128 i128);impl_discrete_steps_integer!(u64 i64;u8 i8,u16 i16,u32 i32,u64 i64,usize isize;u128 i128);impl_discrete_steps_integer!(u128 i128;u8 i8,u16 i16,u32 i32,u64 i64,u128 i128,usize isize;);impl_discrete_steps_integer!(usize isize;u8 i8,u16 i16,u32 i32,u64 i64,usize isize;u128 i128);pub trait RangeBoundsExt{fn start_bound_included_checked(&self)->Option;fn start_bound_excluded_checked(&self)->Option;fn end_bound_included_checked(&self)->Option;fn end_bound_excluded_checked(&self)->Option;fn start_bound_included(&self)->T;fn start_bound_excluded(&self)->T;fn end_bound_included(&self)->T;fn end_bound_excluded(&self)->T;fn start_bound_included_bounded(&self,lb:T)->Optionwhere T:Ord;fn start_bound_excluded_bounded(&self,lb:T)->Optionwhere T:Ord;fn end_bound_included_bounded(&self,ub:T)->Optionwhere T:Ord;fn end_bound_excluded_bounded(&self,ub:T)->Optionwhere T:Ord;fn to_range_checked(&self)->Option>{match(self.start_bound_included_checked(),self.end_bound_excluded_checked()){(Some(start),Some(end))=>Some(start..end),_=>None,}}fn to_range(&self)->Range{self.start_bound_included()..self.end_bound_excluded()}fn to_range_bounded(&self,min:T,max:T)->Option>where T:Ord{Some(self.start_bound_included_bounded(min)?..self.end_bound_excluded_bounded(max)?)}fn to_range_inclusive_checked(&self)->Option>{match(self.start_bound_included_checked(),self.end_bound_included_checked()){(Some(start),Some(end))=>Some(start..=end),_=>None,}}fn to_range_inclusive(&self)->RangeInclusive{self.start_bound_included()..=self.end_bound_included()}fn to_range_inclusive_bounded(&self,min:T,max:T)->Option>where T:Ord{Some(self.start_bound_included_bounded(min)?..=self.end_bound_included_bounded(max)?)}}macro_rules!impl_range_bounds_ext{($($source:ident=>$($target:ident)+);*$(;)?)=>{$($(implRangeBoundsExt<$target>for R where R:RangeBounds<$target>,{fn start_bound_included_checked(&self)->Option<$target>{match self.start_bound(){Bound::Included(x)=>Some(*x),Bound::Excluded(x)=>DiscreteSteps::<$source>::forward_delta_checked(*x),Bound::Unbounded=>Some(Bounded::minimum()),}}fn start_bound_excluded_checked(&self)->Option<$target>{match self.start_bound(){Bound::Included(x)=>DiscreteSteps::<$source>::backward_delta_checked(*x),Bound::Excluded(x)=>Some(*x),Bound::Unbounded=>None,}}fn end_bound_included_checked(&self)->Option<$target>{match self.end_bound(){Bound::Included(x)=>Some(*x),Bound::Excluded(x)=>DiscreteSteps::<$source>::backward_delta_checked(*x),Bound::Unbounded=>Some(Bounded::maximum()),}}fn end_bound_excluded_checked(&self)->Option<$target>{match self.end_bound(){Bound::Included(x)=>DiscreteSteps::<$source>::forward_delta_checked(*x),Bound::Excluded(x)=>Some(*x),Bound::Unbounded=>None,}}fn start_bound_included(&self)->$target{match self.start_bound(){Bound::Included(x)=>*x,Bound::Excluded(x)=>DiscreteSteps::<$source>::forward_delta(*x),Bound::Unbounded=>Bounded::minimum(),}}fn start_bound_excluded(&self)->$target{match self.start_bound(){Bound::Included(x)=>DiscreteSteps::<$source>::backward_delta(*x),Bound::Excluded(x)=>*x,Bound::Unbounded=>DiscreteSteps::<$source>::backward_delta(Bounded::minimum()),}}fn end_bound_included(&self)->$target{match self.end_bound(){Bound::Included(x)=>*x,Bound::Excluded(x)=>DiscreteSteps::<$source>::backward_delta(*x),Bound::Unbounded=>Bounded::maximum(),}}fn end_bound_excluded(&self)->$target{match self.end_bound(){Bound::Included(x)=>DiscreteSteps::<$source>::forward_delta(*x),Bound::Excluded(x)=>*x,Bound::Unbounded=>DiscreteSteps::<$source>::forward_delta(Bounded::maximum()),}}fn start_bound_included_bounded(&self,lb:$target)->Option<$target>where$target:Ord{match self.start_bound(){Bound::Included(x)=>Some(*x).filter(|&x|lb<=x),Bound::Excluded(x)=>DiscreteSteps::<$source>::forward_delta_checked(*x).filter(|&x|lb<=x),Bound::Unbounded=>Some(lb),}}fn start_bound_excluded_bounded(&self,lb:$target)->Option<$target>where$target:Ord{match self.start_bound(){Bound::Included(x)=>DiscreteSteps::<$source>::backward_delta_checked(*x).filter(|&x|lb<=x),Bound::Excluded(x)=>Some(*x).filter(|&x|lb<=x),Bound::Unbounded=>Some(lb),}}fn end_bound_included_bounded(&self,ub:$target)->Option<$target>where$target:Ord{match self.end_bound(){Bound::Included(x)=>Some(*x).filter(|&x|x<=ub),Bound::Excluded(x)=>DiscreteSteps::<$source>::backward_delta_checked(*x).filter(|&x|x<=ub),Bound::Unbounded=>Some(ub),}}fn end_bound_excluded_bounded(&self,ub:$target)->Option<$target>where$target:Ord{match self.end_bound(){Bound::Included(x)=>DiscreteSteps::<$source>::forward_delta_checked(*x).filter(|&x|x<=ub),Bound::Excluded(x)=>Some(*x).filter(|&x|x<=ub),Bound::Unbounded=>Some(ub),}}})+)*};}impl_range_bounds_ext!(u16=>u8 i8 u16 i16;u32=>u32 i32;u64=>u64 i64;u128=>u128 i128;usize=>isize usize;);} // codesnip-guard: zero_one #[cfg_attr(any(),rustfmt::skip)]pub use self::zero_one::{One,Zero};#[cfg_attr(any(),rustfmt::skip)]mod zero_one{pub trait Zero:Sized{fn zero()->Self;#[inline]fn is_zero(&self)->bool where Self:PartialEq{self==&Self::zero()}#[inline]fn set_zero(&mut self){*self=Self::zero();}}pub trait One:Sized{fn one()->Self;#[inline]fn is_one(&self)->bool where Self:PartialEq{self==&Self::one()}#[inline]fn set_one(&mut self){*self=Self::one();}}macro_rules!impl_zero_one{($({$(<$T:ident:$Bound:ident>)?$Trait:ident$method:ident$($t:ty)*,$e:expr})*)=>{$(impl_zero_one!(@impl[$(<$T:$Bound>)?]$Trait$method[$($t)*],$e);)*};(@impl[<$T:ident:$Bound:ident>]$Trait:ident$method:ident[$($t:ty)*],$e:expr)=>{$(impl<$T:$Bound>$Trait for$t{fn$method()->Self{$e}})*};(@impl[]$Trait:ident$method:ident[$($t:ty)*],$e:expr)=>{$(impl$Trait for$t{fn$method()->Self{$e}})*};}impl_zero_one!({Zero zero u8 u16 u32 u64 usize i8 i16 i32 i64 isize u128 i128,0}{Zero zero f32 f64,0.}{One one u8 u16 u32 u64 usize i8 i16 i32 i64 isize u128 i128,1}{One one f32 f64,1.}{Zero zero std::num::Wrapping,Self(T::zero())}{One one std::num::Wrapping,Self(T::one())});}