binius_field/arch/portable/arithmetic/
itoh_tsujii.rs

1// Copyright 2026 The Binius Developers
2
3//! Itoh-Tsujii inversion for the GHASH field `GF(2^128)`.
4//!
5//! For a non-zero `x`, the inverse is `x^(2^128 - 2) = (x^(2^127 - 1))^2`. The exponent `2^127 - 1`
6//! is built up with an addition chain on the powers `beta_k := x^(2^k - 1)`, using the identity
7//!
8//! ```text
9//! beta_{a+b} = (beta_a)^(2^b) * beta_b.
10//! ```
11//!
12//! Squaring `beta_a` repeatedly `b` times (the `x -> x^(2^b)` power map) is an `F_2`-linear
13//! transformation. We precompute each power map as a [`BytewiseLookupTransformation`] (the [Method
14//! of Four Russians]), wrapped into a `Ghash128b -> Ghash128b` transform, and hold them in a
15//! process-wide [`LazyLock`] so the tables are computed once and shared read-only across all
16//! threads.
17//!
18//! [Method of Four Russians]: <https://en.wikipedia.org/wiki/Method_of_Four_Russians>
19
20use std::{array, iter, ops::Mul, sync::LazyLock};
21
22use bytemuck::TransparentWrapper;
23
24use crate::{
25	BinaryField1b, Divisible, ExtensionField, Ghash128b,
26	arch::M128,
27	arithmetic_traits::{InvertOrZero, Square},
28	linear_transformation::{
29		BytewiseLookupTransformation, BytewiseLookupTransformationFactory,
30		InputWrappingTransformationFactory, LinearTransformationFactory,
31		OutputWrappingTransformationFactory, Transformation, WrappingTransformation,
32	},
33};
34
35/// Number of bits in a GHASH element.
36const FIELD_BITS: usize = 128;
37
38/// A precomputed `x -> x^(2^n)` power map as a byte-lookup transform on `Ghash128b`.
39///
40/// The underlying [`BytewiseLookupTransformation`] operates on the underlier `M128`; the input and
41/// output wrappers lift it to a `Ghash128b -> Ghash128b` transform.
42type GhashPowerMap =
43	WrappingTransformation<BytewiseLookupTransformation<M128, M128>, Ghash128b, Ghash128b>;
44
45/// The power maps needed by the Itoh-Tsujii addition chain for the GHASH field.
46///
47/// Each field holds the transform for one power map `x -> x^(2^n)`, for the values of `n` that
48/// appear in the chain (`pow_2_7` is reused for both the `7 -> 14` and `56 -> 63` steps).
49struct GhashPowerMapTables {
50	pow_2_3: GhashPowerMap,
51	pow_2_7: GhashPowerMap,
52	pow_2_14: GhashPowerMap,
53	pow_2_28: GhashPowerMap,
54	pow_2_63: GhashPowerMap,
55}
56
57impl GhashPowerMapTables {
58	fn new() -> Self {
59		Self {
60			pow_2_3: compute_power_map_transform(3),
61			pow_2_7: compute_power_map_transform(7),
62			pow_2_14: compute_power_map_transform(14),
63			pow_2_28: compute_power_map_transform(28),
64			pow_2_63: compute_power_map_transform(63),
65		}
66	}
67}
68
69static GHASH_POWER_MAP_TABLES: LazyLock<GhashPowerMapTables> =
70	LazyLock::new(GhashPowerMapTables::new);
71
72/// Build the byte-lookup transform for the power map `x -> x^(2^n)` over `Ghash128b`.
73///
74/// The power map is the `F_2`-linear transformation whose matrix has one column per input bit
75/// (`compute_power_map_matrix`). [`BytewiseLookupTransformation`] turns that column set into
76/// byte-indexed lookup tables; the input/output wrappers make it accept and return `Ghash128b`.
77fn compute_power_map_transform(n: usize) -> GhashPowerMap {
78	let matrix = compute_power_map_matrix(n);
79	OutputWrappingTransformationFactory::<_, Ghash128b, Ghash128b>::new(
80		InputWrappingTransformationFactory::<_, Ghash128b, M128>::new(
81			BytewiseLookupTransformationFactory,
82		),
83	)
84	.create(&matrix)
85}
86
87/// Compute the matrix of the `F_2`-linear power map `x -> x^(2^n)`.
88///
89/// Column `i` is the image of the `i`-th basis element, i.e. `basis(i)^(2^n)`, obtained by squaring
90/// `n` times.
91fn compute_power_map_matrix(n: usize) -> [Ghash128b; FIELD_BITS] {
92	array::from_fn(|i| {
93		let basis = <Ghash128b as ExtensionField<BinaryField1b>>::basis(i);
94		iter::successors(Some(basis), |basis_pow_2_i| Some(basis_pow_2_i.square()))
95			.nth(n)
96			.expect("closure always returns Some")
97	})
98}
99
100/// Invert each GHASH element (scalar or packed) via the Itoh-Tsujii algorithm.
101///
102/// Zero elements map to zero, matching `InvertOrZero` semantics.
103///
104/// The bound is phrased in terms of the field operations (`Square`, `Mul`) plus
105/// `Divisible<Ghash128b>` rather than `P: PackedField`. `PackedField`'s blanket impl lists
106/// `InvertOrZero` in its where-clause, so requiring it here would form a trait-resolution cycle
107/// when this function backs the `InvertOrZero` impls. `Divisible<Ghash128b>` carries no such
108/// obligation, keeps the function statically GHASH-typed, and is satisfied both by the GHASH packed
109/// fields and (reflexively) by the scalar `Ghash128b`, so the scalar inverts directly
110/// without routing through a packed type.
111pub fn invert_b128<P>(x: P) -> P
112where
113	P: Copy + Square + Mul<Output = P> + Divisible<Ghash128b>,
114{
115	let tables = &*GHASH_POWER_MAP_TABLES;
116
117	// Addition chain for 127: 1, 2, 3, 6, 7, 14, 28, 56, 63, 126, 127.
118	let beta_1 = x;
119	let beta_2 = beta_1.square() * beta_1;
120	let beta_3 = beta_2.square() * beta_1;
121	let beta_6 = pow_2_n(beta_3, &tables.pow_2_3) * beta_3;
122	let beta_7 = beta_6.square() * beta_1;
123	let beta_14 = pow_2_n(beta_7, &tables.pow_2_7) * beta_7;
124	let beta_28 = pow_2_n(beta_14, &tables.pow_2_14) * beta_14;
125	let beta_56 = pow_2_n(beta_28, &tables.pow_2_28) * beta_28;
126	let beta_63 = pow_2_n(beta_56, &tables.pow_2_7) * beta_7;
127	let beta_126 = pow_2_n(beta_63, &tables.pow_2_63) * beta_63;
128	let beta_127 = beta_126.square() * beta_1;
129	// x^(-1) = (x^(2^127 - 1))^2.
130	beta_127.square()
131}
132
133/// Apply the power map `x -> x^(2^n)` to every GHASH scalar of `x`.
134fn pow_2_n<P>(x: P, power_map: &GhashPowerMap) -> P
135where
136	P: Divisible<Ghash128b>,
137{
138	Divisible::<Ghash128b>::from_iter(
139		Divisible::<Ghash128b>::value_iter(x).map(|scalar| power_map.transform(&scalar)),
140	)
141}
142
143/// `InvertOrZero` strategy wrapper backed by the [Itoh-Tsujii](invert_b128) inversion.
144///
145/// This is the single inversion strategy for the GHASH field across every architecture — there is
146/// no carryless-multiply inverse, so the same addition-chain algorithm applies whether the square
147/// and multiply underneath are CLMUL/PMULL-accelerated or software. Each arch type-aliases its
148/// `GhashInvert1x` to this wrapper.
149#[repr(transparent)]
150#[derive(TransparentWrapper)]
151pub struct GhashItohTsujii<T>(T);
152
153impl<P> InvertOrZero for GhashItohTsujii<P>
154where
155	P: Copy + Square + Mul<Output = P> + Divisible<Ghash128b>,
156{
157	#[inline]
158	fn invert_or_zero(self) -> Self {
159		Self::wrap(invert_b128(Self::peel(self)))
160	}
161}
162
163#[cfg(test)]
164mod tests {
165	use proptest::prelude::*;
166
167	use super::*;
168	use crate::{Field, PackedField, PackedGhash1x128b, PackedGhash2x128b};
169
170	#[test]
171	fn test_compute_power_map_matrix_is_squaring() {
172		// The 2^1 power map is just squaring; column i must equal basis(i)^2.
173		let matrix = compute_power_map_matrix(1);
174		for i in 0..FIELD_BITS {
175			let basis = <Ghash128b as ExtensionField<BinaryField1b>>::basis(i);
176			assert_eq!(matrix[i], basis.square());
177		}
178	}
179
180	#[test]
181	fn test_power_map_transform_matches_repeated_squaring() {
182		let power_map = compute_power_map_transform(7);
183		for &raw in &[0u128, 1, 2, 0x87, 0x21ac73a21d46a21badd6747bcdfc5d4d] {
184			let x = Ghash128b::from(raw);
185			let mut expected = x;
186			for _ in 0..7 {
187				expected = expected.square();
188			}
189			assert_eq!(power_map.transform(&x), expected);
190		}
191	}
192
193	#[test]
194	fn test_invert_b128_known_values() {
195		let one = PackedGhash1x128b::broadcast(Ghash128b::ONE);
196		assert_eq!(invert_b128(one), one);
197
198		let zero = PackedGhash1x128b::broadcast(Ghash128b::ZERO);
199		assert_eq!(invert_b128(zero), zero);
200	}
201
202	// `invert_b128` now backs `InvertOrZero` itself, so the multiplicative-inverse property (with
203	// `0 -> 0`) is the independent oracle: given a separately-tested `mul`, `x * x^-1 == 1` fully
204	// characterizes invert-or-zero.
205	proptest! {
206		#[test]
207		fn test_invert_b128_is_multiplicative_inverse_scalar(raw in any::<u128>()) {
208			let x = Ghash128b::from(raw);
209			let inv = invert_b128(x);
210			if x == Ghash128b::ZERO {
211				prop_assert_eq!(inv, Ghash128b::ZERO);
212			} else {
213				prop_assert_eq!(x * inv, Ghash128b::ONE);
214			}
215		}
216
217		#[test]
218		fn test_invert_b128_is_multiplicative_inverse_1x(raw in any::<u128>()) {
219			let scalar = Ghash128b::from(raw);
220			let x = PackedGhash1x128b::broadcast(scalar);
221			let inv = invert_b128(x);
222			if scalar == Ghash128b::ZERO {
223				prop_assert_eq!(inv, x);
224			} else {
225				prop_assert_eq!(x * inv, PackedGhash1x128b::broadcast(Ghash128b::ONE));
226			}
227		}
228
229		#[test]
230		fn test_invert_b128_is_multiplicative_inverse_2x(a in any::<u128>(), b in any::<u128>()) {
231			let x = PackedGhash2x128b::from_scalars([a, b].map(Ghash128b::from));
232			let inv = invert_b128(x);
233			let ones = PackedGhash2x128b::from_scalars(
234				[a, b].map(|raw| {
235					if Ghash128b::from(raw) == Ghash128b::ZERO {
236						Ghash128b::ZERO
237					} else {
238						Ghash128b::ONE
239					}
240				}),
241			);
242			prop_assert_eq!(x * inv, ones);
243		}
244	}
245}
binius_field/arch/portable/arithmetic/itoh_tsujii.rs

binius_field/arch/portable/arithmetic/
itoh_tsujii.rs
