From 3bb280a7879bc3877f534352b1f83a36cddedb6e Mon Sep 17 00:00:00 2001
From: Quaternions <krakow20@gmail.com>
Date: Tue, 27 Aug 2024 16:34:21 -0700
Subject: [PATCH] ruin zeroes

---
 src/zeroes.rs | 68 ++++++++++++++++++++++++++++-----------------------
 1 file changed, 37 insertions(+), 31 deletions(-)

diff --git a/src/zeroes.rs b/src/zeroes.rs
index 0117963..5a16c57 100644
--- a/src/zeroes.rs
+++ b/src/zeroes.rs
@@ -1,41 +1,47 @@
 //find roots of polynomials
 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]
-pub fn zeroes2(a0:Planar64,a1:Planar64,a2:Planar64)->ArrayVec<Planar64,2>{
-	if a2==Planar64::ZERO{
-		return zeroes1(a0,a1);
-	}
-	let radicand=a1.get() as i128*a1.get() as i128-a2.get() as i128*a0.get() as i128*4;
-	if 0<radicand{
-		//start with f64 sqrt
-		//failure case: 2^63 < sqrt(2^127)
-		let planar_radicand=Planar64::raw(unsafe{(radicand as f64).sqrt().to_int_unchecked()});
-		//TODO: one or two newtons
-		//sort roots ascending and avoid taking the difference of large numbers
-		match (Planar64::ZERO<a2,Planar64::ZERO<a1){
-			(true, true )=>[(-a1-planar_radicand)/(a2*2),(a0*2)/(-a1-planar_radicand)].into(),
-			(true, false)=>[(a0*2)/(-a1+planar_radicand),(-a1+planar_radicand)/(a2*2)].into(),
-			(false,true )=>[(a0*2)/(-a1-planar_radicand),(-a1-planar_radicand)/(a2*2)].into(),
-			(false,false)=>[(-a1+planar_radicand)/(a2*2),(a0*2)/(-a1+planar_radicand)].into(),
-		}
-	}else if radicand==0{
-		return ArrayVec::from_iter([a1/(a2*-2)]);
-	}else{
-		return ArrayVec::new_const();
+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(a0:Planar64,a1:Planar64)->ArrayVec<Planar64,2>{
-	if a1==Planar64::ZERO{
-		return ArrayVec::new_const();
+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{
-		let q=((-a0.get() as i128)<<32)/(a1.get() as i128);
-		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();
-		}
+		ArrayVec::from_iter([Ratio::new(-a0,a1)])
 	}
 }