nimber 0.1.1

//! Internal module of the `nimber` crate dealing with finite nimbers.
//!
//! The type [`FiniteNimber`] is exported as a public type by the
//! top-level module. It's a singleton struct wrapping
//! [`FiniteNimberEnum`].

use core::fmt::{Debug, Display, Formatter};
use core::ops::{Add, Div, Mul, Neg, Sub};
use itertools::Itertools;

/// Alias for the `u64` type, used within the implementation as the
/// element type of [`FiniteNimber`] and [`FiniteNimberRef`].
type Word = u64;

/// The largest "level" of any nimber that fits in a [`Word`]. See
/// [`FiniteNimberRef`] for an explanation of levels.
const WORDLEVELS: usize = 6; // 2^{2^6} = 64 = size of Word

/// A type representing a finite nimber.
///
/// This type has no upper bound: it can represent _any_ finite
/// nimber, in principle (limited only by available memory). For this
/// reason, it implements [`Clone`] but not [`Copy`], because cloning
/// it may require a memory allocation. However, _small_ nimbers
/// (fitting in one `u64`) are stored without an extra memory
/// allocation at all, for speed.
///
/// Nimbers can be constructed using the [`From`] trait, either
/// starting from a single [`u64`] for small nimbers, or from a
/// [`Vec`] or slice of [`u64`] for larger ones. In the vector or
/// slice case, the values you provide are taken in a little-endian
/// manner, so that the first `u64` in the list represents the
/// low-order bits of the overall number.
///
/// Most operations on nimbers are performed using the ordinary `+`,
/// `*` and `/` operators (and `-` if you like, although nimbers have
/// characteristic 2, it's the same as `+`). These can be applied to
/// any combination of `FiniteNimber` and `&FiniteNimber`.
///
/// `FiniteNimber` also provides a few extra arithmetic operations as
/// ordinary methods, like [`FiniteNimber::sqrt`].
///
/// To extract the results of a nimber calculation into another type,
/// you can retrieve an `&[u64]` via the [`FiniteNimber::as_slice`]
/// method. You can also format nimbers for output using [`Display`].
///
/// # Examples
///
/// ```
/// use nimber::FiniteNimber;
///
/// // Two different ways to make the same nimber
/// assert_eq!(FiniteNimber::from(0xfedcba9876543210u64),
///            FiniteNimber::from(&[0xfedcba9876543210u64, 0]));
///
/// // Arithmetic operations via std::ops overloading
/// assert_eq!(FiniteNimber::from(0xf0f0) + FiniteNimber::from(0xff00),
///            FiniteNimber::from(0x0ff0)); // addition just looks like XOR
/// assert_eq!(FiniteNimber::from(0xf0f0) - FiniteNimber::from(0xff00),
///            FiniteNimber::from(0x0ff0)); // subtraction is same as addition
/// assert_eq!(FiniteNimber::from(0x8000) * FiniteNimber::from(0x4000),
///            FiniteNimber::from(0xb9c5)); // multiplication is weirder!
/// assert_eq!(FiniteNimber::from(0xb9c5) / FiniteNimber::from(0x8000),
///            FiniteNimber::from(0x4000)); // and division inverts it
///
/// // Extra operations via dedicated FiniteNimber methods
/// assert_eq!(FiniteNimber::from(0x8000).square(),
///            FiniteNimber::from(0xde4a)); // faster way to multiply by self
/// assert_eq!(FiniteNimber::from(0xde4a).sqrt(),
///            FiniteNimber::from(0x8000)); // square root inverts squaring
///
/// // Get the results back out as a slice of u64
/// assert_eq!(FiniteNimber::from(0x8000).square().as_slice(), &[0xde4a]);
/// let big_nimber = FiniteNimber::from(&[0, 0x8000000000000000]);
/// assert_eq!(big_nimber.square().as_slice(),
///            &[0xe3a92850a9218171, 0xde4ae3a94a88a921]);
///
/// // Or format for display, using the * prefix notation from game theory.
/// // This displays large numbers as if they were integers, so the highest-
/// // order digit comes first.
/// assert_eq!(format!("{}", FiniteNimber::from(0x1234)), "*0x1234");
/// assert_eq!(format!("{}", FiniteNimber::from(&[1, 2, 3])),
///            "*0x300000000000000020000000000000001");
/// assert_eq!(format!("{}", big_nimber.square()),
///            "*0xde4ae3a94a88a921e3a92850a9218171");
/// ```
#[derive(Clone, PartialEq, Eq, Hash)]
pub struct FiniteNimber(FiniteNimberEnum);

/// Enumeration forming the guts of [`FiniteNimber`]. Large finite
/// nimbers are stored as a `Vec<u64>`, but small ones are stored
/// directly, which means that the lowest few levels of each recursive
/// algorithm can run without needing any memory allocation.
#[derive(Clone, PartialEq, Eq, Hash)]
enum FiniteNimberEnum {
    /// 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.
    Single(Word),

    /// A nimber of level > 6 is stored in this branch, as a vector of
    /// `u64`. The vector is always non-empty, and its length is
    /// normalised so that the last element is nonzero.
    Vector(Vec<Word>),
}

/// Internal type representing a reference to a [`FiniteNimber`]. Most
/// of the actual arithmetic implementation lives in methods of this
/// type; the trait implementations on `FiniteNimber` itself are
/// mostly small wrappers.
///
/// The key concept in computing with finite nimbers is a _level_.
/// The field of finite nimbers consists of a nested sequence of
/// finite fields of size 2^(2^n): the field of order 2, nested
/// inside the field of order 4, then 16, 256, 65536, etc. A
/// nimber's "level" is the index in this sequence of the smallest
/// 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.
///
/// Most nimber calculation algorithms (except addition, which is
/// exceptionally simple) involve recursing to lower and lower levels,
/// by splitting a nimber at level n into two half-sized nimbers at
/// level n−1, then separately calculating the two half-sized parts
/// of the answer by arithmetic on the smaller nimbers, then joining
/// the two halves back together for output.
///
/// The recursive functions implementing nimber arithmetic will take a
/// 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
/// nimbers, and use that to start the recursion.
///
/// In the descriptions of these recursive algorithms, we'll refer to
/// the variable t, which is the nimber \*2^(2^(n−1)), i.e. the
/// smallest nimber that has level n. This nimber t has the useful
/// property that multiplying it by any smaller nimber works the same
/// way in the strange nimber multiplication as it does in ordinary
/// integers: that is, if a<t and b<t, then a\*t+b looks the same
/// whether you mean integer `*` and `+` or nimber `*` and `+`. So
/// [`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
/// simply.
///
/// Another important constant is `h`, which refers to the nimber
/// whose numerical value is half that of `t`. For example, when
/// splitting a 16-bit nimber into 8-bit nimbers, `t` is \*0x100, and
/// `h` is \*0x80. This value is important because nimber
/// multiplication at _any_ level has the property that t^2 = t+h. So
/// nimber calculations on numbers represented in the form `a*t+b` can
/// 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()`].
///
/// For levels above 6, a nimber looks like a list of `Word`, and a
/// `FiniteNimberRef` stores a `&[Word]`, so that `split()` simply
/// returns two references to a subrange of the same slice. For
/// smaller levels, it's necessary to do bit shifts to construct the
/// output nimber from `split`: for example, \*0x1234 must be split
/// into \*0x12 and \*0x34, and there's no `u64` in the input
/// representation that holds either of those values. Therefore
/// `FiniteNimberRef` must also be able to store an actual `u64`
/// directly.

