// codesnip-guard: RangeLongestIncreasingSubsequence
#[cfg_attr(any(),rustfmt::skip)]pub use self::range_longest_increasing_subsequence::{RangeLongestIncreasingSubsequence,WeightedRangeLongestIncreasingSubsequence};#[cfg_attr(any(),rustfmt::skip)]mod range_longest_increasing_subsequence{use super::{WaveletMatrix,Zero};use std::{cmp::Reverse,collections::BinaryHeap,ops::Add};#[derive(Debug,Clone)]pub struct RangeLongestIncreasingSubsequence{len:usize,positions:WaveletMatrix<usize>}impl RangeLongestIncreasingSubsequence{pub fn new<T>(sequence:&[T])->Self where T:Ord{let mut order:Vec<_>=(0..sequence.len()).collect();order.sort_unstable_by_key(|&i|(&sequence[i],Reverse(i)));let mut permutation=vec![0;sequence.len()];for(rank,i)in order.into_iter().enumerate(){permutation[i]=rank;}Self::from_permutation(&permutation)}pub fn from_permutation(permutation:&[usize])->Self{let len=permutation.len();assert!(len<=(u32::MAX/2)as usize);let mut used=vec![false;len];for&value in permutation{assert!(!used[value]);used[value]=true;}drop(used);let positions=if len<=1<<19{Self::seaweed_blocks(permutation)}else{Self::seaweed_doubling(permutation,1<<16)};Self{len,positions:WaveletMatrix::new(positions)}}pub fn longest_length(&self,left:usize,right:usize)->usize{assert!(left<=right&&right<=self.len);right-left-self.positions.rank_lessthan(right,left..self.len)}#[doc=" Answers follow the input order; intervals are half-open."]pub fn longest_length_batch(&self,queries:impl IntoIterator<Item=(usize,usize)>)->Vec<usize>{let queries=queries.into_iter();let mut result=Vec::with_capacity(queries.size_hint().0);let counts=self.positions.rank_lessthan_batch(queries.map(|(left,right)|{assert!(left<=right&&right<=self.len);result.push(right-left);(right,left..self.len)}));for(length,count)in result.iter_mut().zip(counts){*length-=count;}result}fn seaweed_doubling(permutation:&[usize],block_size:usize)->Vec<usize>{let n=permutation.len();let mut indices=vec![0u32;n];let mut row=vec![0u32;n];let mut column=vec![0u32;n];let mut links=vec![0u32;n];let mut buffer=vec![0u32;n];let mut row_order=vec![0u32;n];let mut column_order=vec![0u32;n];let mut marked=vec![false;n];let mut local=Vec::with_capacity(block_size);let mut offsets:Vec<_>=(0..n).step_by(block_size).collect();for(i,&value)in permutation.iter().enumerate(){let at=&mut offsets[value/block_size];indices[*at]=i as u32;buffer[*at]=(value%block_size)as u32;*at+=1;}for start in(0..n).step_by(block_size){let end=(start+block_size).min(n);local.clear();local.extend(buffer[start..end].iter().map(|&x|x as usize));for(i,j)in Self::seaweed_blocks(&local).into_iter().enumerate(){if j<end-start{links[start+j]=(start+i+1)as u32;}}}let mut width=block_size;while width<n{for mid in(width..n).step_by(2*width){let start=mid-width;let end=(mid+width).min(n);if indices[mid-1]<indices[mid]{continue;}let(mut i,mut j)=(start,mid);for(position,entry)in buffer.iter_mut().enumerate().take(end).skip(start){let index=if j==end||(i<mid&&indices[i]<indices[j]){let x=i;i+=1;x}else{let x=j;j+=1;x};row_order[index]=position as u32;*entry=indices[index];}indices[start..end].copy_from_slice(&buffer[start..end]);column[start..end].fill(end as u32);for i in start..mid{row[row_order[i]as usize]=if links[i]==0{0}else{row_order[(links[i]-1)as usize]+1};column[row_order[i]as usize]=row_order[i];}for i in mid..end{row[row_order[i]as usize]=row_order[i]+1;if links[i]!=0{column[row_order[(links[i]-1)as usize]as usize]=row_order[i];}}for i in start..end{row_order[i]=i as u32;column_order[i]=i as u32;}let mut step=1;while step<end-start{for middle in(start+step..end).step_by(2*step){let left=middle-step;let right=(middle+step).min(end);let(mut i,mut j)=(left,middle);for entry in&mut buffer[left..right]{let index=if j==right||(i<middle&&row[row_order[i]as usize]>row[row_order[j]as usize]){let x=i;i+=1;x}else{let x=j;j+=1;x};*entry=row_order[index];}row_order[left..right].copy_from_slice(&buffer[left..right]);marked[left..middle].fill(true);marked[middle..right].fill(false);let mut x=left;let(mut i,mut j)=(left,middle);for entry in&mut buffer[left..right]{if j==right||(i<middle&&column[column_order[i]as usize]<column[column_order[j]as usize]){if marked[column_order[i]as usize]{links[column_order[i]as usize]=column[column_order[i]as usize];marked[column_order[i]as usize]=false;*entry=column_order[i];i+=1;}else{while!marked[row_order[x]as usize]{x+=1;}links[row_order[x]as usize]=column[column_order[i]as usize];marked[row_order[x]as usize]=false;*entry=row_order[x];i+=1;}}else if row[column_order[j]as usize]>=row[row_order[x]as usize]{links[column_order[j]as usize]=column[column_order[j]as usize];*entry=column_order[j];j+=1;}else{marked[column_order[j]as usize]=true;while!marked[row_order[x]as usize]{x+=1;}links[row_order[x]as usize]=column[column_order[j]as usize];marked[row_order[x]as usize]=false;*entry=row_order[x];j+=1;}}column_order[left..right].copy_from_slice(&buffer[left..right]);column[left..right].copy_from_slice(&links[left..right]);}step*=2;}links[start..end].fill(0);for i in start..end{if column[i]<end as u32{links[column[i]as usize]=row[i];}}}width*=2;}drop((indices,row,column,buffer,row_order,column_order,marked));let mut positions=vec![n;n];for(i,x)in links.into_iter().enumerate(){if x!=0{positions[(x-1)as usize]=i;}}positions}fn seaweed_blocks(permutation:&[usize])->Vec<usize>{let n=permutation.len();let width=(n.isqrt()*3/4).max(1);let blocks=n.div_ceil(width);let mut boundary=vec![0u32;n];let mut positions=vec![n;n];let mut heaps:Vec<BinaryHeap<u32>>=(0..blocks).map(|_|BinaryHeap::with_capacity(width)).collect();let mut dirty=vec![false;blocks];let mut maxima=vec![0;blocks];let mut ends=vec![0;blocks];let mut pending=BinaryHeap::with_capacity(width);let mut active_blocks=0;for(i,&value)in permutation.iter().enumerate(){let block=value/width;active_blocks=active_blocks.max(block+1);let start=block*width;let end=ends[block].max(value+1);let append_only=!dirty[block]&&value>=ends[block];ends[block]=end;if dirty[block]{dirty[block]=false;pending.clear();pending.extend(heaps[block].iter().filter(|&&x|x&1!=0).map(|&x|Reverse(x>>1)));for entry in&mut boundary[start..end]{if pending.peek().is_some_and(|smallest|*entry>smallest.0){std::mem::swap(entry,&mut pending.peek_mut().unwrap().0);}}}boundary[value]=(i+1)as u32;let mut displaced=0;for entry in&mut boundary[value+1..end]{if*entry>displaced{std::mem::swap(entry,&mut displaced);}}for(offset,maximum)in maxima[block+1..active_blocks].iter_mut().enumerate(){if*maximum>displaced{let index=block+1+offset;let mut largest=heaps[index].peek_mut().unwrap();let next=*maximum;*largest=displaced<<1|1;displaced=next;drop(largest);*maximum=heaps[index].peek().copied().unwrap()>>1;dirty[index]=true;}}if displaced!=0{positions[displaced as usize-1]=i;}let heap=&mut heaps[block];if append_only{heap.push(((i+1)as u32)<<1);}else{heap.clear();heap.extend(boundary[start..end].iter().filter(|&&x|x!=0).map(|&x|x<<1));}maxima[block]=heap.peek().copied().unwrap()>>1;}positions}}#[doc=" Weights must be positive. Addition must be exact, associative, commutative,"]#[doc=" and strictly preserve order, with `Zero::zero()` as its identity."]#[derive(Debug,Clone)]pub struct WeightedRangeLongestIncreasingSubsequence<W>{len:usize,nodes:Vec<[(u32,u32);2]>,subsequences:Vec<SubsequenceNode<W>>,roots:Vec<(u32,u32)>}impl<W>WeightedRangeLongestIncreasingSubsequence<W>{pub fn new<T>(sequence:&[T],weights:&[W])->Self where T:Ord,W:Clone+Ord+Zero+Add<Output=W>{let len=sequence.len();assert_eq!(len,weights.len());assert!(len<u32::MAX as usize);for weight in weights{assert!(weight>&W::zero());}let mut result=WeightedRangeLongestIncreasingSubsequenceBuilder{nodes:vec![vec![FunctionNode{left:(0,0),right:(0,0),first:W::zero(),last:W::zero(),}]],subsequences:vec![vec![SubsequenceNode{left:0,right:0,weight:W::zero(),}]],roots:Vec::with_capacity(len+1)};let mut order:Vec<_>=(0..len).collect();order.sort_unstable_by_key(|&i|(&sequence[i],Reverse(i)));let mut ranks=vec![0;len];for(rank,i)in order.into_iter().enumerate(){ranks[i]=rank;}let size=len.next_power_of_two();let mut outer=vec![(0,0);2*size];let mut ends=vec![0;2*size];result.roots.push((0,0));for(i,(&rank,weight))in ranks.iter().zip(weights).enumerate(){let mut position=size+rank;let mut predecessor=(0,0);let mut end=0;while position>1{if position&1!=0&&ends[position-1]>end{predecessor=result.merge(outer[position-1],predecessor,len);end=ends[position-1];}position/=2;}let suffix=result.push_subsequence(SubsequenceNode{left:0,right:i as u32,weight:weight.clone()});let mut position=size+rank;outer[position]=result.append_prefix(predecessor,len,i+1,suffix);ends[position]=i+1;while position>1{let right=position&1!=0;position/=2;ends[position]=i+1;outer[position]=if right{outer[2*position+1]}else{result.merge(outer[2*position],outer[2*position+1],len)};}result.roots.push(outer[1]);if(i+1).is_power_of_two(){result.compact(&mut outer);}}drop(outer);result.compact(&mut[]);let mut nodes=Vec::with_capacity(result.nodes.iter().map(Vec::len).sum());for chunk in result.nodes{nodes.extend(chunk.into_iter().map(|node|[node.left,node.right]));}let mut subsequences=Vec::with_capacity(result.subsequences.iter().map(Vec::len).sum());for chunk in result.subsequences{subsequences.extend(chunk);}Self{len,nodes,subsequences,roots:result.roots}}pub fn maximum_weight(&self,left:usize,right:usize)->W where W:Clone+Zero+Add<Output=W>{assert!(left<=right&&right<=self.len);if left==right{return W::zero();}let(mut tree,mut len,mut position)=(self.roots[right],self.len,left);let mut weight=self.subsequences[tree.1 as usize].weight.clone();while tree.0!=0{let node=self.nodes[tree.0 as usize];let mid=len/2;if position<mid{tree=node[0];len=mid;}else{tree=node[1];position-=mid;len-=mid;}weight=weight+self.subsequences[tree.1 as usize].weight.clone();}weight}#[doc=" Returns the original indices of one optimum in the half-open interval."]pub fn maximum_weight_subsequence(&self,left:usize,right:usize)->Vec<usize>{assert!(left<=right&&right<=self.len);if left==right{return Vec::new();}let(mut tree,mut len,mut position)=(self.roots[right],self.len,left);let mut suffixes=Vec::new();if tree.1!=0{suffixes.push(tree.1);}while tree.0!=0{let node=self.nodes[tree.0 as usize];let mid=len/2;if position<mid{tree=node[0];len=mid;}else{tree=node[1];position-=mid;len-=mid;}if tree.1!=0{suffixes.push(tree.1);}}let mut result=Vec::new();let mut stack=Vec::new();for suffix in suffixes.into_iter().rev(){stack.push(suffix);while let Some(index)=stack.pop(){let node=&self.subsequences[index as usize];if node.left==0{result.push(node.right as usize);}else{stack.push(node.right);stack.push(node.left);}}}result}}struct WeightedRangeLongestIncreasingSubsequenceBuilder<W>{nodes:Vec<Vec<FunctionNode<W>>>,subsequences:Vec<Vec<SubsequenceNode<W>>>,roots:Vec<(u32,u32)>}impl<W>WeightedRangeLongestIncreasingSubsequenceBuilder<W>{fn compact(&mut self,outer:&mut[(u32,u32)])where W:Clone{let node_len=(self.nodes.len()-1)*CHUNK_SIZE+self.nodes.last().unwrap().len();let subsequence_len=(self.subsequences.len()-1)*CHUNK_SIZE+self.subsequences.last().unwrap().len();let mut nodes=vec![0u32;node_len];let mut subsequences=vec![0u32;subsequence_len];for&(tree,suffix)in self.roots.iter().chain(outer.iter()){nodes[tree as usize]=1;subsequences[suffix as usize]=1;}for index in(1..node_len).rev(){if nodes[index]==0{continue;}let node=self.node(index as u32);nodes[node.left.0 as usize]=1;nodes[node.right.0 as usize]=1;subsequences[node.left.1 as usize]=1;subsequences[node.right.1 as usize]=1;}for index in(1..subsequence_len).rev(){if subsequences[index]==0{continue;}let node=self.subsequence(index as u32);if node.left!=0{subsequences[node.left as usize]=1;subsequences[node.right as usize]=1;}}nodes[0]=0;subsequences[0]=0;let mut count=1;for index in 1..subsequence_len{if subsequences[index]==0{continue;}let mut node=self.subsequence(index as u32).clone();if node.left!=0{node.left=subsequences[node.left as usize];node.right=subsequences[node.right as usize];}self.subsequences[count/CHUNK_SIZE][count%CHUNK_SIZE]=node;subsequences[index]=count as u32;count+=1;}self.subsequences.truncate(count.div_ceil(CHUNK_SIZE));self.subsequences.last_mut().unwrap().truncate((count-1)%CHUNK_SIZE+1);count=1;for index in 1..node_len{if nodes[index]==0{continue;}let mut node=self.node(index as u32).clone();node.left=(nodes[node.left.0 as usize],subsequences[node.left.1 as usize]);node.right=(nodes[node.right.0 as usize],subsequences[node.right.1 as usize]);self.nodes[count/CHUNK_SIZE][count%CHUNK_SIZE]=node;nodes[index]=count as u32;count+=1;}self.nodes.truncate(count.div_ceil(CHUNK_SIZE));self.nodes.last_mut().unwrap().truncate((count-1)%CHUNK_SIZE+1);for tree in self.roots.iter_mut().chain(outer.iter_mut()){*tree=(nodes[tree.0 as usize],subsequences[tree.1 as usize]);}}fn node(&self,index:u32)->&FunctionNode<W>{&self.nodes[index as usize/CHUNK_SIZE][index as usize%CHUNK_SIZE]}fn subsequence(&self,index:u32)->&SubsequenceNode<W>{&self.subsequences[index as usize/CHUNK_SIZE][index as usize%CHUNK_SIZE]}fn push_subsequence(&mut self,node:SubsequenceNode<W>)->u32{if self.subsequences.last().unwrap().len()==CHUNK_SIZE{self.subsequences.push(Vec::with_capacity(CHUNK_SIZE));}let index=(self.subsequences.len()-1)*CHUNK_SIZE+self.subsequences.last().unwrap().len();let index=index.try_into().expect("too many subsequence nodes");self.subsequences.last_mut().unwrap().push(node);index}fn concatenate(&mut self,left:u32,right:u32)->u32 where W:Clone+Add<Output=W>{if left==0{return right;}if right==0{return left;}self.push_subsequence(SubsequenceNode{left,right,weight:self.subsequence(left).weight.clone()+self.subsequence(right).weight.clone()})}fn join(&mut self,left:(u32,u32),right:(u32,u32))->(u32,u32)where W:Clone+Add<Output=W>{if left==right&&left.0==0{return left;}let first=self.node(left.0).first.clone()+self.subsequence(left.1).weight.clone();let last=self.node(right.0).last.clone()+self.subsequence(right.1).weight.clone();if self.nodes.last().unwrap().len()==CHUNK_SIZE{self.nodes.push(Vec::with_capacity(CHUNK_SIZE));}let index=(self.nodes.len()-1)*CHUNK_SIZE+self.nodes.last().unwrap().len();let index=index.try_into().expect("too many function nodes");self.nodes.last_mut().unwrap().push(FunctionNode{left,right,first,last});(index,0)}fn append(&mut self,tree:(u32,u32),suffix:u32)->(u32,u32)where W:Clone+Add<Output=W>{(tree.0,self.concatenate(tree.1,suffix))}fn append_prefix(&mut self,tree:(u32,u32),len:usize,end:usize,suffix:u32)->(u32,u32)where W:Clone+Add<Output=W>{if end==0{return tree;}if len<=end{return self.append(tree,suffix);}let node=self.node(tree.0);let(left,right)=(node.left,node.right);let left=self.append(left,tree.1);let right=self.append(right,tree.1);let mid=len/2;let left=self.append_prefix(left,mid,end,suffix);let right=self.append_prefix(right,len-mid,end.saturating_sub(mid),suffix);self.join(left,right)}fn merge(&mut self,low:(u32,u32),high:(u32,u32),len:usize)->(u32,u32)where W:Clone+Ord+Add<Output=W>{if high==(0,0)||low==high{return low;}if low==(0,0){return high;}let low_node=self.node(low.0);let high_node=self.node(high.0);let low_weight=&self.subsequence(low.1).weight;let high_weight=&self.subsequence(high.1).weight;if high_node.last.clone()+high_weight.clone()>low_node.last.clone()+low_weight.clone(){return high;}if high_node.first.clone()+high_weight.clone()<=low_node.first.clone()+low_weight.clone(){return low;}let mid=len/2;let high_wins=self.node(high_node.left.0).last.clone()+self.subsequence(high_node.left.1).weight.clone()+high_weight.clone()>self.node(low_node.left.0).last.clone()+self.subsequence(low_node.left.1).weight.clone()+low_weight.clone();let(low_left,low_right)=(low_node.left,low_node.right);let(high_left,high_right)=(high_node.left,high_node.right);if high_wins{let left=self.append(high_left,high.1);let low_right=self.append(low_right,low.1);let high_right=self.append(high_right,high.1);let right=self.merge(low_right,high_right,len-mid);self.join(left,right)}else{let low_left=self.append(low_left,low.1);let high_left=self.append(high_left,high.1);let left=self.merge(low_left,high_left,mid);let right=self.append(low_right,low.1);self.join(left,right)}}}const CHUNK_SIZE:usize=1<<12;#[derive(Debug,Clone)]struct FunctionNode<W>{left:(u32,u32),right:(u32,u32),first:W,last:W}#[derive(Debug,Clone)]struct SubsequenceNode<W>{left:u32,right:u32,weight:W}}
// 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 BinaryIndexedTree<M>where M:Monoid{n:usize,bit:Vec<M::T>}impl<M>Clone for BinaryIndexedTree<M>where M:Monoid{fn clone(&self)->Self{Self{n:self.n,bit:self.bit.clone()}}}impl<M>Debug for BinaryIndexedTree<M>where M:Monoid<T:Debug>{fn fmt(&self,f:&mut Formatter<'_>)->fmt::Result{f.debug_struct("BinaryIndexedTree").field("n",&self.n).field("bit",&self.bit).finish()}}impl<M>BinaryIndexedTree<M>where 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.n);let mut k=k+1;while k<=self.n{self.bit[k]=M::operate(&self.bit[k],&x);k+=k&(!k+1);}}#[inline]pub fn partition_point_acc<P>(&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}}impl<G:Group>BinaryIndexedTree<G>{#[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<usize>{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<usize>{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 k<WORD_SIZE{self&(1 as$t)<<k!=0}else{false}}fn rank1(&self,k:usize)->usize{const WORD_SIZE:usize=(0 as$t).count_zeros()as usize;if k<WORD_SIZE{(self&!(!(0 as$t)<<k)).count_ones()as usize}else{self.count_ones()as usize}}})*};}impl_rank_select_for_bits!(u8 u16 u32 u64 usize i8 i16 i32 i64 isize u128 i128);#[inline]fn select_word_scalar(mut word:u64,mut rank:usize)->usize{let count=word.count_ones()as usize;debug_assert!(rank<count);if rank<4{for _ in 0..rank{word&=word-1;}return word.trailing_zeros()as usize;}if count-rank<=4{for _ in 0..count-rank-1{word&=!(1<<(u64::BITS as usize-1-word.leading_zeros()as usize));}return u64::BITS as usize-1-word.leading_zeros()as usize;}let mut offset=0;let mut width=u64::BITS as usize/2;while width!=0{let mask=u64::MAX>>(u64::BITS as usize-width);let count=(word&mask).count_ones()as usize;if rank<count{word&=mask;}else{word>>=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<<rank,word).trailing_zeros()as usize}}#[derive(Debug,Clone)]#[repr(C)]pub struct BitVectorBlock{pub bits:u64,pub rank:usize}#[derive(Debug,Clone)]pub struct BitVector{blocks:Vec<BitVectorBlock>,len:usize,sum:usize,select_samples:[Vec<usize>;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<<count)-1)};blocks.push(BitVectorBlock{bits,rank:sum});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)}fn from_blocks(blocks:Vec<BitVectorBlock>,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<<offset)!=0,block.rank+(block.bits&!(u64::MAX<<offset)).count_ones()as usize)}#[inline]fn rank1(&self,k:usize)->usize{self.access_rank1(k).1}fn select1(&self,k:usize)->Option<usize>{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<usize>{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 FromIterator<bool>for BitVector{fn from_iter<I:IntoIterator<Item=bool>>(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: 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<WaveletMatrixQuadBlock>,starts:[usize;4],select_samples:[Vec<usize>;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<usize>;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<u64>{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)<<i;i+=8;}for(j,&value)in chunk[i..].iter().enumerate(){word|=(((value>>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<const SCALE:i32>(base:*const i64,index:__m512i)->__m512i{_mm512_i64gather_epi64::<SCALE>(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::<i64>();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::<i64>();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<T>{len:usize,bit_length:usize,zeros:Vec<usize>,bit_vectors:Vec<BitVector>,quad_vectors:Vec<WaveletMatrixQuadVector>,compress:VecCompress<T>,#[cfg(target_arch="x86_64")]backend:super::SimdBackend}impl<T>WaveletMatrix<T>where T:Ord+Clone{pub fn new(v:Vec<T>)->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<I:Copy>(indices:&[I],bit:impl Fn(I)->bool)->Vec<u64>{indices.chunks(64).map(|chunk|{chunk.iter().enumerate().fold(0,|word,(i,&index)|word|((bit(index)as u64)<<i))}).collect()}fn partition<I:Copy>(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<I:Copy>(v:Vec<T>,code:impl Fn(usize)->I,index:impl Fn(I)->usize,pack:impl Fn(&[I],usize)->Vec<u64>,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<quad_bits&&(d%2==1||d+1==quad_bits){if d%2==1{let low=pack(&indices,d-1);quad_vectors.push(WaveletMatrixQuadVector::from_words(&low,Some(&words),len));}else{quad_vectors.push(WaveletMatrixQuadVector::from_words(&words,None,len));}}let bits=BitVector::from_words(&words,len);let zero_count=bits.rank0(len);if d==0{zeros.push(zero_count);bit_vectors.push(bits);break;}partition(&indices,&words,zero_count,&mut next);zeros.push(zero_count);bit_vectors.push(bits);mem::swap(&mut indices,&mut next);}Self{len,bit_length,zeros,bit_vectors,quad_vectors,compress,#[cfg(target_arch="x86_64")]backend:match super::simd_backend(){super::SimdBackend::Avx512 if!is_x86_feature_detected!("avx512vpopcntdq")=>{super::SimdBackend::Avx2}backend=>backend,}}}pub fn new_with_init<F>(v:Vec<T>,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+1<bits{let block=&vector.blocks[k/64];let high=block.rank[2]as usize+block.rank[3]as usize+(block.hi&!(u64::MAX<<(k%64))).count_ones()as usize;let middle=if digit&2==0{k-high}else{this.zeros[this.level(d+1)]+high};f(d+1,middle,value.clone());}k=vector.starts[digit]+rank;f(d,k,value.clone());}}return this;}for(mut k,value)in v.into_iter().enumerate(){for d in(0..this.bit_length).rev(){let level=this.level(d);let(bit,rank1)=this.bit_vectors[level].access_rank1(k);k=if bit{this.zeros[level]+rank1}else{k-rank1};f(d,k,value.clone());}}this}fn level(&self,d:usize)->usize{self.bit_length-1-d}fn rank1(&self,level:usize,k:usize)->usize{self.bit_vectors[level].rank1(k)}fn reorder<U>(&self,level:usize,current:Vec<U>)->Vec<U>{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<usize>)->Range<usize>{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<const OP:u8>(&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::<OP>(states,count);}return;}self.batch_inner::<OP>(states,count);}#[cfg(target_arch="x86_64")]#[target_feature(enable="popcnt")]unsafe fn batch_popcnt<const OP:u8>(&self,states:&mut[[usize;4];16],count:usize){self.batch_inner::<OP>(states,count);}#[inline(always)]fn batch_inner<const OP:u8>(&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)<<d;}if OP==QUANTILE{state[2]-=if bit{count0}else{0};}if OP==RANK_LESSTHAN{state[3]+=if bit{count0}else{0};}}}}#[doc=" get k-th value"]pub fn access(&self,mut k: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)<<d;k=if bit{self.zeros[level]+rank1}else{k-rank1};}self.compress.values()[idx].clone()}#[doc=" Returns the values at `indices` in input order."]pub fn access_batch(&self,indices:impl IntoIterator<Item=usize>)->Vec<T>{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<self.len);state[0]=index;}self.batch::<ACCESS>(&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>)->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<Item=(T,Range<usize>)>)->Vec<usize>{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::<RANK>(&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<usize>{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<usize>,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)<<d;range.start=if bit{self.zeros[level]+start1}else{start0};range.end=if bit{self.zeros[level]+end1}else{end0};}self.compress.values()[idx].clone()}#[inline(always)]fn quad_quantile(&self,mut range:Range<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<Item=(Range<usize>,usize)>)->Vec<T>{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<range.len());*state=[range.start,range.end,*k,0];}if let[(range,k)]=queries{result.push(self.quantile(range.clone(),*k));continue;}self.batch::<QUANTILE>(&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<usize>,ranks:&[usize])->Vec<T>{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<range.len()));if ranks.is_empty(){return Vec::new();}if ranks[0]==ranks[ranks.len()-1]{return vec![self.quantile(range,ranks[0]);ranks.len()];}let use_batch=ranks.len()<1024;#[cfg(target_arch="x86_64")]let use_batch=use_batch||self.backend==super::SimdBackend::Avx512;if use_batch&&ranks.len()<self.compress.size().div_ceil(16){return self.quantile_batch(ranks.iter().map(|&k|(range.clone(),k)));}let mut result=Vec::with_capacity(ranks.len());if!self.quad_vectors.is_empty(){let mut stack=vec![(range,0,0,0,ranks)];while let Some((range,base,index,level,ranks))=stack.pop(){if level==self.quad_vectors.len(){result.resize(result.len()+ranks.len(),self.compress.values()[index].clone());continue;}if ranks.len()==1{result.push(self.quad_quantile(range,ranks[0]-base,level,index));continue;}let vector=&self.quad_vectors[level];let start=vector.ranks(range.start);let end=vector.ranks(range.end);let mut split=base+range.len();let mut requested=ranks.len();for digit in(0..4).rev(){split-=end[digit]-start[digit];let boundary=ranks[..requested].partition_point(|&k|k<split);if boundary<requested{stack.push((vector.starts[digit]+start[digit]..vector.starts[digit]+end[digit],split,index*4+digit,level+1,&ranks[boundary..requested]));}requested=boundary;}}return result;}let mut stack=vec![(range,0,0,self.bit_length,ranks)];while let Some((range,base,index,bits,ranks))=stack.pop(){if bits==0{result.resize(result.len()+ranks.len(),self.compress.values()[index].clone());continue;}let d=bits-1;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 split=base+end0-start0;let boundary=ranks.partition_point(|&k|k<split);if boundary<ranks.len(){stack.push((self.zeros[level]+start1..self.zeros[level]+end1,split,index|(1<<d),d,&ranks[boundary..]));}if boundary!=0{stack.push((start0..end0,base,index,d,&ranks[..boundary]));}}result}#[doc=" get k-th smallest value out of range"]pub fn quantile_outer(&self,mut range:Range<usize>,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<count{break;}k-=count;digit+=1;}index=index*4+digit;range=vector.starts[digit]+start[digit]..vector.starts[digit]+end[digit];outer=vector.starts[digit]+outer_start[digit]..vector.starts[digit]+outer_end[digit];}return self.compress.values()[index].clone();}let mut idx=0;let mut orange=0..self.len;for d in(0..self.bit_length-self.compress.size().is_power_of_two()as usize).rev(){let level=self.level(d);let range_start1=self.rank1(level,range.start);let range_end1=self.rank1(level,range.end);let outer_start1=self.rank1(level,orange.start);let outer_end1=self.rank1(level,orange.end);let range_start0=range.start-range_start1;let range_end0=range.end-range_end1;let outer_start0=orange.start-outer_start1;let outer_end0=orange.end-outer_end1;let z=(outer_end0-outer_start0)-(range_end0-range_start0);if z<=k{k-=z;idx|=1<<d;range.start=self.zeros[level]+range_start1;range.end=self.zeros[level]+range_end1;orange.start=self.zeros[level]+outer_start1;orange.end=self.zeros[level]+outer_end1;}else{range.start=range_start0;range.end=range_end0;orange.start=outer_start0;orange.end=outer_end0;}}self.compress.values()[idx].clone()}#[doc=" the number of value less than val in range"]pub fn rank_lessthan(&self,val:T,range:Range<usize>)->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<Item=(T,Range<usize>)>)->Vec<usize>{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::<RANK_LESSTHAN>(&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<usize>,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<usize>,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<T>,mut range:Range<usize>)->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<Item=(Range<T>,Range<usize>)>)->Vec<usize>{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<upper{states[i*2]=[range.start,range.end,lower,0];states[i*2+1]=[range.start,range.end,upper,0];}}self.batch::<RANK_LESSTHAN>(&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<F>(&self,val:T,mut range:Range<usize>,mut f:F)where F:FnMut(usize,Range<usize>){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<M>(&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<M::T>=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&current[..self.zeros[level]]{acc=M::operate(&acc,w);prefix.push(acc.clone());}}WaveletMatrixFold{wavelet_matrix:self,prefix,offsets}}pub fn build_point_add<M>(&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(&current[..self.zeros[level]]));}WaveletMatrixPointAdd{wavelet_matrix:self,bits}}}pub struct WaveletMatrixPointAdd<'a,T,M>where T:Ord+Clone,M:AbelianGroup{wavelet_matrix:&'a WaveletMatrix<T>,bits:Vec<BinaryIndexedTree<M>>}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<self.wavelet_matrix.len);for d in(0..self.wavelet_matrix.bit_length).rev(){let level=self.wavelet_matrix.level(d);let(bit,rank1)=self.wavelet_matrix.bit_vectors[level].access_rank1(index);if bit{index=self.wavelet_matrix.zeros[level]+rank1;}else{index-=rank1;self.bits[level].update(index,value.clone());}}}pub fn fold_lessthan(&self,value:T,range:Range<usize>)->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<T>,range:Range<usize>)->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<usize>,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<T>,prefix:Vec<M::T>,offsets:Vec<usize>}impl<'a,T,M>WaveletMatrixFold<'a,T,M>where T:Ord+Clone,M:AbelianGroup{pub fn fold_lessthan(&self,val:T,range:Range<usize>)->M::T{self.fold_lessthan_with_count(val,range).1}pub fn fold_lessthan_with_count(&self,val:T,range:Range<usize>)->(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<T>,range:Range<usize>)->M::T{self.fold_range_with_count(valrange,range).1}pub fn fold_range_with_count(&self,valrange:Range<T>,mut range:Range<usize>)->(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<usize>)->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<const N:usize>(&self,indices:[usize;N],mut ranges:[Range<usize>;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<Item=(T,Range<usize>)>)->Vec<M::T>{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{}impl<S>SemiGroup 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<Item=bool>;fn bits(self)->Self::Iter;}pub trait SignedExpBits{type T:ExpBits;fn neg_and_bits(self)->(bool,Self::T);}pub struct Bits<T>{n:T}macro_rules!impl_exp_bits_for_uint{($($t:ty)*)=>{$(impl Iterator for Bits<$t>{type Item=bool;fn next(&mut self)->Option<bool>{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<E>(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<I>(iter:I)->Self::T where I:IntoIterator<Item=Self::T>{let mut iter=iter.into_iter();if let Some(item)=iter.next(){iter.fold(item,|acc,x|Self::operate(&acc,&x))}else{Self::unit()}}}impl<M>Monoid 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<E>(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}}}impl<G>Group 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{}impl<M>AbelianMonoid for M where M:Monoid+Commutative{}#[doc=" commutative group"]pub trait AbelianGroup:Group+Commutative{}impl<G>AbelianGroup 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{}impl<M>IdempotentMonoid 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: 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 Compressor<T>where Self:FromIterator<T>,T:Ord{fn index_exact(&self,index:&T)->Option<usize>;fn size(&self)->usize;}pub trait OrderedCompressor<T>:Compressor<T>where T:Ord{fn index_lower_bound(&self,index:&T)->usize;}#[derive(Debug,Clone)]pub struct VecCompress<T>{data:Vec<T>}impl<T>VecCompress<T>{pub fn from_sorted_unique(data:Vec<T>)->Self{Self{data}}pub fn values(&self)->&[T]{&self.data}}impl<T>FromIterator<T>for VecCompress<T>where T:Ord{fn from_iter<I>(iter:I)->Self where I:IntoIterator<Item=T>{let mut data:Vec<_>=iter.into_iter().collect();data.sort_unstable();data.dedup();Self{data}}}impl<T>Compressor<T>for VecCompress<T>where T:Ord{fn index_exact(&self,index:&T)->Option<usize>{self.data.binary_search(index).ok()}fn size(&self)->usize{self.data.len()}}impl<T>OrderedCompressor<T>for VecCompress<T>where T:Ord{fn index_lower_bound(&self,index:&T)->usize{self.data.partition_point(|x|x<index)}}#[derive(Clone)]pub struct HashCompress<T>{data:HashMap<T,usize>}impl<T>Debug for HashCompress<T>where T:Debug+Eq+Hash{fn fmt(&self,f:&mut fmt::Formatter<'_>)->fmt::Result{f.debug_struct("HashCompress").field("data",&self.data).finish()}}impl<T>FromIterator<T>for HashCompress<T>where T:Ord+Hash{fn from_iter<I>(iter:I)->Self where I:IntoIterator<Item=T>{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}}}impl<T>Compressor<T>for HashCompress<T>where T:Ord+Hash{fn index_exact(&self,index:&T)->Option<usize>{self.data.get(index).copied()}fn size(&self)->usize{self.data.len()}}}
// 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.}{<T:Zero>Zero zero std::num::Wrapping<T>,Self(T::zero())}{<T:One>One one std::num::Wrapping<T>,Self(T::one())});}