//find roots of polynomials use arrayvec::ArrayVec; use std::cmp::Ordering; use deferred_division::ratio::Ratio; use fixed_wide::typenum::Unsigned; use fixed_wide::fixed::Fixed; use fixed_wide::wide::WideMul; #[inline] pub fn zeroes2<const N:usize,F:Unsigned>(a0:Fixed<N,F>,a1:Fixed<N,F>,a2:Fixed<N,F>)->ArrayVec<Ratio<Fixed<N,F>,Fixed<N,F>>,2> where Fixed<N,F>:WideMul, <Fixed<N,F> as WideMul>::Output:std::ops::Mul<i64>, <Fixed<N,F> as WideMul>::Output:std::ops::Sub<<<Fixed<N,F> as WideMul>::Output as std::ops::Mul<i64>>::Output> { let a2pos=match a2.cmp(&Fixed::<N,F>::ZERO){ Ordering::Greater=>true, Ordering::Equal=>return zeroes1(a0,a1), Ordering::Less=>true, }; let radicand=a1.wide_mul(a1)-a2.wide_mul(a0)*4; match radicand.cmp(&Fixed::<N,F>::ZERO){ Ordering::Greater=>{ //start with f64 sqrt //failure case: 2^63 < sqrt(2^127) let planar_radicand=radicand.sqrt(); //TODO: one or two newtons //sort roots ascending and avoid taking the difference of large numbers match (a2pos,Fixed::<N,F>::ZERO<a1){ (true, true )=>[Ratio::new(-a1-planar_radicand,a2*2),Ratio::new(a0*2,-a1-planar_radicand)].into(), (true, false)=>[Ratio::new(a0*2,-a1+planar_radicand),Ratio::new(-a1+planar_radicand,a2*2)].into(), (false,true )=>[Ratio::new(a0*2,-a1-planar_radicand),Ratio::new(-a1-planar_radicand,a2*2)].into(), (false,false)=>[Ratio::new(-a1+planar_radicand,a2*2),Ratio::new(a0*2,-a1+planar_radicand)].into(), } }, Ordering::Equal=>ArrayVec::from_iter([Ratio::new(a1,a2*-2)]), Ordering::Less=>ArrayVec::new_const(), } } #[inline] pub fn zeroes1<const N:usize,F:Unsigned>(a0:Fixed<N,F>,a1:Fixed<N,F>)->ArrayVec<Ratio<Fixed<N,F>,Fixed<N,F>>,2>{ if a1==Fixed::<N,F>::ZERO{ ArrayVec::new_const() }else{ ArrayVec::from_iter([Ratio::new(-a0,a1)]) } }