p3-binary-field 0.8.0

//! The Cantor basis of the tower: `v_0 = 1` and `v_i² + v_i = v_{i−1}`, taking the root
//! whose lowest bit is clear.
//!
//! `v_0, …, v_{2^k−1}` is an `F_2`-basis of `GF(2^(2^k))`, so the first `2^ℓ` vectors span the
//! additive NTT domain `S_ℓ` and the domain a level sees is the domain every wider level sees.

/// The Cantor basis of `GF(2^128)` in the tower representation, lowest vector first.
///
/// Regenerated by the `regenerate_cantor_basis` test in this module.
pub(crate) const CANTOR_BASIS_128: [u128; 128] = [
    0x00000000000000000000000000000001, // v_0
    0x00000000000000000000000000000002, // v_1
    0x0000000000000000000000000000000c, // v_2
    0x00000000000000000000000000000006, // v_3
    0x00000000000000000000000000000066, // v_4
    0x000000000000000000000000000000be, // v_5
    0x00000000000000000000000000000036, // v_6
    0x000000000000000000000000000000ac, // v_7
    0x00000000000000000000000000001444, // v_8
    0x000000000000000000000000000028f6, // v_9
    0x0000000000000000000000000000c70e, // v_10
    0x000000000000000000000000000061b0, // v_11
    0x000000000000000000000000000067e6, // v_12
    0x0000000000000000000000000000e4e8, // v_13
    0x000000000000000000000000000062f4, // v_14
    0x0000000000000000000000000000cd68, // v_15
    0x0000000000000000000000000110a454, // v_16
    0x0000000000000000000000000220108c, // v_17
    0x0000000000000000000000000cc04ae6, // v_18
    0x00000000000000000000000006608fe0, // v_19
    0x00000000000000000000000066764360, // v_20
    0x000000000000000000000000becba7e2, // v_21
    0x00000000000000000000000036a31b50, // v_22
    0x000000000000000000000000acaabff8, // v_23
    0x00000000000000000000000044c0ca14, // v_24
    0x000000000000000000000000f6f7191e, // v_25
    0x0000000000000000000000000e278cfa, // v_26
    0x000000000000000000000000b04a5fbe, // v_27
    0x000000000000000000000000e6f91854, // v_28
    0x000000000000000000000000e89aa24c, // v_29
    0x000000000000000000000000f49dac5c, // v_30
    0x000000000000000000000000685bb70c, // v_31
    0x0000000000000000000101003bc34582, // v_32
    0x000000000000000000020200396bb03a, // v_33
    0x0000000000000000000c0c00d0b0e8b8, // v_34
    0x0000000000000000000606004e9ad958, // v_35
    0x00000000000000000066660044f0b874, // v_36
    0x000000000000000000bebe00aec9fbde, // v_37
    0x000000000000000000363600635f3dec, // v_38
    0x000000000000000000acac00e7fcaec8, // v_39
    0x000000000000000014444514b62553ae, // v_40
    0x000000000000000028f6f4280229b55a, // v_41
    0x0000000000000000c70e02c700c6a3ac, // v_42
    0x000000000000000061b0b66163f636ca, // v_43
    0x000000000000000067e68067519e8e0c, // v_44
    0x0000000000000000e4e856e4fe5a106c, // v_45
    0x000000000000000062f4c26257122bf8, // v_46
    0x0000000000000000cd68c4cd741f31fc, // v_47
    0x0000000000000000a4547fb4e2455764, // v_48
    0x0000000000000000108ccf306a375460, // v_49
    0x00000000000000004ae6de8a2a81316e, // v_50
    0x00000000000000008fe0feefcd4bc340, // v_51
    0x00000000000000004360a235b3dc7fe0, // v_52
    0x0000000000000000a7e2466cf2223352, // v_53
    0x00000000000000001b5097b817c26446, // v_54
    0x0000000000000000bff88f155221e120, // v_55
    0x0000000000000000ca148c0a63a8ee6a, // v_56
    0x0000000000000000191e39ee81f99646, // v_57
    0x00000000000000008cfadcabc3c3aeee, // v_58
    0x00000000000000005fbe2b152121db4a, // v_59
    0x0000000000000000185473e1430fdc24, // v_60
    0x0000000000000000a24cee38da4d8828, // v_61
    0x0000000000000000ac5c0231925591ac, // v_62
    0x0000000000000000b70c3fecc1c7c4cc, // v_63
    0x000000010001000051ec1382cf362d58, // v_64
    0x0000000200020000988d9fd3a2bef392, // v_65
    0x0000000c000c00002f822c5f61945b52, // v_66
    0x0000000600060000b8f0134234d14098, // v_67
    0x0000006600660000df36d959956bdac0, // v_68
    0x000000be00be00001f2c8d1d900e9b18, // v_69
    0x00000036003600005fe999a3debb9ea6, // v_70
    0x000000ac00ac0000637b1635bdb090a6, // v_71
    0x0000144414440000e3487406d92c5496, // v_72
    0x000028f628f600008db07c530db5036a, // v_73
    0x0000c70ec70e00002e5834cfbf6bb772, // v_74
    0x000061b061b00000d8a7c9649277163c, // v_75
    0x000067e667e60000aaef92cdbf75be70, // v_76
    0x0000e4e8e4e800005336c5936a5e98b8, // v_77
    0x000062f462f40000520751dd3dae2c88, // v_78
    0x0000cd68cd68000043ec54073b43381c, // v_79
    0x0110a454a4550110393fb5bfb75f91c0, // v_80
    0x0220108c108e02205d095f01c84ac1de, // v_81
    0x0cc04ae64aea0cc0f72c19c6d2dc3176, // v_82
    0x06608fe08fe60660139b69c2aa148fc0, // v_83
    0x667643604306667624635608ad1d04a8, // v_84
    0xbecba7e2a75cbecb36d02fddc0ce697c, // v_85
    0x36a31b501b6636a3d3dea74e02b99a0c, // v_86
    0xacaabff8bf54acaae95d8634b11339fe, // v_87
    0x44c0ca14de5044c05d11c1d0ef2f5148, // v_88
    0xf6f7191e31e8f6f72b9f2efbb8c15e7c, // v_89
    0x0e278cfa4bf40e27448196416b37c3e2, // v_90
    0xb04a5fbe3e0eb04a8ec28859c3484036, // v_91
    0xe6f918547fb2e6f9a23a62337d1227ac, // v_92
    0xe89aa24c46a4e89a4d60005a92307dc8, // v_93
    0xf49dac5ccea8f49dd1917996f644f246, // v_94
    0x685bb70c7a64685b5685858ec16125ce, // v_95
    0x3bc34582f5b83ac34f7d352c22a85848, // v_96
    0x396bb03a88013b6bfa8e5b530962b4ec, // v_97
    0xd0b0e8b86564dcb09c4967490a0eaf20, // v_98
    0x4e9ad9583710489a5782d8adc3ed609c, // v_99
    0x44f0b8749c5622f097ff66566c81a23a, // v_100
    0xaec9fbdeb8ce10c9560d64b18f1ea1e4, // v_101
    0x635f3dec44b9555f0d3718d7b4efd9f4, // v_102
    0xe7fcaec8dc834bfc232edf18ad02d7cc, // v_103
    0xb62553ae295cf3312d288469082f3e6c, // v_104
    0x0229b55a94aef601f511844e1cdd037c, // v_105
    0x00c6a3aca2a20201e213bf6d8bd2a1cc, // v_106
    0x63f636ca8719d597dc7886c0904bf926, // v_107
    0x519e8e0cb2bbd1f9fa3035797f5d6f52, // v_108
    0xfe5a106cf17aa8be13062f9ad892a87c, // v_109
    0x57122bf8fe5b9570162c56709eceba34, // v_110
    0x741f31fcf4e0b0d25386f18c7587e27c, // v_111
    0xe245576468d29df19413bb1e4cc9c15a, // v_112
    0x6a375460c586a507598b9ef6cf61f8a4, // v_113
    0x2a81316ed8a2f40b70ae9ec84dfeb48c, // v_114
    0xcd4bc340ab6d33a4566d03b99de7c522, // v_115
    0xb3dc7fe0fb3311e93163772c0d4592da, // v_116
    0xf22233522942b44e1bc51d5f4e30af34, // v_117
    0x17c264468c01807a614eeccd14f35232, // v_118
    0x5221e1208a8add3417755bb752221124, // v_119
    0x63a8ee6aaa1defa259bbf63fde0852da, // v_120
    0x81f996468ed9b817c15dddc5c448b84a, // v_121
    0xc3c3aeeeadd71f687b7440d8015c79a8, // v_122
    0x2121db4a37d50a345c04f31a67e0ea34, // v_123
    0x430fdc246f3330eef0a8ff5700db4fec, // v_124
    0xda4d8828fabe34758244f65913a656c0, // v_125
    0x925591ace16290644c9c0c5214ba545c, // v_126
    0xc1c7c4ccd801fe2b6b65f1fbd173938c, // v_127
];