#[derive(Clone, Copy)]
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.
    Single(Word),

    /// A nimber of level > 6 is stored in this branch, as a reference
    /// to a slice of `u64`. The slice is always non-empty, and its
    /// length is normalised so that the last element is nonzero.
    Slice(&'a [Word]),
}

impl Debug for FiniteNimber {
    /// The debug representation of a `FiniteNimber` is displayed with
    /// parentheses if it's stored internally as a single word without
    /// a subsidiary memory allocation, or square brackets if it's
    /// stored via an auxiliary `Vec<u64>`.
    ///
    /// In the latter case, the elements of the `Vec` are shown
    /// individually, in little-endian order, matching the storage
    /// representation.
    ///
    /// # Examples
    ///
    /// ```
    /// use nimber::FiniteNimber;
    ///
    /// assert_eq!(format!("{:?}", FiniteNimber::from(0x1234)),
    ///            "FiniteNimber(0x0000000000001234)");
    /// assert_eq!(format!("{:?}", FiniteNimber::from(&[0x1234, 0x5678])),
    ///            "FiniteNimber[0x0000000000001234, 0x0000000000005678]");
    /// ```
    fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), core::fmt::Error> {
        match &self.0 {
            FiniteNimberEnum::Single(v) => {
                write!(f, "FiniteNimber(0x{:016x})", v)
            }
            FiniteNimberEnum::Vector(s) => {
                write!(f, "FiniteNimber[")?;
                let mut sep: &str = "";
                for v in s {
                    write!(f, "{}0x{:016x}", sep, v)?;
                    sep = ", ";
                }
                write!(f, "]")
            }
        }
    }
}

impl Debug for FiniteNimberRef<'_> {
    fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), core::fmt::Error> {
        match self {
            FiniteNimberRef::Single(v) => {
                write!(f, "FiniteNimberRef(0x{:16x})", v)
            }
            FiniteNimberRef::Slice(s) => {
                write!(f, "FiniteNimberRef[")?;
                let mut sep: &str = "";
                for v in *s {
                    write!(f, "{}0x{:016x}", sep, v)?;
                    sep = ", ";
                }
                write!(f, "]")
            }
        }
    }
}

impl Display for FiniteNimber {
    /// Format a `FiniteNimber` for user-friendly display.
    ///
    /// The general format is a `*` indicating a nimber, followed by a
    /// hex number, e.g. `*0x1234`. However, for nimbers with value
    /// less than 10, the `0x` is omitted, so in particular you get
    /// `*0` and `*1`.
    ///
    /// At present, the only supported formatting option is the
    /// "alternate" flag. This causes the `0x` never to be omitted, so
    /// that nimbers look more consistent with each other.
    ///
    /// # Examples
    ///
    /// ```
    /// use nimber::FiniteNimber;
    ///
    /// // Some general examples
    /// assert_eq!(format!("{}", FiniteNimber::from(0x1234)), "*0x1234");
    /// let big = FiniteNimber::from(
    ///     &[0xe3a92850a9218171, 0xde4ae3a94a88a921]);
    /// assert_eq!(format!("{}", big), "*0xde4ae3a94a88a921e3a92850a9218171");
    ///
    /// // 0x is omitted for 0,1,...,9, but not for 10
    /// assert_eq!(format!("{}", FiniteNimber::from(0)), "*0");
    /// assert_eq!(format!("{}", FiniteNimber::from(1)), "*1");
    /// assert_eq!(format!("{}", FiniteNimber::from(9)), "*9");
    /// assert_eq!(format!("{}", FiniteNimber::from(10)), "*0xa");
    ///
    /// // The alternate # flag makes the 0x unconditional
    /// assert_eq!(format!("{:#}", FiniteNimber::from(0)), "*0x0");
    /// assert_eq!(format!("{:#}", FiniteNimber::from(1)), "*0x1");
    /// assert_eq!(format!("{:#}", FiniteNimber::from(9)), "*0x9");
    /// ```
    fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), core::fmt::Error> {
        let r = self.to_ref();
        if !f.alternate() && r.level() < WORDLEVELS && r.low_word() < 10 {
            write!(f, "*{}", r.low_word())
        } else {
            write!(f, "*0x")?;
            let mut started = false;
            for v in r.as_slice().iter().rev() {
                if started {
                    write!(f, "{:016x}", v)?;
                } else if *v != 0 {
                    write!(f, "{:x}", v)?;
                    started = true;
                }
            }
            if !started {
                write!(f, "0")?;
            }
            Ok(())
        }
    }
}

impl Default for FiniteNimber {
    /// The default value for a nimber is zero.
    fn default() -> Self {
        Self(FiniteNimberEnum::Single(0))
    }
}

/// Some low-level administrative methods.
impl<'a> FiniteNimberRef<'a> {
    /// The `FiniteNimberRef` representing zero.
    fn zero() -> Self {
        FiniteNimberRef::Single(0)
    }

    /// Return the low-order `Word` of a `FiniteNimberRef`. If you've
    /// determined via [`FiniteNimberRef::level()`] that its level is
    /// at most [`WORDLEVELS`], then this is the _whole_ of the value.
    fn low_word(self) -> Word {
        match self {
            FiniteNimberRef::Single(v) => v,
            FiniteNimberRef::Slice(s) => {
                *(s.first().expect("FiniteNimberRef::Slice is never empty"))
            }
        }
    }

