Use the groestl-style T-table lookup for fast kupyna hashing (#778)

A lot of the galois-field ops were running on each and every column
load. I didn't do anything clever or original, aside from understanding
how the groestl implementation worked, and shamelessly copying as much
from that as possible. Once that cleaned out a big pile of intermediate
functions, I picked some low-hanging fruit with the new `const fn
compute_src_cols`.

Moving all the field xor ops into a groestl-esque T-table sees a
significant 8x speedup. One weird issue I hit was on the long-hashes,
LLVM stops automatically unrolling the loops once they're longer than 8
iterations, which was seeing a major dropoff in throughput. I could've
manually unrolled the whole thing, but chunking it into two separate
8-longs (for a 16-byte long chunk) gives me equivalent compiled code as
both will unroll it, even though it's a bit silly to read.

Benchmarks before this PR:
```
running 8 tests
test kupyna256_10    ... bench:         276.20 ns/iter (+/- 5.70) = 36 MB/s
test kupyna256_100   ... bench:       2,761.53 ns/iter (+/- 44.17) = 36 MB/s
test kupyna256_1000  ... bench:      27,618.72 ns/iter (+/- 235.44) = 36 MB/s
test kupyna256_10000 ... bench:     274,159.70 ns/iter (+/- 3,121.12) = 36 MB/s
test kupyna512_10    ... bench:         391.96 ns/iter (+/- 21.09) = 25 MB/s
test kupyna512_100   ... bench:       3,915.42 ns/iter (+/- 142.03) = 25 MB/s
test kupyna512_1000  ... bench:      39,039.71 ns/iter (+/- 182.78) = 25 MB/s
test kupyna512_10000 ... bench:     389,250.45 ns/iter (+/- 1,132.40) = 25 MB/s
```

After swapping the galois-field xor, subbytes loading, and shifting with
a quick lookup in the T-table:
```
running 8 tests
test kupyna256_10    ... bench:          38.17 ns/iter (+/- 3.66) = 263 MB/s
test kupyna256_100   ... bench:         370.85 ns/iter (+/- 5.36) = 270 MB/s
test kupyna256_1000  ... bench:       3,568.39 ns/iter (+/- 901.09) = 280 MB/s
test kupyna256_10000 ... bench:      36,356.67 ns/iter (+/- 1,409.35) = 275 MB/s
test kupyna512_10    ... bench:          48.18 ns/iter (+/- 0.31) = 208 MB/s
test kupyna512_100   ... bench:         462.74 ns/iter (+/- 21.66) = 216 MB/s
test kupyna512_1000  ... bench:       4,738.49 ns/iter (+/- 250.94) = 211 MB/s
test kupyna512_10000 ... bench:      47,398.92 ns/iter (+/- 1,137.88) = 210 MB/s
```

Then the `const fn computer_src_cols` instead of computing at runtime:
```
running 8 tests
test kupyna256_10    ... bench:          33.61 ns/iter (+/- 0.35) = 303 MB/s
test kupyna256_100   ... bench:         338.50 ns/iter (+/- 24.91) = 295 MB/s
test kupyna256_1000  ... bench:       3,369.35 ns/iter (+/- 85.24) = 296 MB/s
test kupyna256_10000 ... bench:      32,791.09 ns/iter (+/- 315.57) = 304 MB/s
test kupyna512_10    ... bench:          45.69 ns/iter (+/- 2.19) = 222 MB/s
test kupyna512_100   ... bench:         442.11 ns/iter (+/- 26.73) = 226 MB/s
test kupyna512_1000  ... bench:       4,370.52 ns/iter (+/- 55.67) = 228 MB/s
test kupyna512_10000 ... bench:      43,252.35 ns/iter (+/- 252.63) = 231 MB/s
```

Overall a ~8.5x speedup for the 256 hashes, and >9x for the 512!

Author: jkoudys

kupyna/src/lib.rs

@@ -15,6 +15,7 @@ pub mod block_api;
 mod consts;
 mod long;
 mod short;
+mod table;
 pub(crate) mod utils;
 
 use digest::consts::{U28, U32, U48, U64};

kupyna/src/long.rs

@@ -1,53 +1,88 @@
-use crate::utils::{
-    add_constant_plus, add_constant_xor, apply_s_box, mix_columns, read_u64s_be, xor,
-};
-use core::array;
+use crate::table::TABLE;
+use crate::utils::{read_u64s_be, xor};
 
 pub(crate) const COLS: usize = 16;
-const ROUNDS: u64 = 14;
+const ROUNDS: usize = 14;
+
+// Bit shift amounts to extract each byte from a u64
+const BYTE_SHIFTS: [usize; 8] = [56, 48, 40, 32, 24, 16, 8, 0];
+
+// ShiftRows offsets for long variant: rows 0-6 shift by index, row 7 shifts by 11
+const SHIFTS: [usize; 8] = [0, 1, 2, 3, 4, 5, 6, 11];
+
+// Precomputed source columns: SRC_COLS[col][row] = (col + COLS - SHIFTS[row]) % COLS
+const fn compute_src_cols() -> [[usize; 8]; COLS] {
+    let mut result = [[0; 8]; COLS];
+    let mut col = 0;
+    while col < COLS {
+        let mut row = 0;
+        while row < 8 {
+            result[col][row] = (col + COLS - SHIFTS[row]) % COLS;
+            row += 1;
+        }
+        col += 1;
+    }
+    result
+}
+const SRC_COLS: [[usize; 8]; COLS] = compute_src_cols();
 
 pub(crate) fn compress(prev_vector: &mut [u64; COLS], message_block: &[u8; 128]) {
-    // Convert message block from u8 to u64 (column-major order as per paper)
     let message_u64 = read_u64s_be::<128, COLS>(message_block);
     let m_xor_p = xor(*prev_vector, message_u64);
     let t_xor_mp = t_xor_l(m_xor_p);
     let t_plus_m = t_plus_l(message_u64);
     *prev_vector = xor(xor(t_xor_mp, t_plus_m), *prev_vector);
 }
 
-fn t_plus_l(state: [u64; COLS]) -> [u64; COLS] {
-    let mut state = state;
-    for nu in 0..ROUNDS {
-        add_constant_plus(&mut state, nu as usize);
-        apply_s_box(&mut state);
-        state = rotate_rows(state);
-        mix_columns(&mut state);
+/// Compute one output column using T-table lookups
+#[inline(always)]
+fn column(x: &[u64; COLS], col: usize) -> u64 {
+    let mut t = 0u64;
+    for row in 0..8 {
+        let byte = ((x[SRC_COLS[col][row]] >> BYTE_SHIFTS[row]) & 0xFF) as usize;
+        t ^= TABLE[row][byte];
     }
-    state
+    t
 }
 
-fn rotate_rows(state: [u64; COLS]) -> [u64; COLS] {
-    //shift amounts for each row (0-6: row index, 7: special case = 11)
-    const SHIFTS: [usize; 8] = [0, 1, 2, 3, 4, 5, 6, 11];
-
-    array::from_fn(|col| {
-        let rotated_bytes = array::from_fn(|row| {
-            let shift = SHIFTS[row];
-            let src_col = (col + COLS - shift) % COLS;
-            let src_bytes = state[src_col].to_be_bytes();
-            src_bytes[row]
-        });
-        u64::from_be_bytes(rotated_bytes)
-    })
+fn t_plus_l(mut state: [u64; COLS]) -> [u64; COLS] {
+    for round in 0..ROUNDS {
+        // AddConstantPlus
+        for (i, word) in state.iter_mut().enumerate() {
+            *word = word
+                .swap_bytes()
+                .wrapping_add(
+                    0x00F0F0F0F0F0F0F3u64 ^ (((((COLS - i - 1) * 0x10) ^ round) as u64) << 56),
+                )
+                .swap_bytes();
+        }
+        // Fused SubBytes + ShiftRows + MixColumns via T-tables
+        let prev = state;
+        for (col, slot) in state[..8].iter_mut().enumerate() {
+            *slot = column(&prev, col);
+        }
+        for (col, slot) in state[8..].iter_mut().enumerate() {
+            *slot = column(&prev, col + 8);
+        }
+    }
+    state
 }
 
-pub(crate) fn t_xor_l(state: [u64; COLS]) -> [u64; COLS] {
-    let mut state = state;
-    for nu in 0..ROUNDS {
-        add_constant_xor(&mut state, nu as usize);
-        apply_s_box(&mut state);
-        state = rotate_rows(state);
-        mix_columns(&mut state);
+pub(crate) fn t_xor_l(mut state: [u64; COLS]) -> [u64; COLS] {
+    for round in 0..ROUNDS {
+        // AddConstantXor
+        for (i, word) in state.iter_mut().enumerate() {
+            let constant = ((i * 0x10) ^ round) as u64;
+            *word ^= constant << 56;
+        }
+        // Fused SubBytes + ShiftRows + MixColumns via T-tables
+        let prev = state;
+        for (col, slot) in state[..8].iter_mut().enumerate() {
+            *slot = column(&prev, col);
+        }
+        for (col, slot) in state[8..].iter_mut().enumerate() {
+            *slot = column(&prev, col + 8);
+        }
     }
     state
 }

kupyna/src/short.rs

@@ -1,53 +1,80 @@
-use crate::utils::{
-    add_constant_plus, add_constant_xor, apply_s_box, mix_columns, read_u64s_be, xor,
-};
-use core::array;
+use crate::table::TABLE;
+use crate::utils::{read_u64s_be, xor};
 
 pub(crate) const COLS: usize = 8;
-const ROUNDS: u64 = 10;
+const ROUNDS: usize = 10;
+
+// Bit shift amounts to extract each byte from a u64
+const BYTE_SHIFTS: [usize; 8] = [56, 48, 40, 32, 24, 16, 8, 0];
+
+// Precomputed source columns: SRC_COLS[col][row] = (col + COLS - row) % COLS
+// ShiftRows for short variant: row i shifts by i positions
+const fn compute_src_cols() -> [[usize; 8]; COLS] {
+    let mut result = [[0; 8]; COLS];
+    let mut col = 0;
+    while col < COLS {
+        let mut row = 0;
+        while row < 8 {
+            result[col][row] = (col + COLS - row) % COLS;
+            row += 1;
+        }
+        col += 1;
+    }
+    result
+}
+const SRC_COLS: [[usize; 8]; COLS] = compute_src_cols();
 
 pub(crate) fn compress(prev_vector: &mut [u64; COLS], message_block: &[u8; 64]) {
-    // Convert message block from u8 to u64 (column-major order as per paper)
     let message_u64 = read_u64s_be::<64, COLS>(message_block);
     let m_xor_p = xor(*prev_vector, message_u64);
     let t_xor_mp = t_xor_l(m_xor_p);
     let t_plus_m = t_plus_l(message_u64);
     *prev_vector = xor(xor(t_xor_mp, t_plus_m), *prev_vector);
 }
 
-fn t_plus_l(state: [u64; COLS]) -> [u64; COLS] {
-    let mut state = state;
-    for nu in 0..ROUNDS {
-        add_constant_plus(&mut state, nu as usize);
-        apply_s_box(&mut state);
-        state = rotate_rows(state);
-        mix_columns(&mut state);
+/// Compute one output column using T-table lookups
+#[inline(always)]
+fn column(x: &[u64; COLS], col: usize) -> u64 {
+    let mut t = 0u64;
+    for row in 0..8 {
+        let byte = ((x[SRC_COLS[col][row]] >> BYTE_SHIFTS[row]) & 0xFF) as usize;
+        t ^= TABLE[row][byte];
     }
-    state
+    t
 }
 
-fn rotate_rows(state: [u64; COLS]) -> [u64; COLS] {
-    //shift amounts for each row (0-6: row index, 7: special case)
-    const SHIFTS: [usize; 8] = [0, 1, 2, 3, 4, 5, 6, 7];
-
-    array::from_fn(|col| {
-        let rotated_bytes = array::from_fn(|row| {
-            let shift = SHIFTS[row];
-            let src_col = (col + COLS - shift) % COLS;
-            let src_bytes = state[src_col].to_be_bytes();
-            src_bytes[row]
-        });
-        u64::from_be_bytes(rotated_bytes)
-    })
+fn t_plus_l(mut state: [u64; COLS]) -> [u64; COLS] {
+    for round in 0..ROUNDS {
+        // AddConstantPlus
+        for (i, word) in state.iter_mut().enumerate() {
+            *word = word
+                .swap_bytes()
+                .wrapping_add(
+                    0x00F0F0F0F0F0F0F3u64 ^ (((((COLS - i - 1) * 0x10) ^ round) as u64) << 56),
+                )
+                .swap_bytes();
+        }
+        // Fused SubBytes + ShiftRows + MixColumns via T-tables
+        let prev = state;
+        for (col, slot) in state.iter_mut().enumerate() {
+            *slot = column(&prev, col);
+        }
+    }
+    state
 }
 
-pub(crate) fn t_xor_l(state: [u64; COLS]) -> [u64; COLS] {
-    let mut state = state;
-    for nu in 0..ROUNDS {
-        add_constant_xor(&mut state, nu as usize);
-        apply_s_box(&mut state);
-        state = rotate_rows(state);
-        mix_columns(&mut state);
+pub(crate) fn t_xor_l(mut state: [u64; COLS]) -> [u64; COLS] {
+    for round in 0..ROUNDS {
+        // AddConstantXor
+        for (i, word) in state.iter_mut().enumerate() {
+            let constant = ((i * 0x10) ^ round) as u64;
+            *word ^= constant << 56;
+        }
+        // Fused SubBytes + ShiftRows + MixColumns via T-tables
+        let prev = state;
+        for (col, slot) in state.iter_mut().enumerate() {
+            *slot = column(&prev, col);
+        }
     }
     state
 }

kupyna/src/table.rs

@@ -0,0 +1,65 @@
+use crate::consts::{MDS_MATRIX, SBOXES};
+
+/// GF(2^8) multiplication with reduction polynomial x^8 + x^4 + x^3 + x^2 + 1
+const fn gf_multiply(mut x: u8, mut y: u8) -> u8 {
+    const REDUCTION_POLYNOMIAL: u8 = 0x1d; // x^4 + x^3 + x^2 + 1
+
+    let mut r = 0u8;
+    let mut i = 0;
+    while i < 8 {
+        if y & 1 == 1 {
+            r ^= x;
+        }
+        let hbit = x & 0x80;
+        x <<= 1;
+        if hbit != 0 {
+            x ^= REDUCTION_POLYNOMIAL;
+        }
+        y >>= 1;
+        i += 1;
+    }
+    r
+}
+
+/// Generate T-tables that fuse SubBytes + MixColumns
+///
+/// TABLE[row][byte] gives the contribution to an output column when input byte
+/// at position `row` has value `byte`, after applying S-box and MDS multiplication.
+const fn generate_t_table() -> [[u64; 256]; 8] {
+    let mut table = [[0u64; 256]; 8];
+
+    let mut row = 0;
+    while row < 8 {
+        let mut byte = 0;
+        while byte < 256 {
+            // Apply S-box for this row position (S-boxes cycle with period 4)
+            let s = SBOXES[row % 4][byte];
+
+            // Compute contribution to each output row via MDS multiplication
+            let mut out = [0u8; 8];
+            let mut out_row = 0;
+            while out_row < 8 {
+                // Extract MDS coefficient: MDS_MATRIX[out_row] byte at position `row`
+                let mds_coef = (MDS_MATRIX[out_row] >> (8 * (7 - row))) as u8;
+                out[out_row] = gf_multiply(mds_coef, s);
+                out_row += 1;
+            }
+
+            // Pack into u64 (big-endian)
+            table[row][byte] = ((out[0] as u64) << 56)
+                | ((out[1] as u64) << 48)
+                | ((out[2] as u64) << 40)
+                | ((out[3] as u64) << 32)
+                | ((out[4] as u64) << 24)
+                | ((out[5] as u64) << 16)
+                | ((out[6] as u64) << 8)
+                | (out[7] as u64);
+
+            byte += 1;
+        }
+        row += 1;
+    }
+    table
+}
+
+pub(crate) static TABLE: [[u64; 256]; 8] = generate_t_table();