#[cfg(test)]
mod tests {
    extern crate std;

    use std::println;

    use p3_field::PrimeCharacteristicRing;

    use crate::tower::TowerLevel;
    use crate::{
        BinaryField2, BinaryField4, BinaryField8, BinaryField16, BinaryField32, BinaryField64,
        BinaryField128, Gf2, tables,
    };

    /// `v_0 = 1` and `v_i² + v_i = v_{i−1}`, checked in each level's own arithmetic.
    fn check_recurrence<F: TowerLevel>(count: usize) {
        assert_eq!(F::cantor_basis(0), F::ONE);
        for i in 1..count {
            let v = F::cantor_basis(i);
            assert_eq!(v.square() + v, F::cantor_basis(i - 1), "level bit {i}");
        }
    }

    /// `v_0, …, v_{2^k−1}` spans `GF(2^(2^k))`, so `v_i` fits in `2^⌈log2(i+1)⌉` bits, and
    /// every vector past `v_0` is the root of `x² + x = v_{i−1}` whose lowest bit is clear.
    fn check_nesting<F: TowerLevel>(count: usize)
    where
        F::Repr: Into<u128>,
    {
        for i in 0..count {
            let width = (i + 1).next_power_of_two();
            let repr: u128 = F::cantor_basis(i).to_repr().into();
            assert!(
                width >= 128 || repr < 1u128 << width,
                "v_{i} escapes GF(2^{width})"
            );
            assert_eq!(
                repr & 1,
                u128::from(i == 0),
                "v_{i} is not the even representation"
            );
        }
    }