    /// Return a `&[Word]` representing the contents of a
    /// `FiniteNimberRef`. Because a single word might be stored
    /// directly in the `FiniteNimberRef` itself, this method requires
    /// a `&FiniteNimberRef`: it can't work on a `FiniteNimberRef`
    /// passed by value.
    fn as_slice<'b, 'c>(&'b self) -> &'c [Word]
    where
        'a: 'c,
        'b: 'c,
    {
        match self {
            FiniteNimberRef::Single(v) => core::slice::from_ref(v),
            FiniteNimberRef::Slice(s) => s,
        }
    }

    /// Calculate the "level" of a finite nimber.
    fn level(self) -> usize {
        fn word_level(w: Word) -> usize {
            let log_w: u32 =
                (Word::BITS - 1).saturating_sub(w.leading_zeros());
            u32::BITS.saturating_sub(log_w.leading_zeros()) as usize
        }

        fn usize_level(w: usize) -> usize {
            let log_w: u32 = (usize::BITS - 1) - w.leading_zeros();
            log_w as usize
        }

        match self {
            FiniteNimberRef::Single(v) => word_level(v),
            FiniteNimberRef::Slice(s) => {
                let w = s
                    .len()
                    .checked_sub(1)
                    .expect("slice representation must always be non-empty");
                match w {
                    0 => word_level(s[0]),
                    w => usize_level(w) + (WORDLEVELS + 1),
                }
            }
        }
    }

    /// Split a nimber into two half-sized nimbers, with the output
    /// nimbers having level at most `level`. If the output
    /// `FiniteNimberRef` are large enough, this is implemented by
    /// returning sub-slice references into the slice held by the
    /// original one, so this operation is fast.
    ///
    /// The returned tuple from `split()` contains the low half and
    /// then the high half.
    ///
    /// For example, if `x` is the nimber \*0x1234 and you split it at
    /// level 3: `t` is the nimber *(2^(2^3)), or \*0x100, so you can
    /// write `x = xhi * t + xlo` with `xhi` being \*0x12 and `xlo`
    /// being \*0x34. So `x.split(3)` gives (\*0x34,\*0x12) (or rather,
    /// `FiniteNimberRef` objects containing those values).
    ///
    /// If the input nimber has level greater than `level+1`, then
    /// this function will return its low 2^(2^(level+1)) bits,
    /// discarding any higher-order bits. For example, with `x` equal
    /// to \*0x1234 as above, `x.split(2)` will give (\*0x4, \*0x3),
    /// discarding the 0x12 at the top. This probably isn't useful.
    fn split(
        self,
        level: usize,
    ) -> (FiniteNimberRef<'a>, FiniteNimberRef<'a>) {
        match level.checked_sub(WORDLEVELS) {
            None => {
                // If level < WORDLEVELS then we can cast to Word any old way
                let bits = 1 << (level as Word);
                let mask = (1 << bits) - 1;
                let v = self.low_word();
                return (
                    FiniteNimberRef::Single(v & mask),
                    FiniteNimberRef::Single((v >> bits) & mask),
                );
            }
            Some(wordlevel) => match self {
                FiniteNimberRef::Single(v) => {
                    (FiniteNimberRef::Single(v), Self::zero())
                }
                FiniteNimberRef::Slice(s) => {
                    match wordlevel
                        .try_into()
                        .ok()
                        .and_then(|w| 1usize.checked_shl(w))
                    {
                        Some(words) => {
                            let mut iter = s.chunks(words);
                            let lo = iter.next().expect(
                                "FiniteNimberRef::Slice is never empty",
                            );
                            match iter.next() {
                                Some(hi) => (
                                    FiniteNimberRef::Slice(lo),
                                    FiniteNimberRef::Slice(hi),
                                ),
                                None => {
                                    (FiniteNimberRef::Slice(lo), Self::zero())
                                }
                            }
                        }
                        None => (self.clone(), Self::zero()),
                    }
                }
            },
        }
    }

    /// Join two half-sized nimbers into a bigger one, with the input
    /// nimbers having level at most `level`. This can't be done in
    /// general in a way that generates a `FiniteNimberRef` as output:
    /// it requires allocating a new `FiniteNimber`.
    ///
    /// This is the inverse of [`FiniteNimberRef::split()`], and takes
    /// the same `level` value (hence why `level` is the size of the
    /// outputs for `split` and the inputs for `join`).
    ///
    /// For example, if x and y are the nimbers \*0x34 and \*0x12
    /// respectively, then `x.join(y, 3)` is the nimber \*0x1234, and
    /// `x.join(y, 4)` is \*0x120034.
    ///
    /// As with `split`, if the input nimbers are too large, then
    /// their higher bits are discarded. For example, with `x`,`y` as
    /// above, `x.join(y, 2)` gives \*0x24, discarding the high digit
    /// of each input.
    fn join(self, hi: FiniteNimberRef<'_>, level: usize) -> FiniteNimber {
        match level.checked_sub(WORDLEVELS) {
            None => {
                // If level < WORDLEVELS then we can cast to Word any old way
                let bits = 1 << (level as Word);
                let mask = (1 << bits) - 1;
                let vlo = self.low_word();
                let vhi = hi.low_word();
                ((vlo & mask) | ((vhi & mask) << bits)).into()
            }
            Some(wordlevel) => {
                match wordlevel
                    .try_into()
                    .ok()
                    .and_then(|w| 1usize.checked_shl(w))
                {
                    Some(words) => {
                        let slo = self.as_slice();
                        let shi = hi.as_slice();
                        let padding = words.saturating_sub(slo.len());
                        let vec: Vec<_> = slo
                            .iter()
                            .cloned()
                            .take(words)
                            .chain(core::iter::repeat(0).take(padding))
                            .chain(shi.iter().cloned().take(words))
                            .collect();
                        FiniteNimber::from(vec)
                    }
                    None => FiniteNimber::from(self.as_slice()),
                }
            }
        }
    }
}

/// 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),
            FiniteNimberEnum::Vector(v) => match v.len() {
                0 => FiniteNimberRef::zero(),
                _ => FiniteNimberRef::Slice(&v),
            },
        }
    }

    /// 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> {
    fn from(n: &'a FiniteNimber) -> Self {
        n.to_ref()
    }
}

impl From<FiniteNimberRef<'_>> for FiniteNimber {
    fn from(n: FiniteNimberRef) -> Self {
        match n {
            FiniteNimberRef::Single(w) => Self(FiniteNimberEnum::Single(w)),
            FiniteNimberRef::Slice(v) => FiniteNimber::from(v),
        }
    }
}

impl From<Word> for FiniteNimber {
    /// Make a `FiniteNimber` out of a single integer.
    fn from(val: Word) -> FiniteNimber {
        Self(FiniteNimberEnum::Single(val))
    }
}

impl From<Vec<Word>> for FiniteNimber {
    /// Make a `FiniteNimber` out of a list of integers provided as a
    /// `Vec<u64>`. These 64-bit integers are regarded as representing
    /// one big integer, in little-endian order, so that the first
    /// integer in the list is treated as the least significant 64
    /// bits of the represented integer.
    fn from(mut vec: Vec<Word>) -> FiniteNimber {
        match vec
            .iter()
            .enumerate()
            .rev()
            .find(|(_index, word)| **word != 0)
        {
            None => Self(FiniteNimberEnum::Single(0)),
            Some((0, w)) => Self(FiniteNimberEnum::Single(*w)),
            Some((pos, _)) => {
                vec.truncate(pos + 1);
                Self(FiniteNimberEnum::Vector(vec))
            }
        }
    }
}

impl From<&[Word]> for FiniteNimber {
    /// Make a `FiniteNimber` out of a list of integers provided as a
    /// slice reference. These 64-bit integers are regarded as
    /// representing one big integer, in little-endian order, so that
    /// the first integer in the list is treated as the least
    /// significant 64 bits of the represented integer.
    fn from(slice: &[Word]) -> FiniteNimber {
        match slice
            .iter()
            .enumerate()
            .rev()
            .find(|(_index, word)| **word != 0)
        {
            None => Self(FiniteNimberEnum::Single(0)),
            Some((0, w)) => Self(FiniteNimberEnum::Single(*w)),
            Some((pos, _)) => {
                Self(FiniteNimberEnum::Vector(slice[..=pos].into()))
            }
        }
    }
}

