Skip to main content

binius_field/
transpose.rs

1// Copyright 2023-2025 Irreducible Inc.
2// Copyright 2026 The Binius Developers
3
4use binius_utils::checked_arithmetics::checked_log_2;
5
6use super::packed::PackedField;
7use crate::{BinaryField, ExtensionField, PackedExtension};
8
9/// Transpose square blocks of elements within packed field elements in place.
10///
11/// The input elements are interpreted as a rectangular matrix with height `n = 2^n` in row-major
12/// order. This matrix is interpreted as a vector of square matrices of field elements, and each
13/// square matrix is transposed in-place.
14///
15/// # Arguments
16///
17/// * `log_n`: The base-2 logarithm of the dimension of the n x n square matrix. Must be less than
18///   or equal to the base-2 logarithm of the packing width.
19/// * `elems`: The packed field elements, length is a power-of-two multiple of `1 << log_n`.
20///
21/// # Preconditions
22///
23/// * `log_n` must be at most `P::LOG_WIDTH`.
24/// * `elems.len()` must be a power of two and at least `2^log_n`.
25///
26/// A caller whose dimensions are compile-time constants should use the fixed-size form below.
27/// That form unrolls the butterfly and keeps the array in registers.
28pub fn transpose_square_blocks<P: PackedField>(log_n: usize, elems: &mut [P]) {
29	assert!(P::LOG_WIDTH >= log_n, "dimension n of square blocks must divide packing width");
30
31	let size = elems.len();
32	assert!(size.is_power_of_two(), "elems length must be a power of two, got {size}");
33	let log_size = checked_log_2(size);
34	assert!(
35		log_size >= log_n,
36		"elems must have length at least 2^log_n = {}, got {size}",
37		1 << log_n
38	);
39
40	let log_w = log_size - log_n;
41
42	// See Hacker's Delight, Section 7-3.
43	// https://dl.acm.org/doi/10.5555/2462741
44	for i in 0..log_n {
45		for j in 0..1 << (log_n - i - 1) {
46			for k in 0..1 << (log_w + i) {
47				let idx0 = (j << (log_w + i + 1)) | k;
48				let idx1 = idx0 | (1 << (log_w + i));
49
50				let v0 = elems[idx0];
51				let v1 = elems[idx1];
52				let (v0, v1) = v0.interleave(v1, i);
53				elems[idx0] = v0;
54				elems[idx1] = v1;
55			}
56		}
57	}
58}
59
60/// Transposes square blocks of scalars across a fixed-size array of packed elements, in place.
61///
62/// # Overview
63///
64/// The runtime-sized form in this module computes the same permutation.
65/// This form is for a caller whose block dimension and array length are both constants.
66///
67/// Constant sizes let the compiler unroll the butterfly.
68/// The whole array then stays in registers, which is what a caller in a hot loop wants.
69///
70/// # Algorithm
71///
72/// A butterfly network over `LOG_N` rounds, as in Hacker's Delight, Section 7-3.
73/// Round `i` interleaves element pairs `2^(log_w + i)` apart at block granularity `2^i`.
74///
75/// # Preconditions
76///
77/// All three are checked at compile time, so a violating instantiation fails to build:
78///
79/// * The array length must be a power of two.
80/// * The block dimension must not exceed the base-2 log of the array length.
81/// * The block dimension must not exceed the base-2 log of the packed width.
82pub fn transpose_square_blocks_array<P: PackedField, const LOG_N: usize, const S: usize>(
83	elems: &mut [P; S],
84) {
85	const {
86		assert!(LOG_N <= P::LOG_WIDTH, "LOG_N must not exceed the packed width");
87		assert!(LOG_N <= checked_log_2(S), "LOG_N must not exceed the array length");
88	}
89
90	let log_size = checked_log_2(S);
91
92	// Elements per block that stays contiguous through the butterfly.
93	let log_w = log_size - LOG_N;
94
95	for i in 0..LOG_N {
96		for j in 0..1 << (LOG_N - i - 1) {
97			for k in 0..1 << (log_w + i) {
98				// Partner elements for this round, one stride apart.
99				let idx0 = (j << (log_w + i + 1)) | k;
100				let idx1 = idx0 | (1 << (log_w + i));
101
102				// Interleaving at block granularity 2^i swaps the axes one bit at a time.
103				let (v0, v1) = elems[idx0].interleave(elems[idx1], i);
104				elems[idx0] = v0;
105				elems[idx1] = v1;
106			}
107		}
108	}
109}
110
111pub fn square_transforms_extension_field<F, FE>(values: &mut [FE])
112where
113	F: BinaryField,
114	FE: PackedExtension<F>,
115{
116	transpose_square_blocks(FE::Scalar::LOG_DEGREE, FE::cast_bases_mut(values));
117}
118
119#[cfg(test)]
120mod tests {
121	use std::array;
122
123	use proptest::prelude::*;
124	use rand::{SeedableRng, rngs::StdRng};
125
126	use super::*;
127	use crate::{PackedBinaryField64x1b, PackedBinaryField128x1b, PackedField, Random};
128
129	#[test]
130	fn test_transpose_square_blocks_128x1b() {
131		let mut elems = [
132			PackedBinaryField128x1b::from(0x00000000000000000000000000000000u128),
133			PackedBinaryField128x1b::from(0x00000000000000000000000000000000u128),
134			PackedBinaryField128x1b::from(0x00000000000000000000000000000000u128),
135			PackedBinaryField128x1b::from(0x00000000000000000000000000000000u128),
136			PackedBinaryField128x1b::from(0xffffffffffffffffffffffffffffffffu128),
137			PackedBinaryField128x1b::from(0xffffffffffffffffffffffffffffffffu128),
138			PackedBinaryField128x1b::from(0xffffffffffffffffffffffffffffffffu128),
139			PackedBinaryField128x1b::from(0xffffffffffffffffffffffffffffffffu128),
140		];
141		transpose_square_blocks(3, &mut elems);
142
143		let expected = [
144			PackedBinaryField128x1b::from(0xf0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0u128),
145			PackedBinaryField128x1b::from(0xf0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0u128),
146			PackedBinaryField128x1b::from(0xf0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0u128),
147			PackedBinaryField128x1b::from(0xf0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0u128),
148			PackedBinaryField128x1b::from(0xf0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0u128),
149			PackedBinaryField128x1b::from(0xf0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0u128),
150			PackedBinaryField128x1b::from(0xf0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0u128),
151			PackedBinaryField128x1b::from(0xf0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0u128),
152		];
153		assert_eq!(elems, expected);
154	}
155
156	#[test]
157	fn test_transpose_square_blocks_128x1b_multi_row() {
158		let mut elems = [
159			PackedBinaryField128x1b::from(0x00000000000000000000000000000000u128),
160			PackedBinaryField128x1b::from(0x00000000000000000000000000000000u128),
161			PackedBinaryField128x1b::from(0x00000000000000000000000000000000u128),
162			PackedBinaryField128x1b::from(0x00000000000000000000000000000000u128),
163			PackedBinaryField128x1b::from(0xffffffffffffffffffffffffffffffffu128),
164			PackedBinaryField128x1b::from(0xffffffffffffffffffffffffffffffffu128),
165			PackedBinaryField128x1b::from(0xffffffffffffffffffffffffffffffffu128),
166			PackedBinaryField128x1b::from(0xffffffffffffffffffffffffffffffffu128),
167		];
168		transpose_square_blocks(1, &mut elems);
169
170		let expected = [
171			PackedBinaryField128x1b::from(0xaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaau128),
172			PackedBinaryField128x1b::from(0xaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaau128),
173			PackedBinaryField128x1b::from(0xaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaau128),
174			PackedBinaryField128x1b::from(0xaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaau128),
175			PackedBinaryField128x1b::from(0xaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaau128),
176			PackedBinaryField128x1b::from(0xaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaau128),
177			PackedBinaryField128x1b::from(0xaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaau128),
178			PackedBinaryField128x1b::from(0xaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaau128),
179		];
180		assert_eq!(elems, expected);
181	}
182
183	// The fixed-size form exists only to unroll the loop, so it must compute exactly the
184	// permutation the runtime form computes. Pinning them equal is what justifies keeping both.
185	//
186	//     same input -> runtime form   -> A
187	//                -> fixed-size form -> B
188	//     A == B for every block dimension the array length admits
189	fn check_forms_agree<P, const LOG_N: usize, const S: usize>(seed: u64)
190	where
191		P: PackedField + Random,
192	{
193		let mut rng = StdRng::seed_from_u64(seed);
194
195		// Random lanes over many trials cover every one of the S * WIDTH scalar positions.
196		for _ in 0..100 {
197			let input: [P; S] = array::from_fn(|_| P::random(&mut rng));
198
199			let mut runtime = input;
200			transpose_square_blocks(LOG_N, &mut runtime);
201
202			let mut fixed = input;
203			transpose_square_blocks_array::<P, LOG_N, S>(&mut fixed);
204
205			assert_eq!(fixed, runtime, "forms disagree at LOG_N = {LOG_N}, S = {S}");
206		}
207	}
208
209	#[test]
210	fn fixed_size_form_agrees_with_runtime_form() {
211		// Cover both row widths the callers run at, and every block dimension each admits.
212		//
213		//     64 lanes  -> LOG_N up to 6, array length 8 admits up to 3
214		//     128 lanes -> LOG_N up to 7, array length 8 admits up to 3
215		check_forms_agree::<PackedBinaryField64x1b, 0, 8>(0);
216		check_forms_agree::<PackedBinaryField64x1b, 1, 8>(1);
217		check_forms_agree::<PackedBinaryField64x1b, 3, 8>(2);
218		check_forms_agree::<PackedBinaryField128x1b, 3, 8>(3);
219
220		// A block dimension equal to the array length exercises the widest butterfly.
221		check_forms_agree::<PackedBinaryField64x1b, 4, 16>(4);
222		check_forms_agree::<PackedBinaryField128x1b, 5, 32>(5);
223	}
224
225	#[test]
226	fn transpose_exchanges_element_axis_with_low_scalar_bits() {
227		let mut rng = StdRng::seed_from_u64(0);
228
229		// The permutation itself, stated directly rather than through either implementation.
230		// Splitting a scalar position into a high part and its low three bits:
231		//
232		//     input : element r, position 8i + j  =  value at (r, 8i + j)
233		//     output: element j, position 8i + t  =  value at (t, 8i + j)
234		//
235		// So the element index and the low three bits of the position trade places.
236		for _ in 0..100 {
237			let input: [PackedBinaryField128x1b; 8] =
238				array::from_fn(|_| PackedBinaryField128x1b::random(&mut rng));
239			let mut output = input;
240			transpose_square_blocks_array::<_, 3, 8>(&mut output);
241
242			// Read both sides as scalars, so the assertion is about positions and not underliers.
243			let scalars = |elems: &[PackedBinaryField128x1b; 8]| {
244				elems
245					.iter()
246					.map(|e| e.iter().collect::<Vec<_>>())
247					.collect::<Vec<_>>()
248			};
249			let before = scalars(&input);
250			let after = scalars(&output);
251
252			// High part of the position, which the permutation leaves alone.
253			for i in 0..PackedBinaryField128x1b::WIDTH / 8 {
254				// Element of the output, which is the low three bits of the input position.
255				for j in 0..8 {
256					// Element of the input, which becomes the low three bits of the output.
257					for t in 0..8 {
258						assert_eq!(
259							after[j][i * 8 + t],
260							before[t][i * 8 + j],
261							"i={i}, j={j}, t={t}"
262						);
263					}
264				}
265			}
266		}
267	}
268
269	proptest! {
270		#[test]
271		fn transpose_is_an_involution(values in prop::collection::vec(any::<u128>(), 8)) {
272			// Exchanging two axes twice restores the original layout.
273			// This holds for the fixed-size form on any input, so it is a property, not a case.
274			let input: [PackedBinaryField128x1b; 8] =
275				array::from_fn(|i| PackedBinaryField128x1b::from(values[i]));
276
277			let mut roundtrip = input;
278			transpose_square_blocks_array::<_, 3, 8>(&mut roundtrip);
279			transpose_square_blocks_array::<_, 3, 8>(&mut roundtrip);
280
281			prop_assert_eq!(roundtrip, input);
282		}
283	}
284}