    /// The first `count` vectors are linearly independent over `F_2`: reducing each against
    /// the pivots already found must never drive it to zero.
    fn check_independent<F: TowerLevel>(count: usize)
    where
        F::Repr: Into<u128>,
    {
        // `pivots[b]` holds a vector whose highest set bit is `b`, or zero if there is none.
        let mut pivots = [0u128; 128];
        for i in 0..count {
            let mut v: u128 = F::cantor_basis(i).to_repr().into();
            while v != 0 {
                let high = 127 - v.leading_zeros() as usize;
                if pivots[high] == 0 {
                    pivots[high] = v;
                    break;
                }
                v ^= pivots[high];
            }
            assert_ne!(v, 0, "v_{i} depends on the earlier vectors");
        }
    }

    #[test]
    fn recurrence_holds_at_every_level() {
        check_recurrence::<Gf2>(1);
        check_recurrence::<BinaryField2>(2);
        check_recurrence::<BinaryField4>(4);
        check_recurrence::<BinaryField8>(8);
        check_recurrence::<BinaryField16>(16);
        check_recurrence::<BinaryField32>(32);
        check_recurrence::<BinaryField64>(64);
        check_recurrence::<BinaryField128>(128);
    }

    #[test]
    fn basis_is_nested_and_independent() {
        check_nesting::<BinaryField2>(2);
        check_nesting::<BinaryField4>(4);
        check_nesting::<BinaryField8>(8);
        check_nesting::<BinaryField16>(16);
        check_nesting::<BinaryField32>(32);
        check_nesting::<BinaryField64>(64);
        check_nesting::<BinaryField128>(128);

        check_independent::<BinaryField2>(2);
        check_independent::<BinaryField4>(4);
        check_independent::<BinaryField8>(8);
        check_independent::<BinaryField16>(16);
        check_independent::<BinaryField32>(32);
        check_independent::<BinaryField64>(64);
        check_independent::<BinaryField128>(128);
    }