impl<const N: usize> From<&[Word; N]> for FiniteNimber {
    /// Make a `FiniteNimber` out of a list of integers provided as an
    /// array reference. These 64-bit integers are regarded as
    /// representing one big integer, in little-endian order, so that
    /// the first integer in the list is treated as the least
    /// significant 64 bits of the represented integer.
    fn from(array: &[Word; N]) -> FiniteNimber {
        FiniteNimber::from(array as &[Word])
    }
}

/// Macro to implement a binary operation and its associated
/// assignment operator on all meaningful combinations of
/// `FiniteNimber` and `&FiniteNimber`.
macro_rules! impl_binop_wrappers {
    ($trait:tt, $fn:ident, $assigntrait:tt, $assignfn:ident, $comment:literal) => {
        impl $trait<FiniteNimberRef<'_>> for FiniteNimber {
            type Output = FiniteNimber;
            fn $fn(self, other: FiniteNimberRef<'_>) -> FiniteNimber {
                let aref: FiniteNimberRef = (&self).into();
                let bref: FiniteNimberRef = other;
                $trait::$fn(aref, bref)
            }
        }

        impl $trait<FiniteNimber> for FiniteNimberRef<'_> {
            type Output = FiniteNimber;
            fn $fn(self, other: FiniteNimber) -> FiniteNimber {
                let aref: FiniteNimberRef = self;
                let bref: FiniteNimberRef = (&other).into();
                $trait::$fn(aref, bref)
            }
        }

        impl $trait<&FiniteNimber> for &FiniteNimber {
            /// The result of combining two `FiniteNimber` with an
            /// arithmetic operation is another `FiniteNimber`, which
            /// must be newly allocated.
            type Output = FiniteNimber;
            #[doc=$comment]
            fn $fn(self, other: &FiniteNimber) -> FiniteNimber {
                let aref: FiniteNimberRef = self.into();
                let bref: FiniteNimberRef = other.into();
                $trait::$fn(aref, bref)
            }
        }

        impl $trait<&FiniteNimber> for FiniteNimber {
            /// The result of combining two `FiniteNimber` with an
            /// arithmetic operation is another `FiniteNimber`, which
            /// must be newly allocated.
            type Output = FiniteNimber;
            #[doc=$comment]
            fn $fn(self, other: &FiniteNimber) -> FiniteNimber {
                let aref: FiniteNimberRef = (&self).into();
                let bref: FiniteNimberRef = other.into();
                $trait::$fn(aref, bref)
            }
        }

        impl $trait<FiniteNimber> for &FiniteNimber {
            /// The result of combining two `FiniteNimber` with an
            /// arithmetic operation is another `FiniteNimber`, which
            /// must be newly allocated.
            type Output = FiniteNimber;
            #[doc=$comment]
            fn $fn(self, other: FiniteNimber) -> FiniteNimber {
                let aref: FiniteNimberRef = self.into();
                let bref: FiniteNimberRef = (&other).into();
                $trait::$fn(aref, bref)
            }
        }

        impl $trait<FiniteNimber> for FiniteNimber {
            /// The result of combining two `FiniteNimber` with an
            /// arithmetic operation is another `FiniteNimber`, which
            /// must be newly allocated.
            type Output = FiniteNimber;
            #[doc=$comment]
            fn $fn(self, other: FiniteNimber) -> FiniteNimber {
                let aref: FiniteNimberRef = (&self).into();
                let bref: FiniteNimberRef = (&other).into();
                $trait::$fn(aref, bref)
            }
        }

        impl core::ops::$assigntrait<&FiniteNimber> for FiniteNimber {
            #[doc=$comment]
            fn $assignfn(&mut self, other: &FiniteNimber) {
                let aref: FiniteNimberRef = (&self as &FiniteNimber).into();
                let bref: FiniteNimberRef = other.into();
                *self = $trait::$fn(aref, bref);
            }
        }

        impl core::ops::$assigntrait<FiniteNimber> for FiniteNimber {
            #[doc=$comment]
            fn $assignfn(&mut self, other: FiniteNimber) {
                let aref: FiniteNimberRef = (&self as &FiniteNimber).into();
                let bref: FiniteNimberRef = (&other).into();
                *self = $trait::$fn(aref, bref);
            }
        }
    };
}

impl<'a, 'b> Add<FiniteNimberRef<'a>> for FiniteNimberRef<'b> {
    type Output = FiniteNimber;
    fn add(self, other: FiniteNimberRef<'a>) -> FiniteNimber {
        let vec: Vec<_> = self
            .as_slice()
            .iter()
            .cloned()
            .zip_longest(other.as_slice().iter().cloned())
            .map(|pair| pair.reduce(|a, b| a ^ b))
            .collect();
        FiniteNimber::from(vec)
    }
}

impl<'a, 'b> Sub<FiniteNimberRef<'a>> for FiniteNimberRef<'b> {
    type Output = FiniteNimber;
    fn sub(self, other: FiniteNimberRef<'a>) -> FiniteNimber {
        self + other
    }
}

/// Internal functions having to do with multiplication.
impl<'a> FiniteNimberRef<'a> {
    /// Recursively multiply a nimber by the special value
    /// \*(2^(2^level−1)).
    ///
    /// This is a separate method from multiplication for two reasons.
    /// Firstly, it's more efficient than general multiplication
    /// (because if h itself is split into two halves, one of the
    /// halves is zero and can be ignored); secondly, it's a necessary
    /// _subroutine_ of general multiplication.
    ///
    /// When we call this routine from another context, we think of it
    /// as multiplying by h, where h is the value in the squaring
    /// equation t^2 = t+h described above. But _within_ the context
    /// of this method, `level` is one less, so we have our own
    /// smaller value of h, and the value our caller wants us to
    /// multiply by is h\*t, in our terms.
    ///
    /// The formula for multiplication by h\*t is
    ///
    /// ```text
    ///   (ah t + al) h t
    /// = ah h t^2 + al h
    /// = ah h (t + h) + al h
    /// = (ah h) t + (ah h + al) h
    /// ```
    fn mul_by_h(self, level: usize) -> FiniteNimber {
        match level.checked_sub(1) {
            Some(sublevel) => {
                let (lo, hi) = self.split(sublevel);
                hi.mul_by_h(sublevel)
                    .to_ref()
                    .mul_by_h(sublevel)
                    .to_ref()
                    .join(
                        (lo + hi).to_ref().mul_by_h(sublevel).to_ref(),
                        sublevel,
                    )
            }

            // At level 0, h = 1, so multiplying by h is a no-op.
            None => self.into(),
        }
    }

