2023-09-18 20:20:34 +00:00
|
|
|
//find roots of polynomials
|
2023-09-27 09:12:20 +00:00
|
|
|
use crate::integer::Planar64;
|
|
|
|
|
2023-09-22 22:21:59 +00:00
|
|
|
#[inline]
|
2023-09-27 09:12:20 +00:00
|
|
|
pub fn zeroes2(a0:Planar64,a1:Planar64,a2:Planar64) -> Vec<Planar64>{
|
|
|
|
if a2==Planar64::ZERO{
|
2023-09-18 20:20:34 +00:00
|
|
|
return zeroes1(a0, a1);
|
|
|
|
}
|
2023-10-18 23:30:02 +00:00
|
|
|
let radicand=a1.get() as i128*a1.get() as i128-a2.get() as i128*a0.get() as i128*4;
|
2023-09-27 09:12:20 +00:00
|
|
|
if 0<radicand {
|
|
|
|
//start with f64 sqrt
|
|
|
|
let planar_radicand=Planar64::raw(unsafe{(radicand as f64).sqrt().to_int_unchecked()});
|
|
|
|
//TODO: one or two newtons
|
2023-11-30 05:21:10 +00:00
|
|
|
//sort roots ascending and avoid taking the difference of large numbers
|
|
|
|
match (Planar64::ZERO<a2,Planar64::ZERO<a1){
|
|
|
|
(true, true )=>vec![(-a1-planar_radicand)/(a2*2),(a0*2)/(-a1-planar_radicand)],
|
|
|
|
(true, false)=>vec![(a0*2)/(-a1+planar_radicand),(-a1+planar_radicand)/(a2*2)],
|
|
|
|
(false,true )=>vec![(a0*2)/(-a1-planar_radicand),(-a1-planar_radicand)/(a2*2)],
|
|
|
|
(false,false)=>vec![(-a1+planar_radicand)/(a2*2),(a0*2)/(-a1+planar_radicand)],
|
2023-09-18 20:20:34 +00:00
|
|
|
}
|
2023-09-27 09:12:20 +00:00
|
|
|
} else if radicand==0 {
|
|
|
|
return vec![a1/(a2*-2)];
|
2023-09-18 20:20:34 +00:00
|
|
|
} else {
|
|
|
|
return vec![];
|
|
|
|
}
|
|
|
|
}
|
|
|
|
#[inline]
|
2023-09-27 09:12:20 +00:00
|
|
|
pub fn zeroes1(a0:Planar64,a1:Planar64) -> Vec<Planar64> {
|
|
|
|
if a1==Planar64::ZERO{
|
2023-09-18 20:20:34 +00:00
|
|
|
return vec![];
|
|
|
|
} else {
|
|
|
|
return vec![-a0/a1];
|
|
|
|
}
|
|
|
|
}
|