    /// A handful of the smallest vectors are pinned by literal value, independent of the table
    /// itself: `v_0 = 1` and `v_1 = 0b10` are the generator `X_0` of `GF(4)` over `GF(2)` by
    /// definition, and `v_2, v_3` are checked against the byte-level `GF(2^8)` arithmetic tables
    /// (`tables::square`), a code path the tower's own recursive squaring never exercises. A
    /// hand-edited table that flipped a root would still satisfy `check_recurrence` and
    /// `check_nesting` — both only check the table against itself — so this guards against
    /// exactly that: it fails if the table changes without the byte tables changing too.
    #[test]
    fn smallest_vectors_are_pinned() {
        assert_eq!(BinaryField128::cantor_basis(0).to_repr(), 1);
        assert_eq!(BinaryField128::cantor_basis(1).to_repr(), 2);

        let v2 = u8::try_from(BinaryField128::cantor_basis(2).to_repr()).unwrap();
        let v3 = u8::try_from(BinaryField128::cantor_basis(3).to_repr()).unwrap();
        assert_eq!(tables::square(v2) ^ v2, 2, "v_2^2 + v_2 = v_1");
        assert_eq!(tables::square(v3) ^ v3, v2, "v_3^2 + v_3 = v_2");
        assert_eq!(v2 & 1, 0, "v_2 is not the even representation");
        assert_eq!(v3 & 1, 0, "v_3 is not the even representation");
    }

    #[test]
    #[should_panic(expected = "index out of range")]
    fn out_of_range_index_panics() {
        let _ = BinaryField8::cantor_basis(8);
    }

    /// Regenerates `CANTOR_BASIS_128`. The map `x ↦ x² + x` is `F_2`-linear on
    /// `GF(2^128)` with kernel `{0, 1}`, so `v_i` is recovered by solving a 128×128 `F_2`
    /// system and taking the solution whose lowest bit is clear. Run with
    /// `cargo test -p p3-binary-field -- --ignored --nocapture regenerate_cantor_basis`
    /// and paste the printed table into `CANTOR_BASIS_128`.
    #[test]
    #[ignore]
    fn regenerate_cantor_basis() {
        // Column `b` of the Artin–Schreier matrix is `e_b² + e_b` for the bit-basis vector `e_b`.
        let columns: [u128; 128] = core::array::from_fn(|b| {
            let e = BinaryField128::from_repr(1u128 << b);
            (e.square() + e).to_repr()
        });

        let mut basis = [0u128; 128];
        basis[0] = 1;
        for i in 1..128 {
            basis[i] = solve_artin_schreier(&columns, basis[i - 1]);

            // Never print a table the recurrence does not satisfy.
            let v = BinaryField128::from_repr(basis[i]);
            assert_eq!(v.square() + v, BinaryField128::from_repr(basis[i - 1]));
            assert_eq!(basis[i] & 1, 0, "v_{i} is not the even representation");
        }

        for (i, v) in basis.iter().enumerate() {
            println!("    0x{v:032x}, // v_{i}");
        }
    }

    /// Solves `M x = target` over `F_2`, where `M` has the given columns, and returns the
    /// solution whose lowest bit is clear.
    ///
    /// Gauss-Jordan on the columns: `image[c]` is the partially reduced image of the unit
    /// vector `c` and `preimage[c]` the combination of unit vectors that produced it, so
    /// `M · preimage[c] == image[c]` holds throughout. The invariant
    /// `M · solution + rhs == target` is preserved by every step, so `rhs == 0` at the end
    /// means `solution` is a preimage of `target`.
    fn solve_artin_schreier(columns: &[u128; 128], target: u128) -> u128 {
        let mut image = *columns;
        let mut preimage: [u128; 128] = core::array::from_fn(|c| 1u128 << c);

        let mut rhs = target;
        let mut solution = 0u128;
        let mut consumed = 0u128;

        for bit in 0..128 {
            // Take any column not yet used as a pivot that still has this bit set.
            let Some(pivot) =
                (0..128).find(|&c| consumed >> c & 1 == 0 && image[c] >> bit & 1 == 1)
            else {
                continue;
            };
            consumed |= 1 << pivot;

            // Clear the bit everywhere else, so no later step can set it again.
            for c in 0..128 {
                if c != pivot && image[c] >> bit & 1 == 1 {
                    image[c] ^= image[pivot];
                    preimage[c] ^= preimage[pivot];
                }
            }
            if rhs >> bit & 1 == 1 {
                rhs ^= image[pivot];
                solution ^= preimage[pivot];
            }
        }

        assert_eq!(rhs, 0, "target is outside the image of x -> x^2 + x");

        // The kernel is `{0, 1}`, so exactly one of the two solutions has bit 0 clear.
        if solution & 1 == 1 {
            solution ^ 1
        } else {
            solution
        }
    }
}