    /// Recursively calculate the square of a nimber. Underlies the
    /// public method [`FiniteNimber::square`], and is also used as a
    /// subroutine in other arithmetic implementations.
    ///
    /// This is faster than general multiplication because squaring a
    /// sum in any field of characteristic 2 has a particularly simple
    /// formula: instead of (a+b)^2 = a^2+2ab+b^2 as usual, the
    /// cross-term 2ab vanishes, because we're working mod 2. So you
    /// just get (a+b)^2 = a^2 + b^2, which looks like a schoolchild
    /// mistake, but in this context, is perfectly true!
    ///
    /// So the formula for squaring a = ah\*t+al is
    ///
    /// ```text
    ///   (ah t + al)^2
    /// = ah^2 t^2 + al^2
    /// = ah^2 (t + h) + al^2
    /// = (ah^2) t + (ah^2 h + al^2)
    /// ```
    fn square_recurse(self, level: usize) -> FiniteNimber {
        match level.checked_sub(1) {
            Some(sublevel) => {
                let (lo, hi) = self.split(sublevel);
                let lo2 = lo.square_recurse(sublevel);
                let hi2 = hi.square_recurse(sublevel);
                (lo2 + hi2.to_ref().mul_by_h(sublevel))
                    .to_ref()
                    .join(hi2.to_ref(), sublevel)
            }

            // At level 0, both elements of GF(2) square to themselves.
            None => self.into(),
        }
    }
}

impl<'a, 'b> Mul<FiniteNimberRef<'a>> for FiniteNimberRef<'b> {
    type Output = FiniteNimber;

    /// Recursively multiply two nimbers. Underlies the [`Mul`]
    /// implementations on `FiniteNimber` itself.
    ///
    /// If the two input nimbers have the same level, that's the
    /// interesting case. The most obvious formula for multiplying a =
    /// ah\*t+al by b = bh\*t+bl is
    ///
    /// ```text
    ///   (ah t + al)(bh t + bl)
    /// = ah bh t^2 + ah bl t + al bh t + al bl
    /// = ah bh (t + h) + ah bl t + al bh t + al bl
    /// = (ah bh + ah bl + al bh) t + (ah bh h + al bl)
    /// ```
    ///
    /// But this requires four smaller multiplications, which can be
    /// reduced to three by a variant of the Karatsuba optimisation.
    /// First compute the following value, using a single half-sized
    /// multiplication:
    ///
    /// ```text
    /// k = (ah + al)(bh + bl) = ah bh + ah bl + al bh + al bl
    /// ```
    ///
    /// Then the `t` coefficient in the above product formula contains
    /// three of the four terms in `k`. The missing one is `al bl`,
    /// which we needed to compute anyway for the other coefficient.
    ///
    /// So after computing `k`, we use two more half-sized
    /// multiplications to compute `ah bh` and `al bl`, and then
    /// combine those three smaller products using the following
    /// formula, which uses `al bl` twice, and involves nothing but
    /// addition and multiplication by the special value h:
    ///
    /// ```text
    /// (k + al bl) t + (ah bh h + al bl)
    /// ```
    ///
    /// However, if we're given two nimbers to multiply _not_ at the
    /// same level, we can do something much easier: only split up the
    /// larger nimber (say b), and simply multiply each of its halves
    /// by the smaller one in the obvious way:
    ///
    /// ```text
    ///   a(bh t + bl)
    /// = (a bh) t + (a bl)
    /// ```
    fn mul(self, other: FiniteNimberRef<'a>) -> FiniteNimber {
        let slevel = self.level();
        let olevel = other.level();

        // Sort the two inputs so that alevel <= blevel
        let (a, b, alevel, blevel) = if slevel > olevel {
            (other, self, olevel, slevel)
        } else {
            (self, other, slevel, olevel)
        };

        if alevel < blevel {
            // The easy case, where we only need to split up b and
            // recurse twice to the next level down.
            let sublevel = blevel - 1;
            let (blo, bhi) = b.split(sublevel);
            (blo * a).to_ref().join((bhi * a).to_ref(), sublevel)
        } else {
            // The hard case, where we must split up both nimbers and
            // do three full recursive multiplications plus a mul_by_h.
            match alevel.checked_sub(1) {
                Some(sublevel) => {
                    let (alo, ahi) = a.split(sublevel);
                    let (blo, bhi) = b.split(sublevel);
                    let karatsuba = (alo + ahi) * (blo + bhi);
                    let albl = alo * blo;
                    let ahbh = ahi * bhi;
                    (&albl + ahbh.to_ref().mul_by_h(sublevel))
                        .to_ref()
                        .join((karatsuba + &albl).to_ref(), sublevel)
                }

                // At level 0, we're in GF(2), so multiplication looks
                // like bitwise AND.
                None => FiniteNimber::from(a.low_word() & b.low_word()),
            }
        }
    }
}

/// Internal functions having to do with division.
impl<'a> FiniteNimberRef<'a> {
    /// Recursively calculate the multiplicative inverse of a finite
    /// nimber. Underlies the [`Div`] implementations on
    /// `FiniteNimber` itself, which compute a/b by finding the
    /// inverse of b and then multiplying by a.
    ///
    /// We do this by rearranging into a form where we only have to
    /// divide by a nimber at the next smaller level, and recursing.
    ///
    /// There are several ways to derive the same formula for the
    /// inverse of a number; one involves inverting a matrix, and one
    /// involves highfalutin' stuff about Galois conjugates (an
    /// extension of the standard trick for dividing by a complex
    /// number by multiplying top and bottom by its complex
    /// conjugate). I think this derivation is simpler than either:
    ///
    /// We're after a number u such that (uh t + ul)(ah t + al),
    /// multiplied out and reduced via t^2=t+h so that it's a
    /// polynomial in t of degree < 2, comes out to 0t+1. So we need
    /// its t coefficient to be 0. By the unoptimised version of the
    /// product formula, that t coefficient comes to
    ///
    /// ```text
    /// (uh ah + uh al + ul ah) = (ah + al) uh + (ah) ul
    /// ```
    ///
    /// To make that zero, we need to choose any two numbers uh,ul in
    /// the right ratio, that is, such that uh/ul = (ah)/(ah+al). One
    /// very obvious choice is to take uh=ah and ul=ah+al, so that
    /// those two fractions are obviously equal because we're dividing
    /// the exact same pair of numbers on both sides.
    ///
    /// So the product (ah t + ah + al)(ah t + al) has t coefficient
    /// zero. So we can work out what the rest of it comes to, and
    /// divide by that to scale the output to be 1.
    ///
    /// Omitting the boring algebra of multiplying that out, this
    /// leads to the following formula for the inverse of (ah t + al):
    ///
    /// ```text
    /// (ah t + ah + al) / (al (ah + al) + ah^2 h)
    /// ```
    ///
    /// 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) => {
                let (lo, hi) = self.split(sublevel);
                let sum = lo + hi;
                let sq = hi.square_recurse(sublevel);
                let det = lo * sum.to_ref() + sq.to_ref().mul_by_h(sublevel);
                let detinv = det.to_ref().inverse_recurse(sublevel)?;
                Some(
                    (detinv.to_ref() * sum)
                        .to_ref()
                        .join((detinv.to_ref() * hi).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.
            None => match self.low_word() {
                0 => None,
                v => Some(FiniteNimber::from(v)),
            },
        }
    }

    /// Compute the multiplicative inverse of a nimber. Wrapper around
    /// `inverse_recurse` which also chooses the starting level.
    fn inverse(&self) -> Option<FiniteNimber> {
        self.inverse_recurse(self.level())
    }
}

