From dc7af1b333f3cef4014cb612d5a63b3e9321a81c Mon Sep 17 00:00:00 2001
From: Josh Stone <cuviper@gmail.com>
Date: Thu, 21 Feb 2019 16:06:55 -0800
Subject: [PATCH] Add `finv()` and `fdiv()`

These use the floating point `norm()` to avoid cases where the more
generic `inv()` and `div()` might overflow in `norm_sqr()`.

Closes #23.
---
 src/lib.rs | 56 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 56 insertions(+)

diff --git a/src/lib.rs b/src/lib.rs
index 69461fd..b3b3fb4 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -408,6 +408,62 @@ impl<T: Clone + Float> Complex<T> {
         }
         ((one + self).ln() - (one - self).ln()) / two
     }
+
+    /// Returns `1/self` using floating-point operations.
+    ///
+    /// This may be more accurate than the generic `self.inv()` in cases
+    /// where `self.norm_sqr()` would overflow to ∞ or underflow to 0.
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// use num_complex::Complex64;
+    /// let c = Complex64::new(1e300, 1e300);
+    ///
+    /// // The generic `inv()` will overflow.
+    /// assert!(!c.inv().is_normal());
+    ///
+    /// // But we can do better for `Float` types.
+    /// let inv = c.finv();
+    /// assert!(inv.is_normal());
+    /// println!("{:e}", inv);
+    ///
+    /// let expected = Complex64::new(5e-301, -5e-301);
+    /// assert!((inv - expected).norm() < 1e-315);
+    /// ```
+    #[inline]
+    pub fn finv(&self) -> Complex<T> {
+        let norm = self.norm();
+        self.conj() / norm / norm
+    }
+
+    /// Returns `self/other` using floating-point operations.
+    ///
+    /// This may be more accurate than the generic `Div` implementation in cases
+    /// where `other.norm_sqr()` would overflow to ∞ or underflow to 0.
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// use num_complex::Complex64;
+    /// let a = Complex64::new(2.0, 3.0);
+    /// let b = Complex64::new(1e300, 1e300);
+    ///
+    /// // Generic division will overflow.
+    /// assert!(!(a / b).is_normal());
+    ///
+    /// // But we can do better for `Float` types.
+    /// let quotient = a.fdiv(b);
+    /// assert!(quotient.is_normal());
+    /// println!("{:e}", quotient);
+    ///
+    /// let expected = Complex64::new(2.5e-300, 5e-301);
+    /// assert!((quotient - expected).norm() < 1e-315);
+    /// ```
+    #[inline]
+    pub fn fdiv(&self, other: Complex<T>) -> Complex<T> {
+        self * other.finv()
+    }
 }
 
 impl<T: Clone + FloatCore> Complex<T> {