From 5234795895ae475efe637fd193158c0b0ef05890 Mon Sep 17 00:00:00 2001
From: Leo Lahti <leo@lahti.pw>
Date: Tue, 3 Sep 2024 15:01:32 +0300
Subject: [PATCH] Moved division to ComplexField and implemented Div for
 missing types.

More specificially, implemented Div(Assign)<Scale<RhsE>> for
Scale<LhsE>.
---
 faer-entity/src/lib.rs         | 48 ++++++++++++++++++++++------------
 src/complex_native/c32_impl.rs |  8 ++++++
 src/complex_native/c64_impl.rs |  7 +++++
 src/linalg/mat_ops.rs          | 17 ++++++++++++
 src/sparse/linalg/mod.rs       |  9 +++----
 5 files changed, 68 insertions(+), 21 deletions(-)

diff --git a/faer-entity/src/lib.rs b/faer-entity/src/lib.rs
index 8de486c4..bed0096a 100644
--- a/faer-entity/src/lib.rs
+++ b/faer-entity/src/lib.rs
@@ -422,16 +422,18 @@ pub trait ComplexField:
     + core::ops::Add<Self, Output = Self>
     + core::ops::Sub<Self, Output = Self>
     + core::ops::Mul<Self, Output = Self>
+    + core::ops::Div<Self, Output = Self>
     + core::ops::AddAssign<Self>
     + core::ops::SubAssign<Self>
     + core::ops::MulAssign<Self>
