StrafesNET/common
3
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

ruin zeroes

Browse Source
This commit is contained in:
Quaternions 2024-08-27 16:34:21 -07:00
parent 76b1e6af54
commit 3bb280a787
1 changed files with 37 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

68
src/zeroes.rs
View File

@ -1,41 +1,47 @@
//find roots of polynomials //find roots of polynomials
use arrayvec::ArrayVec; use arrayvec::ArrayVec;
use crate::integer::Planar64; 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] #[inline]
pub fn zeroes2(a0:Planar64,a1:Planar64,a2:Planar64)->ArrayVec<Planar64,2>{ 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>
if a2==Planar64::ZERO{ where
return zeroes1(a0,a1); Fixed<N,F>:WideMul,
} <Fixed<N,F> as WideMul>::Output:std::ops::Mul<i64>,
let radicand=a1.get() as i128*a1.get() as i128-a2.get() as i128*a0.get() as i128*4; <Fixed<N,F> as WideMul>::Output:std::ops::Sub<<<Fixed<N,F> as WideMul>::Output as std::ops::Mul<i64>>::Output>
if 0<radicand{ {
//start with f64 sqrt let a2pos=match a2.cmp(&Fixed::<N,F>::ZERO){
//failure case: 2^63 < sqrt(2^127) Ordering::Greater=>true,
let planar_radicand=Planar64::raw(unsafe{(radicand as f64).sqrt().to_int_unchecked()}); Ordering::Equal=>return zeroes1(a0,a1),
//TODO: one or two newtons Ordering::Less=>true,
//sort roots ascending and avoid taking the difference of large numbers };
match (Planar64::ZERO<a2,Planar64::ZERO<a1){ let radicand=a1.wide_mul(a1)-a2.wide_mul(a0)*4;
(true, true )=>[(-a1-planar_radicand)/(a2*2),(a0*2)/(-a1-planar_radicand)].into(), match radicand.cmp(&Fixed::<N,F>::ZERO){
(true, false)=>[(a0*2)/(-a1+planar_radicand),(-a1+planar_radicand)/(a2*2)].into(), Ordering::Greater=>{
(false,true )=>[(a0*2)/(-a1-planar_radicand),(-a1-planar_radicand)/(a2*2)].into(), //start with f64 sqrt
(false,false)=>[(-a1+planar_radicand)/(a2*2),(a0*2)/(-a1+planar_radicand)].into(), //failure case: 2^63 < sqrt(2^127)
} let planar_radicand=radicand.sqrt();
}else if radicand==0{ //TODO: one or two newtons
return ArrayVec::from_iter([a1/(a2*-2)]); //sort roots ascending and avoid taking the difference of large numbers
}else{ match (a2pos,Fixed::<N,F>::ZERO<a1){
return ArrayVec::new_const(); (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] #[inline]
pub fn zeroes1(a0:Planar64,a1:Planar64)->ArrayVec<Planar64,2>{ 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==Planar64::ZERO{ if a1==Fixed::<N,F>::ZERO{
return ArrayVec::new_const(); ArrayVec::new_const()
}else{ }else{
let q=((-a0.get() as i128)<<32)/(a1.get() as i128); ArrayVec::from_iter([Ratio::new(-a0,a1)])
if i64::MIN as i128<=q&&q<=i64::MAX as i128{
return ArrayVec::from_iter([Planar64::raw(q as i64)]);
}else{
return ArrayVec::new_const();
}
} }
} }