/// subfield that contains it.
///
/// In other words:
-/// * the nimbers *0 and *1 have level 0
-/// * the nimbers *2 and *3 have level 1
-/// * the nimbers *4 to *15 inclusive have level 2
-/// * the nimbers *16 to *255 inclusive have level 3, etc.
+/// * the nimbers \*0 and \*1 have level 0
+/// * the nimbers \*2 and \*3 have level 1
+/// * the nimbers \*4 to \*15 inclusive have level 2
+/// * the nimbers \*16 to \*255 inclusive have level 3, etc.
///
/// Most nimber calculation algorithms (except addition, which is
/// exceptionally simple) involve recursing to lower and lower levels,
/// parameter `level` giving the (maximum) level of the input nimbers,
/// which makes it easy to determine their remaining recursion depth.
/// The top-level wrappers in `FiniteNimber` will call
-/// [`FiniteNimberRef::level`] to find the level(s) of their input
+/// [`FiniteNimberRef::level()`] to find the level(s) of their input
/// nimbers, and use that to start the recursion.
///
/// In the descriptions of these recursive algorithms, we'll refer to
/// [`FiniteNimberRef::split()`] writes a nimber `x` in the form
/// `xhi*t+xlo` and returns the tuple `(xlo, xhi)` simply by
/// extracting two substrings of the bits of `x`, and
-/// [`FiniteNimberRef.join()`] performs the reverse operation equally
+/// [`FiniteNimberRef::join()`] performs the reverse operation equally
/// simply.
///
/// Another important constant is `h`, which refers to the nimber
/// usually be performed by multiplying out those polynomials in t,
/// and then reducing to degree < 2 by replacing t^2 with t+h wherever
/// it appears. Multiplication by h is therefore needed a lot: it's
-/// done by a special method [`FiniteNimberRef::mul_by_h`].
+/// done by a special method [`FiniteNimberRef::mul_by_h()`].
///
/// For levels above 6, a nimber looks like a list of `Word`, and a
/// `FiniteNimberRef` stores a `&[Word]`, so that `split()` simply
/// directly.
#[derive(Clone, Copy)]
-pub enum FiniteNimberRef<'a> {
+enum FiniteNimberRef<'a> {
/// A nimber of level 6 or below, that is, with an integer value
/// less than 2^64, is stored in this branch by directly storing
/// that integer.
}
}
+/// Internal helper methods.
impl FiniteNimber {
+ /// Return a [`FiniteNimberRef`] referring to the contents of this
+ /// object.
fn to_ref(&self) -> FiniteNimberRef {
match &self.0 {
FiniteNimberEnum::Single(w) => FiniteNimberRef::Single(*w),
},
}
}
+
+ /// Test a specific bit in the integer representation of this
+ /// nimber.
+ fn test_bit(&self, bit: usize) -> bool {
+ let word = bit >> WORDLEVELS;
+ let shift = (bit as u32) & (Word::BITS - 1);
+ ((self.as_slice().get(word).cloned().unwrap_or(0) >> shift) & 1) != 0
+ }
+
+ /// Sort two nimbers into order of their representing integers.
+ fn sort_pair(a: Self, b: Self) -> (Self, Self) {
+ let misordered = a
+ .as_slice()
+ .iter()
+ .cloned()
+ .zip_longest(b.as_slice().iter().cloned())
+ .rev()
+ .map(|pair| {
+ let (left, right) = pair.or(0, 0);
+ if left != right {
+ Some(left > right)
+ } else {
+ None
+ }
+ })
+ .find_map(|v| v)
+ .unwrap_or(false);
+
+ if misordered {
+ (b, a)
+ } else {
+ (a, b)
+ }
+ }
}
impl<'a> From<&'a FiniteNimber> for FiniteNimberRef<'a> {
/// (ah t + ah + al) / (al (ah + al) + ah^2 h)
/// ```
///
- /// where the denominator has no t in it, so it can be done by
- /// inversion at the next smaller level.
+ /// This still involves a division, but the denominator has no t
+ /// in it, so finding its inverse is an operation at the next
+ /// smaller level: we've halved the size of the division required.
fn inverse_recurse(self, level: usize) -> Option<FiniteNimber> {
match level.checked_sub(1) {
Some(sublevel) => {
}
}
+/// Internal functions having to do with solving quadratic equations.
+impl<'a> FiniteNimberRef<'a> {
+ /// Recursively calculate a solution to a quadratic equation in
+ /// the normalised form x^2 + x + c, for some constant c. The
+ /// constant c is the receiver of the method, that is, you call
+ /// `c.quadratic1_recurse(level)` and get back x.
+ ///
+ /// It's unspecified which of the two solutions to the quadratic
+ /// you get. The other one is obtainable as x+1.
+ ///
+ /// Within a particular subfield of the finite nimbers, the
+ /// expression x^2+x always delivers an answer with the high bit
+ /// clear (easily proved by induction). The caller of this
+ /// function is responsible for ensuring that the input nimber
+ /// does have its high bit clear. That is, bit (2^level-1) should
+ /// be zero. E.g. at level=3, `self` should be a nimber in the
+ /// range *0 to *0x7f, and not *0x80 to *0xff.
+ ///
+ /// (For larger c, you _can_ solve a quadratic of this form, but
+ /// you have to do it by escalating to the next higher level. The
+ /// responsibility for this lies with the caller of this recursive
+ /// function.)
+ ///
+ /// With this condition satisfied, we're trying to find x=xh\*t+xl
+ /// such that x^2+x+c=0, i.e. such that
+ ///
+ /// ```text
+ /// 0 = x^2 + x + c
+ /// = (xh t + xl)^2 + (xh t + xl) + (ch t + cl)
+ /// = xh^2 t^2 + xl^2 + xh t + xl + ch t + cl
+ /// = xh^2 (t + h) + xl^2 + xh t + xl + ch t + cl
+ /// = (xh^2 + xh + ch) t + (xh^2 h + xl^2 + xl + cl)
+ /// ```
+ ///
+ /// Both coefficients of this must be 0. The t coefficient is a
+ /// quadratic in the same normalised form, which we start by
+ /// solving recursively (which we can do, because if c had its
+ /// high bit clear, so does ch).
+ ///
+ /// We can use that quadratic to simplify the other one,
+ /// substituting xh^2 for xh+ch to give
+ ///
+ /// ```text
+ /// xl^2 + xl + ((xh + ch) h + cl)
+ /// ```
+ ///
+ /// This is again a quadratic in normalised form, but we aren't
+ /// guaranteed that the constant term has the high bit clear. But
+ /// we have two choices for the value of xh we use, because the
+ /// first quadratic has two roots, differing by 1. Adjusting xh by
+ /// 1 in the expression ((xh+ch)h + cl) adjusts the whole thing by
+ /// h, i.e. flips the top bit. So there's always one correct
+ /// choice of xh and one wrong one, and we just have to pick the
+ /// right one.
+ fn quadratic1_recurse(self, level: usize) -> FiniteNimber {
+ match level.checked_sub(1) {
+ Some(sublevel) => {
+ let (clo, chi) = self.split(sublevel);
+ let mut xhi = chi.quadratic1_recurse(sublevel);
+ let mut c2 =
+ (xhi.to_ref() + chi).to_ref().mul_by_h(sublevel) + clo;
+ let hbit = (1 << sublevel) - 1;
+ if c2.test_bit(hbit) {
+ let one = FiniteNimber::from(1);
+ c2 += one.to_ref().mul_by_h(sublevel); // FIXME: optimise!
+ xhi += one;
+ }
+ let xlo = c2.to_ref().quadratic1_recurse(sublevel);
+ xlo.to_ref().join(xhi.to_ref(), sublevel)
+ }
+
+ // At level 0, the "high bit clear" precondition means the
+ // input nimber must be 0, so the only possible quadratic
+ // is x^2+x = x(x+1), whose solutions are 0 and 1.
+ // Arbitrarily return 0.
+ None => FiniteNimber::from(0),
+ }
+ }
+}
+
impl FiniteNimber {
/// Allows the binary representation of a `FiniteNimber` to be
/// extracted and examined by client code.
let r = self.to_ref();
r.sqrt_recurse(r.level())
}
+
+ /// Compute the solutions to a quadratic equation in the
+ /// normalised form x^2 + x + c, for some constant c. The constant
+ /// c is the receiver of the method, that is, you call
+ /// `c.quadratic1()` and get back two possible values of x.
+ ///
+ /// The two values of x are returned in numerical order.
+ ///
+ /// # Examples
+ ///
+ /// ```
+ /// use nimber::FiniteNimber;
+ ///
+ /// let c = FiniteNimber::from(0x14142135);
+ /// let (x0, x1) = c.quadratic1();
+ ///
+ /// // The two solutions are different
+ /// assert_ne!(x0, x1);
+ ///
+ /// // Both of them satisfy the original equation
+ /// assert_eq!(x0.square() + &x0 + &c, FiniteNimber::from(0));
+ /// assert_eq!(x1.square() + &x1 + &c, FiniteNimber::from(0));
+ /// ```
+ pub fn quadratic1(&self) -> (FiniteNimber, FiniteNimber) {
+ let r = self.to_ref();
+ let level = r.level();
+ let start_level = if self.test_bit((1 << level) - 1) {
+ // Increase the starting level by 1 to avoid the high bit
+ // being set
+ level + 1
+ } else {
+ level
+ };
+ let x0 = r.quadratic1_recurse(start_level);
+ let x1 = &x0 + FiniteNimber::from(1);
+ Self::sort_pair(x0, x1)
+ }
+
+ /// Compute the solutions to a quadratic equation in the monic
+ /// form x^2 + b\*x + c, for constants b,c.
+ ///
+ /// The two values of x are returned in numerical order.
+ ///
+ /// # Examples
+ ///
+ /// ```
+ /// use nimber::FiniteNimber;
+ ///
+ /// let b = FiniteNimber::from(0x16180339);
+ /// let c = FiniteNimber::from(0x14142135);
+ /// let (x0, x1) = FiniteNimber::quadratic2(&b, &c);
+ ///
+ /// // The two solutions are different
+ /// assert_ne!(x0, x1);
+ ///
+ /// // Both of them satisfy the original equation
+ /// assert_eq!(x0.square() + &b * &x0 + &c, FiniteNimber::from(0));
+ /// assert_eq!(x1.square() + &b * &x1 + &c, FiniteNimber::from(0));
+ /// ```
+ pub fn quadratic2(b: &Self, c: &Self) -> (FiniteNimber, FiniteNimber) {
+ // Strategy: if b != 0, we reduce this to the quadratic1()
+ // case by substituting x = b*y, because then the polynomial
+ //
+ // x^2 + bx + c becomes b^2 y^2 + b^2 y + c
+ //
+ // and dividing off b^2, we get y^2 + y + (c/b^2) = 0. So we
+ // solve that for y, and then recover x by multiplying by b.
+ //
+ // If b = 0, then the equation is just x^2 = c, so we take the
+ // square root of c. It only has one, so the quadratic has a
+ // repeated root.
+ if b == &FiniteNimber::from(0) {
+ let x = c.sqrt();
+ (x.clone(), x)
+ } else {
+ let (x0, x1) = (c / b.square()).quadratic1();
+ Self::sort_pair(x0 * b, x1 * b)
+ }
+ }
+
+ /// Compute the solutions to a quadratic equation in the most
+ /// general form a\*x^2 + b\*x + c, for constants a,b,c. This
+ /// function panics if `a == 0`.
+ ///
+ /// The two values of x are returned in numerical order.
+ ///
+ /// # Examples
+ ///
+ /// ```
+ /// use nimber::FiniteNimber;
+ ///
+ /// let a = FiniteNimber::from(0x17320508);
+ /// let b = FiniteNimber::from(0x16180339);
+ /// let c = FiniteNimber::from(0x14142135);
+ /// let (x0, x1) = FiniteNimber::quadratic3(&a, &b, &c);
+ ///
+ /// // The two solutions are different
+ /// assert_ne!(x0, x1);
+ ///
+ /// // Both of them satisfy the original equation
+ /// assert_eq!(&a * x0.square() + &b * &x0 + &c, FiniteNimber::from(0));
+ /// assert_eq!(&a * x1.square() + &b * &x1 + &c, FiniteNimber::from(0));
+ /// ```
+ pub fn quadratic3(
+ a: &Self,
+ b: &Self,
+ c: &Self,
+ ) -> (FiniteNimber, FiniteNimber) {
+ // Strategy: we reduce this to the quadratic2() case by
+ // dividing through by a.
+ Self::quadratic2(&(b / a), &(c / a))
+ }
}
impl_binop_wrappers!(
])
);
}
+
+ #[test]
+ fn quadratic() {
+ // Simplest possible quadratic1: the roots of x^2+x = x(x+1) are (0,1)
+ assert_eq!(
+ FiniteNimber::from(0).quadratic1(),
+ (FiniteNimber::from(0), FiniteNimber::from(1))
+ );
+
+ // Next-simplest, requiring escalation to the next higher
+ // subfield: the roots of x^2+x+1, a polynomial expressible
+ // entirely in F_2 but irreducible in that field, are (*2,*3),
+ // which are in F_4 \ F_2
+ assert_eq!(
+ FiniteNimber::from(1).quadratic1(),
+ (FiniteNimber::from(2), FiniteNimber::from(3))
+ );
+
+ // Some cases of the equation t^2 = t + h which defines the
+ // extension of each subfield to the next
+ assert_eq!(
+ FiniteNimber::from(0x8).quadratic1(),
+ (FiniteNimber::from(0x10), FiniteNimber::from(0x11))
+ );
+ assert_eq!(
+ FiniteNimber::from(0x80).quadratic1(),
+ (FiniteNimber::from(0x100), FiniteNimber::from(0x101))
+ );
+ assert_eq!(
+ FiniteNimber::from(0x80000000).quadratic1(),
+ (
+ FiniteNimber::from(0x100000000),
+ FiniteNimber::from(0x100000001)
+ )
+ );
+
+ assert_eq!(
+ FiniteNimber::from(0x23456789).quadratic1(),
+ (
+ FiniteNimber::from(0x4aec00e4),
+ FiniteNimber::from(0x4aec00e5)
+ )
+ );
+ assert_eq!(
+ FiniteNimber::from(0xabcdef01).quadratic1(),
+ (
+ FiniteNimber::from(0x15b7ac2e8),
+ FiniteNimber::from(0x15b7ac2e9)
+ )
+ );
+
+ // Degenerate case of quadratic2 with b=0. Should return a
+ // repeated root, the same as just sqrt(c).
+ let c = FiniteNimber::from(0x1234);
+ let (x0, x1) = FiniteNimber::quadratic2(&FiniteNimber::from(0), &c);
+ assert_eq!(x0, x1);
+ assert_eq!(x0, c.sqrt());
+
+ // A large random example for quadratic2
+ let b = FiniteNimber::from(&[
+ 0x6a784d9045190cfe,
+ 0x762e7160f38b4da5,
+ 0xabf7158809cf4f3c,
+ 0x2b7e151628aed2a6,
+ ]);
+ let c = FiniteNimber::from(&[
+ 0x0082efa98ec4e6c8,
+ 0x4a4093822299f31d,
+ 0x313198a2e0370734,
+ 0x3243f6a8885a308d,
+ ]);
+ // With the coefficients one way round, the solutions are the
+ // same size as the inputs
+ let (x0, x1) = FiniteNimber::quadratic2(&b, &c);
+ assert_eq!(
+ (&x0, &x1),
+ (
+ &FiniteNimber::from(&[
+ 0x7392f8dbead36e27,
+ 0x36f5dc49b334524e,
+ 0xee0d9a634c2c0d97,
+ 0x54405edcf7af4140
+ ]),
+ &FiniteNimber::from(&[
+ 0x19eab54bafca62d9,
+ 0x40dbad2940bf1feb,
+ 0x45fa8feb45e342ab,
+ 0x7f3e4bcadf0193e6
+ ])
+ )
+ );
+ assert_eq!(x0.square() + &b * &x0 + &c, FiniteNimber::from(0));
+ assert_eq!(x1.square() + &b * &x1 + &c, FiniteNimber::from(0));
+ // With the coefficients the other way round, the solutions
+ // double in size because we have to upgrade to a larger subfield
+ let (x0, x1) = FiniteNimber::quadratic2(&c, &b);
+ assert_eq!(
+ (&x0, &x1),
+ (
+ &FiniteNimber::from(&[
+ 0xfc7a27ac1e033e3b,
+ 0x1796bec3e478bf4b,
+ 0xb3befb47f55d2ca4,
+ 0x9d9791d73ee591e3,
+ 0x0082efa98ec4e6c8,
+ 0x4a4093822299f31d,
+ 0x313198a2e0370734,
+ 0x3243f6a8885a308d
+ ]),
+ &FiniteNimber::from(&[
+ 0xfcf8c80590c7d8f3,
+ 0x5dd62d41c6e14c56,
+ 0x828f63e5156a2b90,
+ 0xafd4677fb6bfa16e,
+ 0x0082efa98ec4e6c8,
+ 0x4a4093822299f31d,
+ 0x313198a2e0370734,
+ 0x3243f6a8885a308d
+ ])
+ )
+ );
+ assert_eq!(x0.square() + &c * &x0 + &b, FiniteNimber::from(0));
+ assert_eq!(x1.square() + &c * &x1 + &b, FiniteNimber::from(0));
+
+ // Finally, a large test of quadratic3
+ let a = FiniteNimber::from(&[
+ 0x1f86c6a11d0c18e9,
+ 0x41082276bf3a2725,
+ 0x5f39cc0605cedc83,
+ 0x19e3779b97f4a7c1,
+ ]);
+ let (x0, x1) = FiniteNimber::quadratic3(&a, &b, &c);
+ assert_eq!(
+ (&x0, &x1),
+ (
+ &FiniteNimber::from(&[
+ 0x9158f8e5521e4d7f,
+ 0x14039613d6356d58,
+ 0x32c757ff646c969f,
+ 0xe491100b296059ec
+ ]),
+ &FiniteNimber::from(&[
+ 0xcfe6ae8db299c447,
+ 0x5bcae8a4aab7a309,
+ 0x28fb6b4774cab164,
+ 0xfcb0c6016cb46b22
+ ])
+ )
+ );
+ assert_eq!(&a * x0.square() + &b * &x0 + &c, FiniteNimber::from(0));
+ assert_eq!(&a * x1.square() + &b * &x1 + &c, FiniteNimber::from(0));
+ }
}
~ian / nimber.git / commitdiff