// codesnip-guard: discrete_logarithm
#[cfg_attr(any(),rustfmt::skip)]pub use self::discrete_logarithm::{discrete_logarithm,discrete_logarithm_prime_mod};#[cfg_attr(any(),rustfmt::skip)]mod discrete_logarithm{use super::{BarrettReduction,PrimeList,Xorshift,check_primitive_root,gcd,lcm,modinv,prime_factors,primitive_root};use std::{cell::UnsafeCell,collections::HashMap};fn pow(x:u64,mut y:u64,br:&BarrettReduction<u128>)->u64{let mut x=x as u128;let mut z:u128=1;while y>0{if y&1==1{z=br.rem(z*x);}x=br.rem(x*x);y>>=1;}z as u64}fn solve_linear_congruence(a:u64,b:u64,m:u64)->Option<(u64,u64)>{let g=gcd(a,m);if!b.is_multiple_of(g){return None;}let(a,b,m)=(a/g,b/g,m/g);Some(((b as u128*modinv(a,m)as u128%m as u128)as _,m))}fn solve_linear_congruences<I>(abm:I)->Option<(u64,u64)>where I:IntoIterator<Item=(u64,u64,u64)>{let mut x=0u64;let mut m0=1u64;for(a,b,m)in abm{let mut b=b+m-a*x%m;if b>=m{b-=m;}let a=a*m0;let g=gcd(a,m);if!b.is_multiple_of(g){return None;}let(a,b,m)=(a/g,b/g,m/g);x+=(b as u128*modinv(a,m)as u128%m as u128*m0 as u128)as u64;m0*=m;}Some((x,m0))}#[derive(Debug)]struct IndexCalculus{primes:PrimeList,br_primes:Vec<BarrettReduction<u64>>,ic:HashMap<u64,IndexCalculusWithPrimitiveRoot>}impl IndexCalculus{fn new()->Self{Self{primes:PrimeList::new(2),br_primes:Default::default(),ic:Default::default()}}fn discrete_logarithm(&mut self,a:u64,b:u64,p:u64)->Option<(u64,u64)>{let lim=((((p as f64).log2()*(p as f64).log2().log2()).sqrt()/2.0+1.).exp2()*0.9)as u32;self.primes.reserve(lim);let prime_count=self.primes.primes_lte(lim).count();self.br_primes.extend(self.primes.primes_lte(lim).skip(self.br_primes.len()).map(|p|BarrettReduction::<u64>::new(p.into())));let br_primes=&self.br_primes[..prime_count];self.ic.entry(p).or_insert_with(| |IndexCalculusWithPrimitiveRoot::new(p,br_primes)).discrete_logarithm(a,b,br_primes)}}const A:[u32;150]=[62,61,60,60,59,58,58,58,57,56,56,56,56,55,55,55,54,54,54,53,53,53,53,52,52,52,52,52,52,51,50,50,50,50,49,49,49,48,48,48,48,48,47,47,47,47,47,47,47,47,47,47,47,47,47,47,45,42,42,41,41,41,41,41,41,41,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,38,38,38,38,38,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,22,22,22,22,22,22,22,22,22,22];fn factorize_smooth(mut x:u64,row:&mut[u64],br_primes:&[BarrettReduction<u64>])->bool{for(j,(&br,r))in br_primes.iter().zip(row).enumerate(){*r=0;loop{let(div,rem)=br.div_rem(x);if rem!=0{break;}*r+=1;x=div;}if j<150&&x>=(1u64<<A[j]){break;}}x==1}#[derive(Debug)]struct QdrtPowPrec{br_qdrt:BarrettReduction<u64>,p0:Vec<u64>,p1:Vec<u64>,p2:Vec<u64>,p3:Vec<u64>}impl QdrtPowPrec{fn new(a:u64,ord:u64,br:&BarrettReduction<u128>)->Self{let qdrt=(ord as f64).powf(0.25).ceil()as u64;let br_qdrt=BarrettReduction::<u64>::new(qdrt);let mut p0=Vec::with_capacity(qdrt as usize);let mut p1=Vec::with_capacity(qdrt as usize);let mut p2=Vec::with_capacity(qdrt as usize);let mut p3=Vec::with_capacity(qdrt as usize);let mut acc=1u64;for _ in 0..qdrt{p0.push(acc);acc=br.rem(acc as u128*a as u128)as u64;}let a=acc;acc=1;for _ in 0..qdrt{p1.push(acc);acc=br.rem(acc as u128*a as u128)as u64;}let a=acc;acc=1;for _ in 0..qdrt{p2.push(acc);acc=br.rem(acc as u128*a as u128)as u64;}let a=acc;acc=1;for _ in 0..qdrt{p3.push(acc);acc=br.rem(acc as u128*a as u128)as u64;}Self{br_qdrt,p0,p1,p2,p3}}fn pow(&self,mut k:u64,br:&BarrettReduction<u128>)->u64{let(a,b)=self.br_qdrt.div_rem(k);let mut x=self.p0[b as usize];k=a;if k>0{let(a,b)=self.br_qdrt.div_rem(k);x=br.rem(x as u128*self.p1[b as usize]as u128)as u64;k=a;}if k>0{let(a,b)=self.br_qdrt.div_rem(k);x=br.rem(x as u128*self.p2[b as usize]as u128)as u64;k=a;}if k>0{let(_,b)=self.br_qdrt.div_rem(k);x=br.rem(x as u128*self.p3[b as usize]as u128)as u64;}x}}fn index_calculus_for_primitive_root(p:u64,ord:u64,br_primes:&[BarrettReduction<u64>],prec:&QdrtPowPrec)->Vec<u64>{let br_ord=BarrettReduction::<u128>::new(ord as u128);let mul=|x:u64,y:u64|br_ord.rem(x as u128*y as u128)as u64;let sub=|x:u64,y:u64|if x<y{x+ord-y}else{x-y};let pc=br_primes.len();let mut mat:Vec<Vec<u64>>=vec![];let mut rows:Vec<Vec<u64>>=vec![];let mut rng=Xorshift::default();let br=BarrettReduction::<u128>::new(p as u128);for i in 0..pc{for ri in 0usize..{let mut row=vec![0u64;pc+1];let mut kk=rng.rand(ord-1)+1;let mut gkk=prec.pow(kk,&br);let mut k=kk;let mut gk=gkk;while ri>=rows.len(){row[pc]=k;if factorize_smooth(gk,&mut row,br_primes){rows.push(row);break;}if k+kk<ord{k+=kk;gk=br.rem(gk as u128*gkk as u128)as u64;}else{kk=rng.rand(ord-1)+1;gkk=prec.pow(kk,&br);k=kk;gk=gkk;}}let row=&mut rows[ri];for j in 0..i{if row[j]!=0{let b=mul(modinv(mat[j][j],ord),row[j]);for(r,a)in row[j..].iter_mut().zip(&mat[j][j..]){*r=sub(*r,mul(*a,b));}}assert_eq!(row[j],0);}if gcd(row[i],ord)==1{let last=rows.len()-1;rows.swap(ri,last);mat.push(rows.pop().unwrap());break;}}}for i in(0..pc).rev(){for j in i+1..pc{mat[i][pc]=sub(mat[i][pc],mul(mat[i][j],mat[j][pc]));}mat[i][pc]=mul(mat[i][pc],modinv(mat[i][i],ord));}(0..pc).map(|i|mat[i][pc]).collect()}#[derive(Debug)]struct IndexCalculusWithPrimitiveRoot{p:u64,ord:u64,prec:QdrtPowPrec,coeff:Vec<u64>}impl IndexCalculusWithPrimitiveRoot{fn new(p:u64,br_primes:&[BarrettReduction<u64>])->Self{let ord=p-1;let g=primitive_root(p);let br=BarrettReduction::<u128>::new(p as u128);let prec=QdrtPowPrec::new(g,ord,&br);let coeff=index_calculus_for_primitive_root(p,ord,br_primes,&prec);Self{p,ord,prec,coeff}}fn index_calculus(&self,a:u64,br_primes:&[BarrettReduction<u64>])->Option<u64>{let p=self.p;let ord=self.ord;let br=BarrettReduction::<u128>::new(p as u128);let a=br.rem(a as _)as u64;if a==1{return Some(0);}if p==2{return None;}let mut rng=Xorshift::new();let mut row=vec![0u64;br_primes.len()];let mut kk=rng.rand(ord-1)+1;let mut gkk=self.prec.pow(kk,&br);let mut k=kk;let mut gk=br.rem(gkk as u128*a as u128)as u64;loop{if factorize_smooth(gk,&mut row,br_primes){let mut res=ord-k;for(&c,&r)in self.coeff.iter().zip(&row){for _ in 0..r{res+=c;if res>=ord{res-=ord;}}}return Some(res);}if k+kk<ord{k+=kk;gk=br.rem(gk as u128*gkk as u128)as u64;}else{kk=rng.rand(ord-1)+1;gkk=self.prec.pow(kk,&br);k=kk;gk=br.rem(gkk as u128*a as u128)as u64;}}}fn discrete_logarithm(&self,a:u64,b:u64,br_primes:&[BarrettReduction<u64>])->Option<(u64,u64)>{let p=self.p;let ord=self.ord;let br=BarrettReduction::<u128>::new(p as u128);let a=br.rem(a as _)as u64;let b=br.rem(b as _)as u64;if a==0{return if b==0{Some((1,1))}else{None};}if b==0{return None;}let x=self.index_calculus(a,br_primes)?;let y=self.index_calculus(b,br_primes)?;solve_linear_congruence(x,y,ord)}}thread_local!(static IC:UnsafeCell<IndexCalculus> =UnsafeCell::new(IndexCalculus::new()););pub fn discrete_logarithm_prime_mod(a:u64,b:u64,p:u64)->Option<u64>{IC.with(|ic|unsafe{&mut*ic.get()}.discrete_logarithm(a,b,p)).map(|t|t.0)}#[doc=" a^x ≡ b (mod n), a has order p^e"]fn pohlig_hellman_prime_power_order(a:u64,b:u64,n:u64,p:u64,e:u32)->Option<u64>{let br=BarrettReduction::<u128>::new(n as u128);let mul=|x:u64,y:u64|br.rem(x as u128*y as u128)as u64;let block_size=(p as f64).sqrt().ceil()as u64;let mut baby=HashMap::<u64,u64>::new();let g=pow(a,p.pow(e-1),&br);let mut xj=1;for j in 0..block_size{baby.entry(xj).or_insert(j);xj=mul(xj,g);}let xi=modinv(xj,n);let mut t=0u64;for k in 0..e{let mut h=pow(mul(modinv(pow(a,t,&br),n),b),p.pow(e-1-k),&br);let mut ok=false;for i in(0..block_size*block_size).step_by(block_size as usize){if let Some(j)=baby.get(&h){t+=(i+j)*p.pow(k);ok=true;break;}h=mul(h,xi);}if!ok{return None;}}Some(t)}#[doc=" a^x ≡ b (mod p^e)"]fn discrete_logarithm_prime_power(a:u64,b:u64,p:u64,e:u32)->Option<(u64,u64)>{assert_ne!(p,0);assert_ne!(e,0);let n=p.pow(e);assert!(a<n);assert!(b<n);assert_eq!(gcd(a,p),1);if p==1{return Some((0,1));}if a==0{return if b==0{Some((1,1))}else{None};}if b==0{return None;}if e==1{return IC.with(|ic|unsafe{&mut*ic.get()}.discrete_logarithm(a,b,p));}let br=BarrettReduction::<u128>::new(n as _);if p==2{if e>=3{if a%4==1&&b%4!=1{return None;}let aa=if a%4==1{a}else{n-a};let bb=if b%4==1{b}else{n-b};let g=5;let ord=n/4;let x=pohlig_hellman_prime_power_order(g,aa,n,p,e-2)?;let y=pohlig_hellman_prime_power_order(g,bb,n,p,e-2)?;let t=solve_linear_congruence(x,y,ord)?;match(a%4==1,b%4==1){(true,true)=>Some(t),(false,true)if t.0%2==0=>Some((t.0,lcm(t.1,2))),(false,false)if t.0%2==1=>Some((t.0,lcm(t.1,2))),(false,false)if a==b=>Some((1,lcm(t.1,2))),_=>None,}}else if a==1{if b==1{Some((0,1))}else{None}}else{assert_eq!(a,3);if b==1{Some((0,2))}else if b==3{Some((1,2))}else{None}}}else{let ord=n-n/p;let pf_ord=prime_factors(ord);let g=(2..).find(|&g|check_primitive_root(g,ord,&br,&pf_ord)).unwrap();let mut pf_p=prime_factors(p-1);pf_p.push((p,e-1));let mut abm=vec![];for(q,c)in pf_p{let m=q.pow(c);let d=ord/m;let gg=pow(g,d,&br);let aa=pow(a,d,&br);let bb=pow(b,d,&br);let x=pohlig_hellman_prime_power_order(gg,aa,n,q,c)?;let y=pohlig_hellman_prime_power_order(gg,bb,n,q,c)?;abm.push((x,y,m));}solve_linear_congruences(abm)}}#[doc=" a^x ≡ b (mod n)"]pub fn discrete_logarithm(a:u64,b:u64,n:u64)->Option<u64>{let a=a%n;let b=b%n;let d=2.max(64-n.leading_zeros()as u64);let mut pw=1%n;for i in 0..d{if pw==b{return Some(i);}pw=(pw as u128*a as u128%n as u128)as u64;}let g=gcd(pw,n);if!b.is_multiple_of(g){return None;}let n=n/g;let b=(b as u128*modinv(pw,n)as u128%n as u128)as u64;let pf=prime_factors(n);let mut abm=vec![];for(p,e)in pf{let q=p.pow(e);let x=discrete_logarithm_prime_power(a%q,b%q,p,e)?;abm.push((1,x.0,x.1));}solve_linear_congruences(abm).map(|x|x.0+d)}}
// codesnip-guard: BarrettReduction
#[cfg_attr(any(),rustfmt::skip)]pub use self::barrett_reduction::{BarrettReduction,Barrettable};#[cfg_attr(any(),rustfmt::skip)]mod barrett_reduction{use super::{One,Zero};use std::ops::{Add,Mul,Sub};#[derive(Debug,Clone,Copy)]pub struct BarrettReduction<T>{m:T,im:T}impl<T>BarrettReduction<T>where T:Barrettable{pub fn new(m:T)->Self{Self{m,im:T::inv_mod_approx(m)}}pub const fn new_with_im(m:T,im:T)->Self{Self{m,im}}pub const fn get_mod(&self)->T{self.m}pub fn div_rem(&self,a:T)->(T,T){T::barrett_reduce(a,self.m,self.im)}pub fn div(&self,a:T)->T{self.div_rem(a).0}pub fn rem(&self,a:T)->T{self.div_rem(a).1}}pub trait Barrettable:Sized+Copy+PartialOrd+Zero+One+Add<Output=Self>+Sub<Output=Self>+Mul<Output=Self>{fn inv_mod_approx(m:Self)->Self;fn div_approx(self,im:Self)->Self;fn barrett_reduce(self,m:Self,im:Self)->(Self,Self){if m==Self::one(){return(self,Self::zero());}let q=self.div_approx(im);let r=self-q*m;if m<=r{(q+Self::one(),r-m)}else{(q,r)}}}impl Barrettable for u32{fn inv_mod_approx(m:Self)->Self{!0/m}fn div_approx(self,im:Self)->Self{((self as u64*im as u64)>>32)as u32}}impl Barrettable for u64{fn inv_mod_approx(m:Self)->Self{!0/m}fn div_approx(self,im:Self)->Self{((self as u128*im as u128)>>64)as u64}}impl Barrettable for u128{fn inv_mod_approx(m:Self)->Self{!0/m}fn div_approx(self,im:Self)->Self{const MASK64:u128=0xffff_ffff_ffff_ffff;let au=self>>64;let ad=self&MASK64;let imu=im>>64;let imd=im&MASK64;let res=au*imu;let x=((ad*imd)>>64)+au*imd;let(x,c)=x.overflowing_add(ad*imu);res+((c as u128)<<64)+(x>>64)}}}
// codesnip-guard: PrimeList
#[cfg_attr(any(),rustfmt::skip)]pub use self::prime_list::{PrimeList,with_prime_list};#[cfg_attr(any(),rustfmt::skip)]mod prime_list{use std::{cell::UnsafeCell,mem::replace,ops::Range};const WHEEL_PRIMES:[u32;4]=[2,3,5,7];const PERIOD:u32=210;const COPRIME:usize=48;const SQRT_THRESHOLD:u32=1<<16;const fn coprime_to_wheel(x:u32)->bool{!x.is_multiple_of(2)&&!x.is_multiple_of(3)&&!x.is_multiple_of(5)&&!x.is_multiple_of(7)}const fn residues()->[u8;COPRIME]{let mut result=[0;COPRIME];let mut i=1;let mut j=0;while i<PERIOD{if coprime_to_wheel(i){result[j]=i as u8;j+=1;}i+=2;}result}const RESIDUES:[u8;COPRIME]=residues();const fn states()->[u8;PERIOD as usize]{let mut result=[0;PERIOD as usize];let mut i=0;let mut j=0;while i<PERIOD{result[i as usize]=j;if coprime_to_wheel(i){j+=1;}i+=1;}result}const STATES:[u8;PERIOD as usize]=states();const fn additions()->[u8;PERIOD as usize]{let mut result=[0;PERIOD as usize];let mut i=0;while i<PERIOD{let mut add=1;while!coprime_to_wheel(i+add){add+=1;}result[i as usize]=add as u8;i+=1;}result}const ADDITIONS:[u8;PERIOD as usize]=additions();const fn gaps()->[u8;COPRIME]{let mut result=[0;COPRIME];let mut i=0;while i<COPRIME{result[i]=ADDITIONS[RESIDUES[i]as usize];i+=1;}result}const GAPS:[u8;COPRIME]=gaps();const fn to_ordinal(x:u32)->u32{x/PERIOD*COPRIME as u32+STATES[(x%PERIOD)as usize]as u32}const fn to_value(x:u32)->u32{x/COPRIME as u32*PERIOD+RESIDUES[x as usize%COPRIME]as u32}const fn ordinal_to_value()->[u16;256]{let mut result=[0;256];let mut i=0;while i<result.len(){result[i]=to_value(i as u32)as u16;i+=1;}result}const ORDINAL_TO_VALUE:[u16;256]=ordinal_to_value();const fn sqrt_bits()->[u64;SQRT_THRESHOLD as usize/128]{let mut result=[!0;SQRT_THRESHOLD as usize/128];let ordinal=to_ordinal(1)as usize;result[ordinal/64]&=!(1<<(ordinal%64));let mut i=RESIDUES[1]as u32;while to_ordinal(i*i)<SQRT_THRESHOLD/2{let ordinal=to_ordinal(i)as usize;if result[ordinal/64]>>(ordinal%64)&1!=0{let mut k=i;while to_ordinal(i*k)<SQRT_THRESHOLD/2{let ordinal=to_ordinal(i*k)as usize;result[ordinal/64]&=!(1<<(ordinal%64));k+=ADDITIONS[(k%PERIOD)as usize]as u32;}}i+=ADDITIONS[(i%PERIOD)as usize]as u32;}result}const SQRT_BITS:[u64;SQRT_THRESHOLD as usize/128]=sqrt_bits();const fn count_sqrt_primes()->usize{let mut result=0;let mut i=RESIDUES[1]as u32;while i<SQRT_THRESHOLD{let ordinal=to_ordinal(i)as usize;result+=(SQRT_BITS[ordinal/64]>>(ordinal%64)&1)as usize;i+=ADDITIONS[(i%PERIOD)as usize]as u32;}result}const SQRT_PRIME_COUNT:usize=count_sqrt_primes();const fn sqrt_primes()->[u32;SQRT_PRIME_COUNT]{let mut result=[0;SQRT_PRIME_COUNT];let mut i=RESIDUES[1]as u32;let mut j=0;while i<SQRT_THRESHOLD{let ordinal=to_ordinal(i)as usize;if SQRT_BITS[ordinal/64]>>(ordinal%64)&1!=0{result[j]=i;j+=1;}i+=ADDITIONS[(i%PERIOD)as usize]as u32;}result}static SQRT_PRIMES:[u32;SQRT_PRIME_COUNT]=sqrt_primes();struct Wheel{mask:Vec<u64>,product:u32}impl Wheel{fn new(primes:&[u32],product:u32)->Self{let mut mask=vec![!0;to_ordinal(product)as usize/64];for&p in primes{let mut k=1;while p*k<product{let ordinal=to_ordinal(p*k)as usize;mask[ordinal/64]&=!(1<<(ordinal%64));k+=ADDITIONS[(k%PERIOD)as usize]as u32;}}Self{mask,product}}}fn make_wheels()->(Vec<Wheel>,usize){const MAX_WHEEL_SIZE:u32=1<<20;const BASE:u32=(PERIOD*64)>>(WHEEL_PRIMES.len()-2);let mut product=BASE;let mut current=vec![];let mut wheels=vec![];for(i,&p)in SQRT_PRIMES.iter().enumerate(){if product*p>MAX_WHEEL_SIZE{wheels.push(Wheel::new(&current,product));current.clear();current.push(p);product=BASE*p;if product>MAX_WHEEL_SIZE{return(wheels,i);}}else{current.push(p);product*=p;}}unreachable!()}fn sieve_dense(bits:&mut[u64],l:u32,r:u32,wheel:&Wheel){let mut left=l as usize/64;let right=(r as usize).div_ceil(64);while left+wheel.mask.len()<=right{for(value,&mask)in bits[left..left+wheel.mask.len()].iter_mut().zip(&wheel.mask){*value&=mask;}left+=wheel.mask.len();}for(value,&mask)in bits[left..right].iter_mut().zip(&wheel.mask){*value&=mask;}}fn ordinal_steps()->Vec<[u32;COPRIME*2]>{SQRT_PRIMES.iter().map(|&p|{let mut result=[0;COPRIME*2];let mut last=to_ordinal(p);for i in 0..COPRIME{let next=to_ordinal(p*(RESIDUES[i]as u32+GAPS[i]as u32));result[i]=next-last;result[i+COPRIME]=next-last;last=next;}result}).collect()}fn sieve_sparse(bits:&mut[u64],mut left:u32,right:u32,prime_index:usize,mut state:u8,steps:&[[u32;COPRIME*2]])->(u32,u8){let p=SQRT_PRIMES[prime_index];while left+p*COPRIME as u32<=right{for _ in 0..COPRIME{let ordinal=left as usize;bits[ordinal/64]&=!(1<<(ordinal%64));left+=steps[prime_index][state as usize];state+=1;}state-=COPRIME as u8;}while left<right{let ordinal=left as usize;bits[ordinal/64]&=!(1<<(ordinal%64));left+=steps[prime_index][state as usize];state+=1;}if state>=COPRIME as u8{state-=COPRIME as u8;}(left,state)}#[derive(Debug,Clone)]pub struct PrimeList{bits:Vec<u64>,bit_len:usize,max_n:u32,prime_count:usize}impl Default for PrimeList{fn default()->Self{Self{bits:vec![],bit_len:0,max_n:1,prime_count:0}}}impl PrimeList{pub fn new(max_n:u32)->Self{let mut self_:Self=Default::default();self_.reserve(max_n);self_}pub fn primes(&self)->PrimeListIter<'_>{self.primes_lte(self.max_n)}pub fn len(&self)->usize{self.prime_count}pub fn is_empty(&self)->bool{self.prime_count==0}pub fn primes_lte(&self,n:u32)->PrimeListIter<'_>{assert!(n<=self.max_n,"expected `n={} <= {}`",n,self.max_n);let bit_len=to_ordinal(n.saturating_add(1))as usize;let words=&self.bits[..self.bits.len().min(bit_len.div_ceil(64))];let last_mask=if bit_len.is_multiple_of(64){!0}else{(1<<(bit_len%64))-1};let(front_word,middle_words,back_word)=match words{[]=>(0,words,0),[word]=>(*word&last_mask,&words[1..],0),[front,middle@..,back]=>(*front,middle,*back&last_mask),};let back_word_ordinal=words.len().saturating_sub(1)as u32*64;PrimeListIter{wheel_indices:0..WHEEL_PRIMES.partition_point(|&p|p<=n)as u8,middle_words,front_word_base:0,front_word_phase:0,front_word,back_word_base:back_word_ordinal/COPRIME as u32*PERIOD,back_word_phase:(back_word_ordinal%COPRIME as u32)as u8,back_word}}pub fn is_prime(&self,n:u32)->bool{assert!(n<=self.max_n,"expected `n={} <= {}`",n,self.max_n);if WHEEL_PRIMES.contains(&n){true}else if!coprime_to_wheel(n){false}else{let ordinal=to_ordinal(n)as usize;ordinal<self.bit_len&&self.bits[ordinal/64]>>(ordinal%64)&1!=0}}pub fn trial_division(&self,n:u64)->PrimeListTrialDivision<'_>{let bound=u64::from(self.max_n).pow(2);assert!(n<=bound,"expected `n={} <= {}`",n,bound);PrimeListTrialDivision{primes:self.primes(),n}}pub fn prime_factors(&self,n:u64)->Vec<(u64,u32)>{self.trial_division(n).collect()}pub fn count_divisors(&self,n:u64)->u64{let mut divisor_cnt=1u64;for(_,cnt)in self.trial_division(n){divisor_cnt*=cnt as u64+1;}divisor_cnt}pub fn divisors(&self,n:u64)->Vec<u64>{let mut d=vec![1u64];for(p,c)in self.trial_division(n){let k=d.len();let mut acc=1;for _ in 0..c{acc*=p;for i in 0..k{d.push(d[i]*acc);}}}d.sort_unstable();d}#[doc=" Extends the prime list up to `max_n`."]pub fn reserve(&mut self,max_n:u32){if max_n<=self.max_n||max_n<2{return;}let limit=max_n.saturating_add(1);self.bit_len=to_ordinal(limit)as usize;if limit<=SQRT_THRESHOLD{self.bits=SQRT_BITS[..self.bit_len.div_ceil(64)].to_vec();}else{self.bits=vec![!0;self.bit_len.div_ceil(64)];let(wheels,medium_primes_begin)=make_wheels();const DENSE_BLOCK:u32=1<<25;for start in(0..limit).step_by(DENSE_BLOCK as usize){let right=start.saturating_add(DENSE_BLOCK).min(limit);for wheel in&wheels{let left=start/wheel.product*wheel.product;sieve_dense(&mut self.bits,to_ordinal(left),to_ordinal(right),wheel);}}let steps=ordinal_steps();let mut positions:Vec<_>=SQRT_PRIMES.iter().map(|&p|to_ordinal(p*p)).collect();let mut states:Vec<_>=SQRT_PRIMES.iter().map(|&p|STATES[(p%PERIOD)as usize]).collect();const SPARSE_BLOCK:u32=1<<22;for start in(0..limit).step_by(SPARSE_BLOCK as usize){let right=to_ordinal(start.saturating_add(SPARSE_BLOCK).min(limit));for i in medium_primes_begin..SQRT_PRIME_COUNT{(positions[i],states[i])=sieve_sparse(&mut self.bits,positions[i],right,i,states[i],&steps);}}for(value,&sqrt_bits)in self.bits.iter_mut().zip(&SQRT_BITS){*value=sqrt_bits;}}self.prime_count=WHEEL_PRIMES.partition_point(|&p|p<=max_n);if let Some((&last,rest))=self.bits.split_last(){self.prime_count+=rest.iter().map(|word|word.count_ones()as usize).sum::<usize>();let last_mask=if self.bit_len.is_multiple_of(64){!0}else{(1<<(self.bit_len%64))-1};self.prime_count+=(last&last_mask).count_ones()as usize;}self.max_n=max_n;}}#[derive(Clone,Debug)]pub struct PrimeListIter<'a>{wheel_indices:Range<u8>,middle_words:&'a[u64],front_word_base:u32,front_word_phase:u8,front_word:u64,back_word_base:u32,back_word_phase:u8,back_word:u64}impl PrimeListIter<'_>{#[inline(always)]fn load_front_word(&mut self)->bool{if let Some((&word,words))=self.middle_words.split_first(){self.middle_words=words;if self.front_word_phase==32{self.front_word_base+=PERIOD*2;self.front_word_phase=0;}else{self.front_word_base+=PERIOD;self.front_word_phase+=16;}self.front_word=word;true}else if self.back_word!=0{self.front_word_base=self.back_word_base;self.front_word_phase=self.back_word_phase;self.front_word=replace(&mut self.back_word,0);true}else{false}}#[inline(always)]fn load_back_word(&mut self)->bool{if let Some((&word,words))=self.middle_words.split_last(){self.middle_words=words;if self.back_word_phase==0{self.back_word_base-=PERIOD*2;self.back_word_phase=32;}else{self.back_word_base-=PERIOD;self.back_word_phase-=16;}self.back_word=word;true}else if self.front_word!=0{self.back_word_base=self.front_word_base;self.back_word_phase=self.front_word_phase;self.back_word=replace(&mut self.front_word,0);true}else{false}}}impl Iterator for PrimeListIter<'_>{type Item=u32;fn next(&mut self)->Option<Self::Item>{if let Some(index)=self.wheel_indices.next(){return Some(WHEEL_PRIMES[index as usize]);}loop{if self.front_word!=0{let bit=self.front_word.trailing_zeros();self.front_word&=self.front_word-1;return Some(self.front_word_base+ORDINAL_TO_VALUE[(self.front_word_phase+bit as u8)as usize]as u32);}if!self.load_front_word(){return None;}}}fn nth(&mut self,mut n:usize)->Option<Self::Item>{if n==0{return self.next();}for index in self.wheel_indices.by_ref(){if n==0{return Some(WHEEL_PRIMES[index as usize]);}n-=1;}loop{let count=self.front_word.count_ones()as usize;if n<count{for _ in 0..n{self.front_word&=self.front_word-1;}return self.next();}n-=count;self.front_word=0;if!self.load_front_word(){return None;}}}}impl DoubleEndedIterator for PrimeListIter<'_>{#[inline]fn next_back(&mut self)->Option<Self::Item>{loop{if self.back_word!=0{let bit=63-self.back_word.leading_zeros();self.back_word-=1<<bit;return Some(self.back_word_base+ORDINAL_TO_VALUE[(self.back_word_phase+bit as u8)as usize]as u32);}if!self.load_back_word(){return self.wheel_indices.next_back().map(|index|WHEEL_PRIMES[index as usize]);}}}fn nth_back(&mut self,mut n:usize)->Option<Self::Item>{if n==0{return self.next_back();}loop{let count=self.back_word.count_ones()as usize;if n<count{for _ in 0..n{let bit=63-self.back_word.leading_zeros();self.back_word-=1<<bit;}return self.next_back();}n-=count;self.back_word=0;if!self.load_back_word(){for index in self.wheel_indices.by_ref().rev(){if n==0{return Some(WHEEL_PRIMES[index as usize]);}n-=1;}return None;}}}}#[derive(Debug,Clone)]pub struct PrimeListTrialDivision<'p>{primes:PrimeListIter<'p>,n:u64}impl Iterator for PrimeListTrialDivision<'_>{type Item=(u64,u32);fn next(&mut self)->Option<Self::Item>{if self.n<=1{return None;}for p in self.primes.by_ref(){let p=u64::from(p);if p*p>self.n{break;}if self.n.is_multiple_of(p){let mut cnt=1u32;self.n/=p;while self.n.is_multiple_of(p){cnt+=1;self.n/=p;}return Some((p,cnt));}}if self.n>1{return Some((replace(&mut self.n,1),1));}None}}pub fn with_prime_list<F>(max_n:u32,f:F)where F:FnOnce(&PrimeList){thread_local!(static PRIME_LIST:UnsafeCell<PrimeList> =Default::default());PRIME_LIST.with(|cell|{unsafe{let pl=&mut*cell.get();pl.reserve(max_n);f(pl);};});}}
// codesnip-guard: Xorshift
#[cfg_attr(any(),rustfmt::skip)]pub use self::xorshift::Xorshift;#[cfg_attr(any(),rustfmt::skip)]mod xorshift{use std::{collections::hash_map::RandomState,hash::{BuildHasher,Hasher},iter::repeat_with};#[derive(Clone,Debug)]pub struct Xorshift{y:u64}impl Xorshift{pub fn new_with_seed(seed:u64)->Self{Xorshift{y:seed}}pub fn new()->Self{Xorshift::new_with_seed(RandomState::new().build_hasher().finish())}pub fn rand64(&mut self)->u64{self.y^=self.y<<5;self.y^=self.y>>17;self.y^=self.y<<11;self.y}pub fn rand(&mut self,k:u64)->u64{self.rand64()%k}pub fn rands(&mut self,k:u64,n:usize)->Vec<u64>{repeat_with(| |self.rand(k)).take(n).collect()}pub fn randf(&mut self)->f64{const UPPER_MASK:u64=0x3FF0_0000_0000_0000;const LOWER_MASK:u64=0x000F_FFFF_FFFF_FFFF;let x=self.rand64();let tmp=UPPER_MASK|(x&LOWER_MASK);let result:f64=f64::from_bits(tmp);f64::from_bits(f64::to_bits(result-1.0)^(x>>63))}pub fn gen_bool(&mut self,p:f64)->bool{self.randf()<p}pub fn shuffle<T>(&mut self,slice:&mut[T]){let mut n=slice.len();while n>1{let i=self.rand(n as _)as usize;n-=1;slice.swap(i,n);}}}impl Default for Xorshift{fn default()->Self{Xorshift::new_with_seed(0x2b99_2ddf_a232_49d6)}}}
// codesnip-guard: extgcd
#[cfg_attr(any(),rustfmt::skip)]pub fn extgcd(mut a:i64,mut b:i64)->(i64,i64,i64){let(mut u,mut v,mut x,mut y)=(1,0,0,1);while a!=0{let k=b/a;x-=k*u;y-=k*v;b-=k*a;std::mem::swap(&mut x,&mut u);std::mem::swap(&mut y,&mut v);std::mem::swap(&mut b,&mut a);}(b,x,y)}
// codesnip-guard: gcd
#[cfg_attr(any(),rustfmt::skip)]#[doc=" binary gcd"]pub fn gcd(mut a:u64,mut b:u64)->u64{if a==0{return b;}if b==0{return a;}let u=a.trailing_zeros();let v=b.trailing_zeros();a>>=u;b>>=v;let k=std::cmp::min(u,v);while a!=b{if a<b{std::mem::swap(&mut a,&mut b);}a-=b;a>>=a.trailing_zeros();}a<<k}
// codesnip-guard: lcm
#[cfg_attr(any(),rustfmt::skip)]pub fn lcm(a:u64,b:u64)->u64{if a==0&&b==0{return 0;}a/gcd(a,b)*b}
// codesnip-guard: miller_rabin
#[cfg_attr(any(),rustfmt::skip)]pub use self::miller_rabin::{miller_rabin,miller_rabin_with_br};#[cfg_attr(any(),rustfmt::skip)]mod miller_rabin{use super::BarrettReduction;macro_rules!impl_test_mr{($name:ident,$ty:ty,$upty:ty)=>{fn$name(n:$ty,br:&BarrettReduction<$upty>,a:$ty)->bool{if br.rem(a as$upty)==0{return true;}let d=n-1;let k=d.trailing_zeros();let mut d=d>>k;let mut y={let mut a=a as$upty;let mut y:$upty=1;while d>0{if d&1==1{y=br.rem(y*a);}a=br.rem(a*a);d>>=1;}y as$ty};if y==1||y==n-1{true}else{for _ in 0..k-1{y=br.rem(y as$upty*y as$upty)as$ty;if y==n-1{return true;}}false}}};}impl_test_mr!(test_mr32,u32,u64);impl_test_mr!(test_mr64,u64,u128);#[doc=" <http://miller-rabin.appspot.com/>"]macro_rules!impl_mr{($name:ident,$test:ident,$ty:ty,$upty:ty,[$($th:expr,[$($a:expr),+]),+],|$n:ident,$br:ident|$last:expr)=>{fn$name($n:$ty,$br:&BarrettReduction<$upty>)->bool{$(if$n>=$th{return$($test($n,$br,$a))&&+})+$last}};}impl_mr!(mr32,test_mr32,u32,u64,[316349281,[2,7,61],49141,[11000544,31481107]],|n,br|test_mr32(n,br,921211727));impl_mr!(mr64,test_mr64,u64,u128,[585226005592931977,[2,325,9375,28178,450775,9780504,1795265022],7999252175582851,[2,123635709730000,9233062284813009,43835965440333360,761179012939631437,1263739024124850375],55245642489451,[2,4130806001517,149795463772692060,186635894390467037,3967304179347715805],350269456337,[2,141889084524735,1199124725622454117,11096072698276303650],1050535501,[4230279247111683200,14694767155120705706,16641139526367750375]],|n,br|mr32(n as u32,&BarrettReduction::<u64>::new(n)));pub fn miller_rabin_with_br(n:u64,br:&BarrettReduction<u128>)->bool{if n.is_multiple_of(2){return n==2;}if n.is_multiple_of(3){return n==3;}if n.is_multiple_of(5){return n==5;}if n.is_multiple_of(7){return n==7;}if n<121{n>2}else{mr64(n,br)}}pub fn miller_rabin(n:u64)->bool{n>=2&&miller_rabin_with_br(n,&BarrettReduction::<u128>::new(n as u128))}}
// codesnip-guard: modinv
#[cfg_attr(any(),rustfmt::skip)]#[doc=" Requires 0 < a < m < 2^63 and gcd(a, m) == 1."]pub fn modinv(a:u64,m:u64)->u64{let(mut a,mut b)=(a as i64,m as i64);let(mut u,mut x)=(1,0);while a!=0{let k=b/a;x-=k*u;b-=k*a;std::mem::swap(&mut x,&mut u);std::mem::swap(&mut b,&mut a);}(if x<0{x+m as i64}else{x})as _}
// codesnip-guard: prime_factors
#[cfg_attr(any(),rustfmt::skip)]pub use self::prime_factors::{divisors,prime_factors,prime_factors_flatten};#[cfg_attr(any(),rustfmt::skip)]mod prime_factors{use super::{BarrettReduction,Xorshift,gcd,miller_rabin_with_br};struct MontgomeryReduction64{modulus:u64,inverse:u64}impl MontgomeryReduction64{fn new(modulus:u64)->Self{let mut inverse=modulus;for _ in 0..6{inverse=inverse.wrapping_mul(2u64.wrapping_sub(modulus.wrapping_mul(inverse)));}Self{modulus,inverse}}fn sub(&self,lhs:u64,rhs:u64)->u64{let(value,borrow)=lhs.overflowing_sub(rhs);value.wrapping_add((borrow as u64).wrapping_neg()&self.modulus)}fn mul(&self,lhs:u64,rhs:u64)->u64{let product=lhs as u128*rhs as u128;let(value,borrow)=((product>>64)as u64).overflowing_sub((((product as u64).wrapping_mul(self.inverse)as u128*self.modulus as u128)>>64)as u64);value.wrapping_add((borrow as u64).wrapping_neg()&self.modulus)}}fn find_factor(n:u64)->Option<u64>{let br=BarrettReduction::<u128>::new(n as u128);if miller_rabin_with_br(n,&br){return None;}let mr=MontgomeryReduction64::new(n);let mut rng=Xorshift::default();let(mut y0,mut c)=(0,n-1);loop{let(mut x,mut y,mut ys,mut g,mut q,mut r,mut k)=(0,y0,0,1,1,1,0);while g==1&&r<=1<<20{x=y;while k<r&&g==1{ys=y;for _ in 0..1024 .min(r-k){y=mr.sub(mr.mul(y,y),c);q=mr.mul(q,mr.sub(x,y));}g=gcd(q,n);k+=1024;}k=r;r<<=1;}if g==n{g=1;y=ys;while g==1{y=mr.sub(mr.mul(y,y),c);g=gcd(mr.sub(x,y),n);}}if g!=1&&g!=n{return Some(g);}y0=((rng.rand64()as u128*(n-2)as u128)>>64)as u64+2;c=((rng.rand64()as u128*(n-1)as u128)>>64)as u64+1;}}pub fn prime_factors_flatten(mut n:u64)->Vec<u64>{if n==0{return vec![];}let k=n.trailing_zeros();let mut res=vec![2;k as usize];n>>=k;while n.is_multiple_of(3){res.push(3);n/=3;}if n!=1{let mut c=vec![n];while let Some(n)=c.pop(){if let Some(m)=find_factor(n){c.push(m);c.push(n/m);}else{res.push(n);}}}res.sort_unstable();res}pub fn prime_factors(n:u64)->Vec<(u64,u32)>{let mut res=Vec::new();for a in prime_factors_flatten(n){if let Some((p,len))=res.last_mut()&&p==&a{*len+=1;continue;}res.push((a,1));}res}pub fn divisors(n:u64)->Vec<u64>{let mut d=vec![1u64];for(p,c)in prime_factors(n){let k=d.len();let mut acc=1;for _ in 0..c{acc*=p;for i in 0..k{d.push(d[i]*acc);}}}d.sort_unstable();d}}
// codesnip-guard: primitive_root
#[cfg_attr(any(),rustfmt::skip)]pub use self::primitive_root::{check_primitive_root,primitive_root};#[cfg_attr(any(),rustfmt::skip)]mod primitive_root{use super::{BarrettReduction,Xorshift,prime_factors};pub fn primitive_root(p:u64)->u64{if p==2{return 1;}let phi=p-1;let pf=prime_factors(phi);let br=BarrettReduction::<u128>::new(p as _);for g in 2..=3.min(p-1){if check_primitive_root(g,phi,&br,&pf){return g;}}let mut rng=Xorshift::default();loop{let g=((rng.rand64()as u128*(p-2)as u128)>>64)as u64+2;if check_primitive_root(g,phi,&br,&pf){return g;}}}pub fn check_primitive_root(g:u64,phi:u64,br:&BarrettReduction<u128>,pf:&[(u64,u32)])->bool{pf.iter().all(|&(q,_)|{let mut g=g as u128;let mut k=phi/q;let mut r:u128=1;while k>0{if k&1==1{r=br.rem(r*g);}g=br.rem(g*g);k>>=1;}r>1})}}
// 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())});}