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

binius_field/arch/portable/arithmetic/
itoh_tsujii.rs