impl<'a, 'b> Div<FiniteNimberRef<'a>> for FiniteNimberRef<'b> {
    type Output = FiniteNimber;
    fn div(self, other: FiniteNimberRef<'a>) -> FiniteNimber {
        let inverse = other.inverse().expect("Division by zero");
        self * inverse.to_ref()
    }
}

/// Internal functions having to do with square root.
impl<'a> FiniteNimberRef<'a> {
    /// Recursively calculate the square root of a nimber. Underlies
    /// the public method [`FiniteNimber::sqrt`].
    ///
    /// To derive a formula for square root, we invert the formula for
    /// squaring. We want r such that r^2 = a. The documentation for
    /// [`FiniteNimberRef::square_recurse`] derives the squaring
    /// formula as
    ///
    /// ```text
    /// (rh t + rl)^2 = (rh^2) t + (rh^2 h + rl^2)
    /// ```
    ///
    /// Setting this equal to the input `ah t + al` and equating
    /// coefficients, we get two equations
    ///
    /// ```text
    /// rh^2          = ah
    /// rh^2 h + rl^2 = al
    /// ```
    ///
    /// Simplifying the second equation using the first and
    /// rearranging gives
    ///
    /// ```text
    /// rh^2 = ah
    /// rl^2 = ah h + al
    /// ```
    ///
    /// and so we can compute each half of `r` by calculating the RHS
    /// of one of those equations and taking its square root.
    fn sqrt_recurse(self, level: usize) -> FiniteNimber {
        match level.checked_sub(1) {
            Some(sublevel) => {
                let (lo, hi) = self.split(sublevel);
                let hi_root = hi.sqrt_recurse(sublevel);
                let sum = hi.mul_by_h(sublevel) + lo;
                let sum_root = sum.to_ref().sqrt_recurse(sublevel);
                sum_root.to_ref().join(hi_root.to_ref(), sublevel)
            }

            // At level 0, both elements of GF(2) square-root to themselves.
            None => self.into(),
        }
    }
}

/// 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.
    ///
    /// The returned slice of `u64` is arranged in little-endian
    /// order, exactly the same as you would pass to `FiniteNimber::from`.
    ///
    /// The slice always has at least one element. The last element in
    /// the slice (representing the highest-order 64-bit chunk) is
    /// always nonzero, except when the nimber is 0, in which case the
    /// slice has a single element 0.
    ///
    /// # Examples
    ///
    /// ```
    /// use nimber::FiniteNimber;
    ///
    /// assert_eq!(FiniteNimber::default().as_slice(), &[0]);
    /// assert_eq!(FiniteNimber::from(1).as_slice(), &[1]);
    /// assert_eq!(FiniteNimber::from(&[1, 2, 0, 0, 0]).as_slice(), &[1, 2]);
    /// ```
    pub fn as_slice(&self) -> &[Word] {
        match &self.0 {
            FiniteNimberEnum::Single(w) => core::slice::from_ref(w),
            FiniteNimberEnum::Vector(v) => &v,
        }
    }

    /// Compute the multiplicative inverse of a nimber, i.e. the same
    /// as `1/x`. Returns an `Option` so that it can return `None` in
    /// the case of division by zero.
    ///
    /// Nimber multiplication is faster than division, so if you're
    /// dividing more than one nimber by the same value, it's faster
    /// to compute its inverse once and then multiply all the input
    /// values by that.
    pub fn inverse(&self) -> Option<FiniteNimber> {
        self.to_ref().inverse()
    }

    /// Compute the square of a nimber, just like multiplying it by
    /// itself in the ordinary way. This method is faster than the
    /// general multiplication algorithm performed by the `Mul` trait,
    /// so if you _know_ you're squaring a number, it helps to use
    /// this method instead of writing `&x * &x`.
    ///
    /// # Examples
    ///
    /// ```
    /// use nimber::FiniteNimber;
    ///
    /// let x = FiniteNimber::from(&[0x3141592653589793, 0x2718281828459045]);
    /// assert_eq!(x.square(), &x * &x);
    /// ```
    pub fn square(&self) -> FiniteNimber {
        let r = self.to_ref();
        r.square_recurse(r.level())
    }

    /// Compute the square root of a nimber. Every nimber has a unique
    /// square root, and the square root of a finite nimber is finite.
    /// So this function can't fail, and doesn't need to return a list
    /// or a tuple or anything more complicated than a single nimber.
    ///
    /// # Examples
    ///
    /// ```
    /// use nimber::FiniteNimber;
    ///
    /// let x = FiniteNimber::from(&[0x3141592653589793, 0x2718281828459045]);
    /// assert_eq!(x.square().sqrt(), x);
    /// assert_eq!(x.sqrt().square(), x);
    /// ```
    pub fn sqrt(&self) -> 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. Returns
    /// `None` 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).unwrap();
    ///
    /// // 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,
    ) -> Option<(FiniteNimber, FiniteNimber)> {
        // Strategy: we reduce this to the quadratic2() case by
        // dividing through by a.
        let ainv = a.inverse()?;
        Some(Self::quadratic2(&(b * &ainv), &(c * &ainv)))
    }
}

impl_binop_wrappers!(
    Add,
    add,
    AddAssign,
    add_assign,
    r##"

Addition of finite nimbers behaves like bitwise XOR on unsigned
integers.

"##
);
impl_binop_wrappers!(
    Sub,
    sub,
    SubAssign,
    sub_assign,
    r##"

Subtraction of finite nimbers behaves just like addition, i.e. like
bitwise XOR on unsigned integers. This is due to the field of nimbers
having characteristic 2, i.e. it is such that `x + x = 0` for any `x`,
so that `-x` is the same as `x`.

"##
);
impl_binop_wrappers!(
    Mul,
    mul,
    MulAssign,
    mul_assign,
    r##"

Multiplication of finite nimbers has the property that if the two
input nimbers both fitted in 2^n bits, then so does the output.

"##
);
impl_binop_wrappers!(
    Div,
    div,
    DivAssign,
    div_assign,
    r##"

Division of finite nimbers has the property that if the two input
nimbers both fitted in 2^n bits, then so does the output.

Dividing by the zero nimber causes a panic.

"##
);

impl Neg for FiniteNimber {
    /// The result of negating a `FiniteNimber` is another
    /// `FiniteNimber`.
    type Output = Self;

    /// Since nimbers have characteristic 2, negating a nimber leaves
    /// it unchanged. So the unary - operator has no effect.
    fn neg(self) -> Self {
        self
    }
}

impl<'a> Neg for &'a FiniteNimber {
    /// Since negation of nimbers is the identity function, this is
    /// the one arithmetic operation where we can return a reference,
    /// namely the same reference passed in.
    type Output = &'a FiniteNimber;

