sqrt: closed loop over bit shift

This commit is contained in:
Quaternions 2024-08-29 18:29:04 -07:00
parent cd1aa26293
commit 61aad93f8d

View File

@ -69,30 +69,18 @@ impl<const CHUNKS:usize,Frac:Unsigned> Fixed<CHUNKS,Frac>
//sqrt
//0110.0000
//pow2 = 0100.0000
let mut pow2=Self{
bits:BInt::<CHUNKS>::ONE.shl(((((CHUNKS as i32*64-Frac::I32-(self.bits.leading_zeros() as i32)+1)>>1)+Frac::I32) as u32).saturating_sub(1)),
frac:PhantomData,
};
let max_shift=((((CHUNKS as i32*64-Frac::I32-(self.bits.leading_zeros() as i32)+1)>>1)+Frac::I32) as u32).saturating_sub(1);
let mut result=Self::ZERO;
//cheat to make the types match
let wide_self=self.wide_mul(Fixed::<CHUNKS,Frac>::ONE);
loop{
//TODO: closed loop over bit shift exponent rather than pow2
if pow2==Self::ZERO{
break result;
for shift in (0..=max_shift).rev(){
let new_result=result|(Fixed::<CHUNKS,Frac>::EPSILON<<shift);
if new_result.wide_mul(new_result)<=wide_self{
result=new_result;
}
//TODO: flip a single bit instead of adding a power of 2
let new_result=result+pow2;
//note that the implicit truncation in the multiply
//means that the algorithm can return a result which squares to a number greater than the input.
match wide_self.cmp(&new_result.wide_mul(new_result)){
core::cmp::Ordering::Less=>(),
core::cmp::Ordering::Equal=>break new_result,
core::cmp::Ordering::Greater=>result=new_result,
}
pow2>>=1;
}
result
}
pub fn sqrt(self)->Self{
if self<Self::ZERO{