+    + core::ops::DivAssign<Self>
 {
     type Real: RealField;
     type Simd: SimdCtx;
     type ScalarSimd: SimdCtx;
     type PortableSimd: SimdCtx;
 
-    /// Converts `value` from `f64` to `Self`.  
+    /// Converts `value` from `f64` to `Self`.
     /// The conversion may be lossy when converting to a type with less precision.
     fn faer_from_f64(value: f64) -> Self;
 
@@ -441,6 +443,8 @@ pub trait ComplexField:
     fn faer_sub(self, rhs: Self) -> Self;
     /// Returns `self * rhs`.
     fn faer_mul(self, rhs: Self) -> Self;
+    /// Returns `self / rhs`.
+    fn faer_div(self, rhs: Self) -> Self;
 
     /// Returns `-self`.
     fn faer_neg(self) -> Self;
@@ -706,8 +710,6 @@ pub trait RealField:
     fn faer_min_positive_sqrt() -> Self;
     fn faer_min_positive_sqrt_inv() -> Self;
 
-    fn faer_div(self, rhs: Self) -> Self;
-
     fn faer_usize_to_index(a: usize) -> IndexFor<Self>;
     fn faer_index_to_usize(a: IndexFor<Self>) -> usize;
     fn faer_max_index() -> IndexFor<Self>;
@@ -788,6 +790,11 @@ impl ComplexField for f32 {
         self * rhs
     }
 
+    #[inline(always)]
+    fn faer_div(self, rhs: Self) -> Self {
+        self / rhs
+    }
+
     #[inline(always)]
     fn faer_neg(self) -> Self {
         -self
@@ -1117,6 +1124,11 @@ impl ComplexField for f64 {
         self * rhs
     }
 
+    #[inline(always)]
+    fn faer_div(self, rhs: Self) -> Self {
+        self / rhs
+    }
+
     #[inline(always)]
     fn faer_neg(self) -> Self {
         -self
@@ -1429,11 +1441,6 @@ impl RealField for f32 {
         Self::MIN_POSITIVE
     }
 
-    #[inline(always)]
-    fn faer_div(self, rhs: Self) -> Self {
-        self / rhs
-    }
-
     #[inline(always)]
     fn faer_usize_to_index(a: usize) -> IndexFor<Self> {
         a as _
@@ -1563,10 +1570,6 @@ impl RealField for f64 {
     fn faer_zero_threshold() -> Self {
         Self::MIN_POSITIVE
     }
-    #[inline(always)]
-    fn faer_div(self, rhs: Self) -> Self {
-        self / rhs
-    }
 
     #[inline(always)]
     fn faer_usize_to_index(a: usize) -> IndexFor<Self> {
@@ -2159,6 +2162,18 @@ impl<E: RealField> ComplexField for Complex<E> {
         }
     }
 
+    #[inline(always)]
+    fn faer_div(self, rhs: Self) -> Self {
+        Self {
+            // (self.re * rhs.re + self.im * rhs.im)/(rhs.re^2 + rhs.im^2)
+            re: Self::Real::faer_div(Self::Real::faer_add(self.re.faer_mul(rhs.re), self.im.faer_mul(rhs.im)),
+                rhs.faer_abs2()),
+            // (self.im * rhs.re - self.re * rhs.im)/(rhs.re^2 + rhs.im^2)
+            im: Self::Real::faer_div(Self::Real::faer_sub(self.im.faer_mul(rhs.re), self.re.faer_mul(rhs.im)),
+                rhs.faer_abs2()),
+        }
+    }
+
     #[inline(always)]
     fn faer_neg(self) -> Self {
         Self {
@@ -2912,10 +2927,6 @@ impl RealField for Symbolic {
         Self
     }
 
-    #[inline(always)]
-    fn faer_div(self, _rhs: Self) -> Self {
-        Self
-    }
 
     #[inline(always)]
     fn faer_usize_to_index(a: usize) -> Self::Index {
@@ -3069,6 +3080,11 @@ impl ComplexField for Symbolic {
         Self
     }
 
+    #[inline(always)]
+    fn faer_div(self, _rhs: Self) -> Self {
+        Self
+    }
+
     #[inline(always)]
     fn faer_neg(self) -> Self {
         Self
diff --git a/src/complex_native/c32_impl.rs b/src/complex_native/c32_impl.rs
index 4f7e17b6..7013eaa5 100644
--- a/src/complex_native/c32_impl.rs
+++ b/src/complex_native/c32_impl.rs
@@ -646,6 +646,14 @@ impl ComplexField for c32 {
         }
     }
 
+    #[inline(always)]
+    fn faer_div(self, rhs: Self) -> Self {
+        Self {
+            re: (self.re * rhs.re + self.im * rhs.im) / (rhs.re.powi(2) + rhs.im.powi(2)),
+            im: (self.im * rhs.re - self.re * rhs.im) / (rhs.re.powi(2) + rhs.im.powi(2)),
+        }
+    }
+
     #[inline(always)]
     fn faer_neg(self) -> Self {
         Self {
diff --git a/src/complex_native/c64_impl.rs b/src/complex_native/c64_impl.rs
index 2867acac..aa4945a1 100644
--- a/src/complex_native/c64_impl.rs
+++ b/src/complex_native/c64_impl.rs
@@ -646,6 +646,13 @@ impl ComplexField for c64 {
             im: self.re * rhs.im + self.im * rhs.re,
         }
     }
+    #[inline(always)]
+    fn faer_div(self, rhs: Self) -> Self {
+        Self {
+            re: (self.re * rhs.re + self.im * rhs.im) / (rhs.re.powi(2) + rhs.im.powi(2)),
+            im: (self.im * rhs.re - self.re * rhs.im) / (rhs.re.powi(2) + rhs.im.powi(2)),
+        }
+    }
 
     #[inline(always)]
     fn faer_neg(self) -> Self {
diff --git a/src/linalg/mat_ops.rs b/src/linalg/mat_ops.rs
index 3caae0e7..3b8a22c5 100644
--- a/src/linalg/mat_ops.rs
+++ b/src/linalg/mat_ops.rs
@@ -1166,6 +1166,17 @@ impl<E: ComplexField, LhsE: Conjugate<Canonical = E>, RhsE: Conjugate<Canonical
     }
 }
 
+impl<E: ComplexField, LhsE: Conjugate<Canonical = E>, RhsE: Conjugate<Canonical = E>>
+    Div<Scale<RhsE>> for Scale<LhsE>
+{
+    type Output = Scale<E>;
+
+    #[inline]
+    fn div(self, rhs: Scale<RhsE>) -> Self::Output {
+        Scale(self.0.canonicalize().faer_div(rhs.0.canonicalize()))
+    }
+}
+
 impl<LhsE: ComplexField, RhsE: Conjugate<Canonical = LhsE>> MulAssign<Scale<RhsE>> for Scale<LhsE> {
     #[inline]
     fn mul_assign(&mut self, rhs: Scale<RhsE>) {
@@ -1186,6 +1197,12 @@ impl<LhsE: ComplexField, RhsE: Conjugate<Canonical = LhsE>> SubAssign<Scale<RhsE
         self.0 = self.0.faer_sub(rhs.0.canonicalize())
     }
 }
+impl<LhsE: ComplexField, RhsE: Conjugate<Canonical = LhsE>> DivAssign<Scale<RhsE>> for Scale<LhsE> {
+    #[inline]
+    fn div_assign(&mut self, rhs: Scale<RhsE>) {
+        self.0 = self.0.faer_div(rhs.0.canonicalize())
+    }
+}
 
 impl<E: ComplexField, LhsE: Conjugate<Canonical = E>, RhsE: Conjugate<Canonical = E>>
     Mul<MatRef<'_, RhsE>> for MatRef<'_, LhsE>
diff --git a/src/sparse/linalg/mod.rs b/src/sparse/linalg/mod.rs
index 8af488e3..d48535ad 100644
--- a/src/sparse/linalg/mod.rs
+++ b/src/sparse/linalg/mod.rs
@@ -850,11 +850,6 @@ pub(crate) mod qd {
                 Self::MIN_POSITIVE
             }
 
-            #[inline(always)]
-            fn faer_div(self, rhs: Self) -> Self {
-                self / rhs
-            }
-
             #[inline(always)]
             fn faer_usize_to_index(a: usize) -> Self::Index {
                 a as _
@@ -1014,6 +1009,10 @@ pub(crate) mod qd {
             fn faer_mul(self, rhs: Self) -> Self {
                 self * rhs
             }
+            #[inline(always)]
+            fn faer_div(self, rhs: Self) -> Self {
+                self / rhs
+            }
 
             #[inline(always)]
             fn faer_neg(self) -> Self {