    /// Since nimbers have characteristic 2, negating a nimber leaves
    /// it unchanged. So the unary - operator has no effect.
    fn neg(self) -> Self {
        self
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn split() {
        let a = FiniteNimber::from(&[
            0x8786858483828180,
            0x9796959493929190,
            0xa7a6a5a4a3a2a1a0,
            0xb7b6b5b4b3b2b1b0,
            0xc7c6c5c4c3c2c1c0,
            0xd7d6d5d4d3d2d1d0,
        ]);

        let (loref, hiref) = a.to_ref().split(4);
        let lo: FiniteNimber = loref.into();
        let hi: FiniteNimber = hiref.into();
        assert_eq!(lo, FiniteNimber::from(&[0x8180]));
        assert_eq!(hi, FiniteNimber::from(&[0x8382]));

        let (loref, hiref) = a.to_ref().split(5);
        let lo: FiniteNimber = loref.into();
        let hi: FiniteNimber = hiref.into();
        assert_eq!(lo, FiniteNimber::from(&[0x83828180]));
        assert_eq!(hi, FiniteNimber::from(&[0x87868584]));

        let (loref, hiref) = a.to_ref().split(6);
        let lo: FiniteNimber = loref.into();
        let hi: FiniteNimber = hiref.into();
        assert_eq!(lo, FiniteNimber::from(&[0x8786858483828180]));
        assert_eq!(hi, FiniteNimber::from(&[0x9796959493929190]));

        let (loref, hiref) = a.to_ref().split(7);
        let lo: FiniteNimber = loref.into();
        let hi: FiniteNimber = hiref.into();
        assert_eq!(
            lo,
            FiniteNimber::from(&[0x8786858483828180, 0x9796959493929190])
        );
        assert_eq!(
            hi,
            FiniteNimber::from(&[0xa7a6a5a4a3a2a1a0, 0xb7b6b5b4b3b2b1b0,])
        );

        let (loref, hiref) = a.to_ref().split(8);
        let lo: FiniteNimber = loref.into();
        let hi: FiniteNimber = hiref.into();
        assert_eq!(
            lo,
            FiniteNimber::from(&[
                0x8786858483828180,
                0x9796959493929190,
                0xa7a6a5a4a3a2a1a0,
                0xb7b6b5b4b3b2b1b0
            ])
        );
        assert_eq!(
            hi,
            FiniteNimber::from(&[0xc7c6c5c4c3c2c1c0, 0xd7d6d5d4d3d2d1d0])
        );

        let (loref, hiref) = a.to_ref().split(9);
        let lo: FiniteNimber = loref.into();
        let hi: FiniteNimber = hiref.into();
        assert_eq!(
            lo,
            FiniteNimber::from(&[
                0x8786858483828180,
                0x9796959493929190,
                0xa7a6a5a4a3a2a1a0,
                0xb7b6b5b4b3b2b1b0,
                0xc7c6c5c4c3c2c1c0,
                0xd7d6d5d4d3d2d1d0
            ])
        );
        assert_eq!(hi, FiniteNimber::from(&[0]));

        let a = FiniteNimber::from(&[0xFF]);

        let (loref, hiref) = a.to_ref().split(0);
        let lo: FiniteNimber = loref.into();
        let hi: FiniteNimber = hiref.into();
        assert_eq!(lo, FiniteNimber::from(&[1]));
        assert_eq!(hi, FiniteNimber::from(&[1]));

        let (loref, hiref) = a.to_ref().split(1);
        let lo: FiniteNimber = loref.into();
        let hi: FiniteNimber = hiref.into();
        assert_eq!(lo, FiniteNimber::from(&[3]));
        assert_eq!(hi, FiniteNimber::from(&[3]));

        let (loref, hiref) = a.to_ref().split(2);
        let lo: FiniteNimber = loref.into();
        let hi: FiniteNimber = hiref.into();
        assert_eq!(lo, FiniteNimber::from(&[0xf]));
        assert_eq!(hi, FiniteNimber::from(&[0xf]));
    }

    #[test]
    fn join() {
        let a = FiniteNimber::from(0x5);
        let b = FiniteNimber::from(0xa);
        assert_eq!(a.to_ref().join(b.to_ref(), 0), FiniteNimber::from(0b01));
        assert_eq!(b.to_ref().join(a.to_ref(), 0), FiniteNimber::from(0b10));
        assert_eq!(a.to_ref().join(b.to_ref(), 1), FiniteNimber::from(0b1001));
        assert_eq!(b.to_ref().join(a.to_ref(), 1), FiniteNimber::from(0b0110));
        assert_eq!(a.to_ref().join(b.to_ref(), 2), FiniteNimber::from(0xa5));
        assert_eq!(b.to_ref().join(a.to_ref(), 2), FiniteNimber::from(0x5a));
        assert_eq!(
            a.to_ref().join(b.to_ref(), 5),
            FiniteNimber::from(0xa00000005)
        );
        assert_eq!(
            a.to_ref().join(b.to_ref(), 6),
            FiniteNimber::from(&[0x5, 0xa])
        );
        assert_eq!(
            a.to_ref().join(b.to_ref(), 8),
            FiniteNimber::from(&[0x5, 0, 0, 0, 0xa])
        );
    }

    #[test]
    fn levels() {
        assert_eq!(FiniteNimber::from(&[]).to_ref().level(), 0);
        assert_eq!(FiniteNimber::from(&[0]).to_ref().level(), 0);
        assert_eq!(FiniteNimber::from(&[1]).to_ref().level(), 0);
        assert_eq!(FiniteNimber::from(&[2]).to_ref().level(), 1);
        assert_eq!(FiniteNimber::from(&[3]).to_ref().level(), 1);
        assert_eq!(FiniteNimber::from(&[4]).to_ref().level(), 2);
        assert_eq!(FiniteNimber::from(&[0xf]).to_ref().level(), 2);
        assert_eq!(FiniteNimber::from(&[0x10]).to_ref().level(), 3);
        assert_eq!(FiniteNimber::from(&[0xff]).to_ref().level(), 3);
        assert_eq!(FiniteNimber::from(&[0x0100]).to_ref().level(), 4);
        assert_eq!(FiniteNimber::from(&[0xffff]).to_ref().level(), 4);
        assert_eq!(FiniteNimber::from(&[0x00010000]).to_ref().level(), 5);
        assert_eq!(FiniteNimber::from(&[0xffffffff]).to_ref().level(), 5);
        assert_eq!(
            FiniteNimber::from(&[0x0000000100000000]).to_ref().level(),
            6
        );
        assert_eq!(
            FiniteNimber::from(&[0xffffffffffffffff]).to_ref().level(),
            6
        );

        assert_eq!(FiniteNimber::from(&[0x55, 0, 0, 0]).to_ref().level(), 3);
        assert_eq!(FiniteNimber::from(&[0x55, 1, 0, 0]).to_ref().level(), 7);
        assert_eq!(FiniteNimber::from(&[0x55, 1, 1, 0]).to_ref().level(), 8);
        assert_eq!(FiniteNimber::from(&[0x55, 1, 1, 1]).to_ref().level(), 8);
        assert_eq!(
            FiniteNimber::from(&[0x55, 1, 1, 1, 1]).to_ref().level(),
            9
        );
        assert_eq!(
            FiniteNimber::from(&[1, 1, 1, 1, 1, 1, 1, 1])
                .to_ref()
                .level(),
            9
        );
        assert_eq!(
            FiniteNimber::from(&[1, 1, 1, 1, 1, 1, 1, 1, 1])
                .to_ref()
                .level(),
            10
        );
    }

    #[test]
    fn addition() {
        let a = FiniteNimber::from(0b1100);
        let b = FiniteNimber::from(0b1010);
        assert_eq!(a + b, FiniteNimber::from(0b0110));

        let a = FiniteNimber::from(&[1, 1, 0, 0]);
        let b = FiniteNimber::from(&[1, 0, 1, 0]);
        assert_eq!(a + b, FiniteNimber::from(&[0, 1, 1 /*, 0 */]));

        let a = FiniteNimber::from(&[1, 1, 1, 1]);
        let b = FiniteNimber::from(2);
        assert_eq!(a + b, FiniteNimber::from(&[3, 1, 1, 1]));
    }

    #[test]
    fn mul_by_h() {
        let a = FiniteNimber::from(1);
        assert_eq!(a.to_ref().mul_by_h(0), FiniteNimber::from(1));
        assert_eq!(a.to_ref().mul_by_h(1), FiniteNimber::from(2));
        assert_eq!(a.to_ref().mul_by_h(2), FiniteNimber::from(8));
        assert_eq!(a.to_ref().mul_by_h(3), FiniteNimber::from(128));
        assert_eq!(
            a.to_ref().mul_by_h(8),
            FiniteNimber::from(&[0, 0, 0, 0x8000000000000000])
        );

        assert_eq!(
            a.to_ref().mul_by_h(3).to_ref().mul_by_h(3),
            FiniteNimber::from(0xde)
        );
        assert_eq!(
            a.to_ref().mul_by_h(8).to_ref().mul_by_h(8),
            FiniteNimber::from(&[
                0xa92181714b010a1a,
                0x4a88a921e2208b6b,
                0xe3a92850a9218171,
                0xde4ae3a94a88a921
            ])
        );
    }

    #[test]
    fn mul() {
        assert_eq!(
            FiniteNimber::from(0xd) * FiniteNimber::from(0xb),
            FiniteNimber::from(0x5)
        );

        assert_eq!(
            FiniteNimber::from(0x123456789abcdef0) * FiniteNimber::from(0x76),
            FiniteNimber::from(0x84946b8cf1e11ef9)
        );

        assert_eq!(
            FiniteNimber::from(&[
                0x3f84d5b5b5470917,
                0xc0ac29b7c97c50dd,
                0xbe5466cf34e90c6c,
                0x452821e638d01377,
                0x082efa98ec4e6c89,
                0xa4093822299f31d0,
                0x13198a2e03707344,
                0x243f6a8885a308d3,
            ]) * FiniteNimber::from(&[
                0x4f7c7b5757f59584,
                0xda06c80abb1185eb,
                0xf4bf8d8d8c31d763,
                0x324e7738926cfbe5,
                0xa784d9045190cfef,
                0x62e7160f38b4da56,
                0xbf7158809cf4f3c7,
                0xb7e151628aed2a6a,
            ]),
            FiniteNimber::from(&[
                0x72890f944be121d8,
                0xe6429b02014feb2e,
                0x2454070e4408eff8,
                0x298c025ec5ac4190,
                0xba6895a109cfcf6d,
                0xcfb08c013e9dd19c,
                0x5eeb02eaf7ae9ea0,
                0x59cbe194a9599171,
            ])
        );
    }

    #[test]
    fn square() {
        assert_eq!(
            FiniteNimber::from(0x80).square(),
            FiniteNimber::from(0xde)
        );
        assert_eq!(
            FiniteNimber::from(0x8000).square(),
            FiniteNimber::from(0xde4a)
        );

        assert_eq!(
            FiniteNimber::from(&[
                0x3f84d5b5b5470917,
                0xc0ac29b7c97c50dd,
                0xbe5466cf34e90c6c,
                0x452821e638d01377,
                0x082efa98ec4e6c89,
                0xa4093822299f31d0,
                0x13198a2e03707344,
                0x243f6a8885a308d3,
            ])
            .square(),
            FiniteNimber::from(&[
                0xd0e74a0945d35342,
                0xfa91473032b5438a,
                0xe17bdc047300f99d,
                0x1a51cba4adb8ddb9,
                0xdc06e76f54f372c0,
                0xe74d7b7c65542a19,
                0xffe69bcb391c628b,
                0x32ae6e49dcd65156,
            ])
        );
    }

    #[test]
    fn div() {
        assert_eq!(
            FiniteNimber::from(1) / FiniteNimber::from(2),
            FiniteNimber::from(3)
        );

        assert_eq!(
            FiniteNimber::from(&[
                0x72890f944be121d8,
                0xe6429b02014feb2e,
                0x2454070e4408eff8,
                0x298c025ec5ac4190,
                0xba6895a109cfcf6d,
                0xcfb08c013e9dd19c,
                0x5eeb02eaf7ae9ea0,
                0x59cbe194a9599171,
            ]) / FiniteNimber::from(&[
                0x3f84d5b5b5470917,
                0xc0ac29b7c97c50dd,
                0xbe5466cf34e90c6c,
                0x452821e638d01377,
                0x082efa98ec4e6c89,
                0xa4093822299f31d0,
                0x13198a2e03707344,
                0x243f6a8885a308d3,
            ]),
            FiniteNimber::from(&[
                0x4f7c7b5757f59584,
                0xda06c80abb1185eb,
                0xf4bf8d8d8c31d763,
                0x324e7738926cfbe5,
                0xa784d9045190cfef,
                0x62e7160f38b4da56,
                0xbf7158809cf4f3c7,
                0xb7e151628aed2a6a,
            ])
        );
    }

    #[test]
    #[should_panic]
    fn div0() {
        // Deliberately divide by zero, and ensure we get the expected
        // panic.
        let _ = FiniteNimber::from(0x1234) / FiniteNimber::from(0);
    }

    #[test]
    fn sqrt() {
        assert_eq!(FiniteNimber::from(0xde).sqrt(), FiniteNimber::from(0x80));
        assert_eq!(
            FiniteNimber::from(0xde4a).sqrt(),
            FiniteNimber::from(0x8000)
        );

        assert_eq!(
            FiniteNimber::from(&[
                0xd0e74a0945d35342,
                0xfa91473032b5438a,
                0xe17bdc047300f99d,
                0x1a51cba4adb8ddb9,
                0xdc06e76f54f372c0,
                0xe74d7b7c65542a19,
                0xffe69bcb391c628b,
                0x32ae6e49dcd65156,
            ])
            .sqrt(),
            FiniteNimber::from(&[
                0x3f84d5b5b5470917,
                0xc0ac29b7c97c50dd,
                0xbe5466cf34e90c6c,
                0x452821e638d01377,
                0x082efa98ec4e6c89,
                0xa4093822299f31d0,
                0x13198a2e03707344,
                0x243f6a8885a308d3,
            ])
        );
    }

    #[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());

        // Degenerate case of quadratic3 with a=0, checking it returns None.
        assert_eq!(
            FiniteNimber::quadratic3(&FiniteNimber::from(0), &c, &c),
            None
        );

        // 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).unwrap();
        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));
    }
}