implement zeroes

This commit is contained in:
Quaternions 2023-10-11 21:39:23 -07:00
parent a8f82a14a9
commit 69712847e3
2 changed files with 39 additions and 14 deletions

View File

@ -1,5 +1,5 @@
//integer units //integer units
#[derive(Clone,Copy,Hash,Debug)] #[derive(Clone,Copy,Hash,PartialEq,PartialOrd,Debug)]
pub struct Time(i64); pub struct Time(i64);
impl Time{ impl Time{
pub const ZERO:Self=Self(0); pub const ZERO:Self=Self(0);
@ -228,13 +228,20 @@ pub struct Unit64Mat3{
///[-1.0,1.0] = [-2^32,2^32] ///[-1.0,1.0] = [-2^32,2^32]
#[derive(Clone,Copy,Hash)] #[derive(Clone,Copy,Hash,PartialEq,PartialOrd)]
pub struct Planar64(i64); pub struct Planar64(i64);
impl Planar64{ impl Planar64{
pub const ZERO:Self=Self(0);
pub const ONE:Self=Self(2<<32); pub const ONE:Self=Self(2<<32);
pub fn int(num:i32)->Self{ pub fn int(num:i32)->Self{
Self(Self::ONE.0*num as i64) Self(Self::ONE.0*num as i64)
} }
pub fn raw(num:i64)->Self{
Self(num)
}
pub fn get(&self)->i64{
self.0
}
pub fn from_ratio(num:i64,den:std::num::NonZeroU64)->Self{ pub fn from_ratio(num:i64,den:std::num::NonZeroU64)->Self{
Self(Self::ONE.0*num/den.get() as i64) Self(Self::ONE.0*num/den.get() as i64)
} }
@ -244,6 +251,13 @@ impl Into<f32> for Planar64{
self.0 as f32/(2<<32) as f32 self.0 as f32/(2<<32) as f32
} }
} }
impl std::ops::Neg for Planar64{
type Output=Planar64;
#[inline]
fn neg(self)->Self::Output{
Planar64(-self.0)
}
}
impl std::ops::Add<Planar64> for Planar64{ impl std::ops::Add<Planar64> for Planar64{
type Output=Planar64; type Output=Planar64;
#[inline] #[inline]
@ -251,6 +265,13 @@ impl std::ops::Add<Planar64> for Planar64{
Planar64(self.0+rhs.0) Planar64(self.0+rhs.0)
} }
} }
impl std::ops::Sub<Planar64> for Planar64{
type Output=Planar64;
#[inline]
fn sub(self, rhs: Self) -> Self::Output {
Planar64(self.0-rhs.0)
}
}
impl std::ops::Mul<i64> for Planar64{ impl std::ops::Mul<i64> for Planar64{
type Output=Planar64; type Output=Planar64;
#[inline] #[inline]

View File

@ -1,26 +1,30 @@
//find roots of polynomials //find roots of polynomials
use crate::integer::Planar64;
#[inline] #[inline]
pub fn zeroes2(a0:f32,a1:f32,a2:f32) -> Vec<f32>{ pub fn zeroes2(a0:Planar64,a1:Planar64,a2:Planar64) -> Vec<Planar64>{
if a2==0f32{ if a2==Planar64::ZERO{
return zeroes1(a0, a1); return zeroes1(a0, a1);
} }
let mut radicand=a1*a1-4f32*a2*a0; let mut radicand=a1.get() as i128*a1.get() as i128-a2.get() as i128*a0.get() as i128*4;
if 0f32<radicand { if 0<radicand {
radicand=radicand.sqrt(); //start with f64 sqrt
if 0f32<a2 { let planar_radicand=Planar64::raw(unsafe{(radicand as f64).sqrt().to_int_unchecked()});
return vec![(-a1-radicand)/(2f32*a2),(-a1+radicand)/(2f32*a2)]; //TODO: one or two newtons
if Planar64::ZERO<a2 {
return vec![(-a1-planar_radicand)/(a2*2),(-a1+planar_radicand)/(a2*2)];
} else { } else {
return vec![(-a1+radicand)/(2f32*a2),(-a1-radicand)/(2f32*a2)]; return vec![(-a1+planar_radicand)/(a2*2),(-a1-planar_radicand)/(a2*2)];
} }
} else if radicand==0f32 { } else if radicand==0 {
return vec![-a1/(2f32*a2)]; return vec![a1/(a2*-2)];
} else { } else {
return vec![]; return vec![];
} }
} }
#[inline] #[inline]
pub fn zeroes1(a0:f32,a1:f32) -> Vec<f32> { pub fn zeroes1(a0:Planar64,a1:Planar64) -> Vec<Planar64> {
if a1==0f32{ if a1==Planar64::ZERO{
return vec![]; return vec![];
} else { } else {
return vec![-a0/a1]; return vec![-a0/a1];