Unfinished first draft of division

Still to do: checking for division by 0, and speeding it up using a
dedicated square() method.

diff --git a/src/finitenimber.rs b/src/finitenimber.rs
index 42c67a97f6cb3c66a0c213ced24164e06e75a7f3..c74e3a1d98763896cd5d8d57ae10a3741b41abe4 100644 (file)
--- a/src/finitenimber.rs
+++ b/src/finitenimber.rs
@@ -1,7 +1,7 @@
 use core::borrow::Borrow;
 use core::cmp::max;
 use core::fmt::{Debug, Display, Formatter};
-use core::ops::{Add, Mul, Sub};
+use core::ops::{Add, Div, Mul, Sub};
 use itertools::Itertools;
 
 pub type Word = u64; // element type of the vectors we use
@@ -388,9 +388,47 @@ impl<'a, 'b> Mul<FiniteNimberRef<'a>> for FiniteNimberRef<'b> {
     }
 }
 
+impl<'a> FiniteNimberRef<'a> {
+    fn inverse_recurse(self, level: usize) -> FiniteNimber {
+        match level.checked_sub(1) {
+            Some(sublevel) => {
+                let (alo, ahi) = self.split(sublevel);
+                let sum = alo + ahi;
+                let sq = ahi * ahi; // todo! dedicated squaring function
+                let det = alo * sum.to_ref() + sq.to_ref().mul_by_h(sublevel);
+                let detinv = det.to_ref().inverse_recurse(sublevel);
+                (detinv.to_ref() * sum)
+                    .to_ref()
+                    .join((detinv.to_ref() * ahi).to_ref(), sublevel)
+            }
+
+            // At level 0, we're in GF(2), so the only invertible
+            // number is 1, whose inverse is the same as itself.
+            //
+            // We also map 0 to 0 in the base case. For some
+            // applications this makes vague sense (in particular,
+            // within a given subfield of order q, raising to the
+            // power q-2 is an extension of inversion which maps 0 to
+            // 0). The wrapping implementation of the Div trait will
+            // check for division by zero.
+            None => FiniteNimber::from(self.low_word()),
+        }
+    }
+}
+
+impl<'a, 'b> Div<FiniteNimberRef<'a>> for FiniteNimberRef<'b> {
+    type Output = FiniteNimber;
+    fn div(self, other: FiniteNimberRef<'a>) -> FiniteNimber {
+        // FIXME: how _do_ we check for division by zero?
+        let level = max(self.level(), other.level());
+        self.mul_recurse(other.inverse_recurse(level).to_ref(), level)
+    }
+}
+
 impl_binop_wrappers!(Add, add, AddAssign, add_assign);
 impl_binop_wrappers!(Sub, sub, SubAssign, sub_assign);
 impl_binop_wrappers!(Mul, mul, MulAssign, mul_assign);
+impl_binop_wrappers!(Div, div, DivAssign, div_assign);
 
 #[cfg(test)]
 mod tests {
@@ -663,4 +701,44 @@ mod tests {
             ])
         );
     }
+
+    #[test]
+    fn div() {
+        assert_eq!(
+            FiniteNimber::from(1) / FiniteNimber::from(2),
+            FiniteNimber::from(3)
+        );
+
+        assert_eq!(
+            FiniteNimber::from(vec![
+                0x72890f944be121d8,
+                0xe6429b02014feb2e,
+                0x2454070e4408eff8,
+                0x298c025ec5ac4190,
+                0xba6895a109cfcf6d,
+                0xcfb08c013e9dd19c,
+                0x5eeb02eaf7ae9ea0,
+                0x59cbe194a9599171,
+            ]) / FiniteNimber::from(vec![
+                0x3f84d5b5b5470917,
+                0xc0ac29b7c97c50dd,
+                0xbe5466cf34e90c6c,
+                0x452821e638d01377,
+                0x082efa98ec4e6c89,
+                0xa4093822299f31d0,
+                0x13198a2e03707344,
+                0x243f6a8885a308d3,
+            ]),
+            FiniteNimber::from(vec![
+                0x4f7c7b5757f59584,
+                0xda06c80abb1185eb,
+                0xf4bf8d8d8c31d763,
+                0x324e7738926cfbe5,
+                0xa784d9045190cfef,
+                0x62e7160f38b4da56,
+                0xbf7158809cf4f3c7,
+                0xb7e151628aed2a6a,
+            ])
+        );
+    }
 }