From bff51ac6a06490a185eda93c9be052cd56286fcd Mon Sep 17 00:00:00 2001 From: Simon Tatham Date: Mon, 14 Apr 2025 21:01:22 +0100 Subject: [PATCH] Quadratic-solving. --- src/finitenimber.rs | 402 ++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 392 insertions(+), 10 deletions(-) diff --git a/src/finitenimber.rs b/src/finitenimber.rs index bb1b79e..19a0116 100644 --- a/src/finitenimber.rs +++ b/src/finitenimber.rs @@ -118,10 +118,10 @@ enum FiniteNimberEnum { /// 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, @@ -134,7 +134,7 @@ enum FiniteNimberEnum { /// 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 @@ -147,7 +147,7 @@ enum FiniteNimberEnum { /// [`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 @@ -159,7 +159,7 @@ enum FiniteNimberEnum { /// 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 @@ -172,7 +172,7 @@ enum FiniteNimberEnum { /// 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. @@ -488,7 +488,10 @@ impl<'a> FiniteNimberRef<'a> { } } +/// 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), @@ -498,6 +501,40 @@ impl FiniteNimber { }, } } + + /// 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> { @@ -885,8 +922,9 @@ impl<'a> 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 { match level.checked_sub(1) { Some(sublevel) => { @@ -969,6 +1007,86 @@ impl<'a> FiniteNimberRef<'a> { } } +/// 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. @@ -1034,6 +1152,118 @@ impl FiniteNimber { 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!( @@ -1486,4 +1716,156 @@ mod tests { ]) ); } + + #[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)); + } } -- 2.30.2