Initial commit

.gitignore

@@ -0,0 +1,2 @@
+/target
+*.sage.py

Cargo.lock

@@ -0,0 +1,25 @@
+# This file is automatically @generated by Cargo.
+# It is not intended for manual editing.
+version = 4
+
+[[package]]
+name = "autocfg"
+version = "1.5.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c08606f8c3cbf4ce6ec8e28fb0014a2c086708fe954eaa885384a6165172e7e8"
+
+[[package]]
+name = "gwrizienn"
+version = "0.1.0"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "num-traits"
+version = "0.2.19"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841"
+dependencies = [
+ "autocfg",
+]

Cargo.toml

@@ -0,0 +1,7 @@
+[package]
+name = "gwrizienn"
+version = "0.1.0"
+edition = "2024"
+
+[dependencies]
+num-traits = "0.2.19"

rustfmt.toml

@@ -0,0 +1,9 @@
+hard_tabs = true
+newline_style = "Unix"
+imports_granularity = "Crate"
+
+unstable_features = true
+format_code_in_doc_comments = true
+format_macro_bodies = true
+format_macro_matchers = true
+format_strings = true

src/lib.rs

@@ -0,0 +1,7 @@
+pub mod matrix;
+pub mod ntt;
+pub mod poly;
+//pub mod r7681;
+pub mod ring;
+pub mod tuple;
+pub mod vector;

src/matrix.rs

@@ -0,0 +1 @@
+

src/ntt.rs

@@ -0,0 +1,103 @@
+use num_traits::Zero;
+use std::ops::{Add, AddAssign, Sub, SubAssign};
+
+#[derive(Debug, Clone, Eq, PartialEq)]
+pub struct NttDomain<E>(pub E);
+
+pub trait Ntt: Sized {
+	type Output: Sized;
+	fn ntt(self) -> Self::Output;
+}
+
+pub trait NttInv: Sized {
+	type Output: Sized;
+	fn ntt_inv(self) -> Self::Output;
+}
+
+impl<E> Add<NttDomain<E>> for NttDomain<E>
+where
+	E: Add<E, Output = E>,
+{
+	type Output = NttDomain<E>;
+	fn add(self, rhs: NttDomain<E>) -> Self::Output {
+		NttDomain(self.0 + rhs.0)
+	}
+}
+
+impl<'a, E> Add<&'a NttDomain<E>> for NttDomain<E>
+where
+	E: 'a + Add<&'a E, Output = E>,
+{
+	type Output = NttDomain<E>;
+	fn add(self, rhs: &'a NttDomain<E>) -> Self::Output {
+		NttDomain(self.0 + &rhs.0)
+	}
+}
+
+impl<E> Sub<NttDomain<E>> for NttDomain<E>
+where
+	E: Sub<E, Output = E>,
+{
+	type Output = NttDomain<E>;
+	fn sub(self, rhs: NttDomain<E>) -> Self::Output {
+		NttDomain(self.0 - rhs.0)
+	}
+}
+
+impl<'a, E> Sub<&'a NttDomain<E>> for NttDomain<E>
+where
+	E: 'a + Sub<&'a E, Output = E>,
+{
+	type Output = NttDomain<E>;
+	fn sub(self, rhs: &'a NttDomain<E>) -> Self::Output {
+		NttDomain(self.0 - &rhs.0)
+	}
+}
+
+impl<E> Zero for NttDomain<E>
+where
+	E: Zero,
+{
+	fn zero() -> Self {
+		NttDomain(E::zero())
+	}
+	fn is_zero(&self) -> bool {
+		self.0.is_zero()
+	}
+}
+
+impl<E> AddAssign for NttDomain<E>
+where
+	E: AddAssign,
+{
+	fn add_assign(&mut self, rhs: Self) {
+		self.0 += rhs.0;
+	}
+}
+
+impl<'a, E> AddAssign<&'a Self> for NttDomain<E>
+where
+	E: 'a + AddAssign<&'a E>,
+{
+	fn add_assign(&mut self, rhs: &'a Self) {
+		self.0 += &rhs.0;
+	}
+}
+
+impl<E> SubAssign for NttDomain<E>
+where
+	E: SubAssign,
+{
+	fn sub_assign(&mut self, rhs: Self) {
+		self.0 -= rhs.0;
+	}
+}
+
+impl<'a, E> SubAssign<&'a Self> for NttDomain<E>
+where
+	E: 'a + SubAssign<&'a E>,
+{
+	fn sub_assign(&mut self, rhs: &'a Self) {
+		self.0 -= &rhs.0;
+	}
+}

src/poly.rs

@@ -0,0 +1,357 @@
+use crate::{ntt::*, ring::Ring};
+
+use num_traits::Zero;
+
+pub trait Poly<T, const N: usize> {
+	type Element: Ring;
+
+	fn modulus() -> T;
+	fn new(v: [Self::Element; N]) -> Self;
+}
+
+pub(crate) const fn bitrev(x: usize, n: usize) -> usize {
+	x.reverse_bits() >> (usize::BITS - n.ilog2())
+}
+
+#[macro_export]
+macro_rules! impl_poly_ring {
+	($Rq:ident, $N:literal, $Zq:ident, $T:ident, $TT:ident, $q:expr) => {
+		/// Element of (Z/qZ)[x]/(x^N+1)
+		#[derive(Clone, Copy, Debug, Eq, Hash, PartialEq)]
+		pub struct $Rq(pub [$Zq; $N]);
+
+		impl $crate::poly::Poly<$T, $N> for $Rq {
+			type Element = $Zq;
+			fn modulus() -> $T {
+				$q
+			}
+			fn new(v: [$Zq; $N]) -> Self {
+				Self(v)
+			}
+		}
+
+		impl std::ops::Add for $Rq {
+			type Output = Self;
+			fn add(mut self, rhs: Self) -> Self {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x += *y);
+				self
+			}
+		}
+
+		impl std::ops::Add<&Self> for $Rq {
+			type Output = Self;
+			fn add(mut self, rhs: &Self) -> Self {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x += *y);
+				self
+			}
+		}
+
+		impl std::ops::AddAssign for $Rq {
+			fn add_assign(&mut self, rhs: Self) {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x += *y);
+			}
+		}
+
+		impl std::ops::AddAssign<&Self> for $Rq {
+			fn add_assign(&mut self, rhs: &Self) {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x += *y);
+			}
+		}
+
+		impl std::ops::Sub for $Rq {
+			type Output = Self;
+			fn sub(mut self, rhs: Self) -> Self {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x -= *y);
+				self
+			}
+		}
+
+		impl std::ops::Sub<&Self> for $Rq {
+			type Output = Self;
+			fn sub(mut self, rhs: &Self) -> Self {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x -= *y);
+				self
+			}
+		}
+
+		impl std::ops::SubAssign for $Rq {
+			fn sub_assign(&mut self, rhs: Self) {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x -= *y);
+			}
+		}
+
+		impl std::ops::SubAssign<&Self> for $Rq {
+			fn sub_assign(&mut self, rhs: &Self) {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x -= *y);
+			}
+		}
+	};
+}
+
+#[macro_export]
+macro_rules! poly {
+	($Rq:ident, $N:literal, $Zq:ident, $T:ident, $TT:ident, $q:expr, $zeta:expr) => {
+		$crate::impl_poly_ring!($Rq, $N, $Zq, $T, $TT, $q);
+		impl $Rq {
+			pub const fn zetas() -> [$T; $N] {
+				let mut zetas = [0; $N];
+				let mut i: usize = 1;
+				let mut z: $TT = $zeta;
+				while i < $N {
+					zetas[$crate::poly::bitrev(i, $N)] = z as $T;
+					z = z * $zeta % $q;
+					i += 1;
+				}
+				zetas
+			}
+			pub const fn one() -> Self {
+				let mut x = Self([$Zq(0); $N]);
+				x.0[0].0 = 1;
+				x
+			}
+			pub fn is_one(&self) -> bool {
+				self.0[0].0 == 1 && self.0.iter().skip(1).all(|x| x.0 == 0)
+			}
+			/*/// Slow multiplication, for testing
+			fn naive_mul(&self, rhs: &Self) -> Self {
+				let mut buf = [$Zq(0); 2*$N];
+				for (i, a) in self.0.iter().enumerate() {
+					for (j, b) in self.0.iter().enumerate() {
+						buf[i+j] = a*b;
+					}
+				}
+				buf
+			}*/
+		}
+
+		impl num_traits::Zero for $Rq {
+			fn zero() -> Self {
+				Self([$Zq(0); $N])
+			}
+			fn is_zero(&self) -> bool {
+				self.0.iter().all(num_traits::Zero::is_zero)
+			}
+		}
+
+		/*impl num_traits::One for $Rq {
+			fn one() -> Self {
+				let mut x = Self([$Zq(0); $N]);
+				x.0[0].0 = 1;
+				x
+			}
+			fn is_one(&self) -> bool {
+				self.0[0].0 == 1 && self.0.iter().skip(1).all(|x| x.0 == 0)
+			}
+		}*/
+
+		impl $crate::ntt::Ntt for $Rq {
+			type Output = $crate::ntt::NttDomain<$Rq>;
+			fn ntt(mut self) -> $crate::ntt::NttDomain<$Rq> {
+				let mut m = 0;
+				let mut len = $N / 2;
+				while len >= 1 {
+					for start in (0..$N).step_by(2 * len) {
+						m += 1;
+						let z = $Zq($Rq::zetas()[m]);
+						for j in start..start + len {
+							let t = z * self.0[j + len];
+							self.0[j + len] = self.0[j] - t;
+							self.0[j] += t;
+						}
+					}
+					len /= 2;
+				}
+				$crate::ntt::NttDomain(self)
+			}
+		}
+
+		impl $crate::ntt::NttInv for $crate::ntt::NttDomain<$Rq> {
+			type Output = $Rq;
+			fn ntt_inv(mut self) -> $Rq {
+				let mut m = $N;
+				let mut len = 1;
+				while len < $N {
+					for start in (0..$N).step_by(2 * len) {
+						m -= 1;
+						let z = $Zq($q - $Rq::zetas()[m]);
+						for j in start..start + len {
+							let t = self.0.0[j];
+							self.0.0[j] += self.0.0[j + len];
+							self.0.0[j + len] = z * (t - self.0.0[j + len]);
+						}
+					}
+					len *= 2;
+				}
+				self.0.0.iter_mut().for_each(|wj| *wj *= $Zq($N).inverse());
+				self.0
+			}
+		}
+
+		impl std::ops::Mul for $crate::ntt::NttDomain<$Rq> {
+			type Output = Self;
+			fn mul(mut self, rhs: Self) -> Self {
+				self.0
+					.0
+					.iter_mut()
+					.zip(rhs.0.0.iter())
+					.for_each(|(x, y)| *x *= *y);
+				self
+			}
+		}
+
+		impl std::ops::Mul<&Self> for $crate::ntt::NttDomain<$Rq> {
+			type Output = Self;
+			fn mul(mut self, rhs: &Self) -> Self {
+				self.0
+					.0
+					.iter_mut()
+					.zip(rhs.0.0.iter())
+					.for_each(|(x, y)| *x *= *y);
+				self
+			}
+		}
+
+		impl std::ops::MulAssign for $crate::ntt::NttDomain<$Rq> {
+			fn mul_assign(&mut self, rhs: Self) {
+				self.0
+					.0
+					.iter_mut()
+					.zip(rhs.0.0.iter())
+					.for_each(|(x, y)| *x *= *y);
+			}
+		}
+
+		impl std::ops::MulAssign<&Self> for $crate::ntt::NttDomain<$Rq> {
+			fn mul_assign(&mut self, rhs: &Self) {
+				self.0
+					.0
+					.iter_mut()
+					.zip(rhs.0.0.iter())
+					.for_each(|(x, y)| *x *= *y);
+			}
+		}
+
+		impl num_traits::Inv for $crate::ntt::NttDomain<$Rq> {
+			type Output = Self;
+			fn inv(mut self) -> Self {
+				self.0.0.iter_mut().for_each(|x| *x = x.inverse());
+				self
+			}
+		}
+	};
+}
+
+#[macro_export]
+macro_rules! poly2 {
+	($Rq:ident, $Rr:ident, $R2q:ident, $N:literal, $Zq:ident, $Z2q:ident, $T:ident, $TT:ident, $TI:ident, $q:expr, $r:expr) => {
+		$crate::impl_poly_ring!($R2q, $N, $Z2q, $T, $TT, 2*$q);
+		impl $R2q {
+			pub const fn new(v: [$Z2q; $N]) -> Self {
+				Self(v)
+			}
+		}
+		impl $crate::ntt::Ntt for $R2q {
+			type Output = $crate::tuple::Tuple<$crate::ntt::NttDomain<$Rr>, $crate::ntt::NttDomain<$Rq>>;
+			fn ntt(self) -> Self::Output {
+				let mut a = $Rr::zero();
+				let mut b = $Rq::zero();
+				for ((x_i, a_i), b_i) in self.0.iter().zip(a.0.iter_mut()).zip(b.0.iter_mut()) {
+					a_i.0 = (x_i.0 % 2) as _;
+					b_i.0 = (x_i.0 % $q) as _;
+				}
+				$crate::tuple::Tuple(a.ntt(), b.ntt())
+			}
+		}
+		impl $crate::ntt::NttInv for $crate::tuple::Tuple<$crate::ntt::NttDomain<$Rr>, $crate::ntt::NttDomain<$Rq>> {
+			type Output = $R2q;
+			fn ntt_inv(self) -> Self::Output {
+				let mut x = $R2q::new([$Z2q::zero(); $N]);
+				for ((x_i, a_i), b_i) in x.0.iter_mut().zip(self.0.ntt_inv().0.iter()).zip(self.1.ntt_inv().0.iter()) {
+					let mut tmp = ((a_i.0 % 2) as $TI *$q - ($q/2*2*b_i.0) as $TI) % (2*$q);
+					if tmp < 0 {
+						tmp += 2*$q;
+					}
+					x_i.0 = tmp as $T;
+				}
+				x
+			}
+		}
+	};
+}
+
+#[cfg(test)]
+mod test {
+	use super::*;
+
+	use crate::ring;
+
+	use num_traits::Inv;
+
+	crate::ring!(Z50177, u32, u64, i64, 50177);
+	poly!(R50177, 256, Z50177, u32, u64, 50177, 66);
+	crate::ring!(Z7681, u16, u32, i32, 7681);
+	poly!(R7681, 256, Z7681, u16, u32, 7681, 62);
+	poly2!(R50177, R7681, R100354, 256, Z50177, Z100354, u32, u64, i64, 50177, 7681);
+	crate::ring!(Z100354, u32, u64, i64, 100354);
+
+	#[test]
+	fn test_bitrev() {
+		assert_eq!(bitrev(0b10101100, 256), 0b00110101);
+		assert_eq!(bitrev(0b01010101, 256), 0b10101010);
+	}
+
+	#[test]
+	fn test_basic() {
+		assert_eq!(R50177::modulus(), 50177);
+	}
+
+	#[test]
+	fn test_ntt() {
+		let a = R50177::new([Z50177::new(42); 256]);
+		let b = R50177::new([Z50177::new(123); 256]);
+		assert_eq!(a, a.clone().ntt().ntt_inv());
+		assert_eq!(a + b, (a.clone().ntt() + b.clone().ntt()).ntt_inv());
+	}
+
+	#[test]
+	fn test_ntt_inverse() {
+		let a = R50177::new([Z50177::new(42); 256]).ntt();
+		assert!((a.clone().inv() * a).ntt_inv().is_one());
+	}
+
+	#[test]
+	fn test_ntt_inverse_composite() {
+		let a = R100354::new([Z100354::new(41); 256]);
+		let mut b = R100354::new([Z100354::new(0); 256]);
+		b.0[0].0 = 2;
+		let c = R100354::new([Z100354::new(82); 256]);
+		assert_eq!(a.clone().ntt().ntt_inv(), a);
+		assert_eq!(b.clone().ntt().ntt_inv(), b);
+		assert_eq!(c.clone().ntt().ntt_inv(), c);
+		assert_eq!((a.ntt() * b.ntt()).ntt_inv(), c);
+	}
+}

src/r7681.rs

@@ -0,0 +1,232 @@
+crate::ring!(Z7681, u16, u32, i32, 7681);
+
+/// Element of (Z/2Z)[x]/(x^N+1)
+		#[derive(Clone, Copy, Debug, Eq, Hash, PartialEq)]
+		pub struct R7681<const N: usize>(pub [Z7681; N]);
+
+		impl<const N: usize> crate::poly::Poly<u16, N> for R7681<N> {
+			type Element = Z7681;
+			fn modulus() -> u16 {
+				7681
+			}
+			fn new(v: [Z7681; N]) -> Self {
+				Self(v)
+			}
+		}
+
+		impl<const N: usize> R7681<N> {
+			pub const fn one() -> Self {
+				let mut x = Self([Z7681(0); N]);
+				x.0[0].0 = 1;
+				x
+			}
+			pub fn is_one(&self) -> bool {
+				self.0[0].0 == 1 && self.0.iter().skip(1).all(|x| x.0 == 0)
+			}
+			pub const fn zetas() -> [u16; N] {
+				let mut zetas = [0; N];
+				let mut i: usize = 1;
+				let mut z = 62;
+				while i < N {
+					zetas[crate::poly::bitrev(i, N)] = z;
+					z = (z as u32 * 62 % 7681) as u16;
+					i += 1;
+				}
+				zetas
+			}
+		}
+
+		impl<const N: usize> std::ops::Add for R7681<N> {
+			type Output = Self;
+			fn add(mut self, rhs: Self) -> Self {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x += *y);
+				self
+			}
+		}
+
+		impl<const N: usize> std::ops::Add<&Self> for R7681<N> {
+			type Output = Self;
+			fn add(mut self, rhs: &Self) -> Self {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x += *y);
+				self
+			}
+		}
+
+		impl<const N: usize> std::ops::AddAssign for R7681<N> {
+			fn add_assign(&mut self, rhs: Self) {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x += *y);
+			}
+		}
+
+		impl<const N: usize> std::ops::AddAssign<&Self> for R7681<N> {
+			fn add_assign(&mut self, rhs: &Self) {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x += *y);
+			}
+		}
+
+		impl<const N: usize> std::ops::Sub for R7681<N> {
+			type Output = Self;
+			fn sub(mut self, rhs: Self) -> Self {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x -= *y);
+				self
+			}
+		}
+
+		impl<const N: usize> std::ops::Sub<&Self> for R7681<N> {
+			type Output = Self;
+			fn sub(mut self, rhs: &Self) -> Self {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x -= *y);
+				self
+			}
+		}
+
+		impl<const N: usize> std::ops::SubAssign for R7681<N> {
+			fn sub_assign(&mut self, rhs: Self) {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x -= *y);
+			}
+		}
+
+		impl<const N: usize> std::ops::SubAssign<&Self> for R7681<N> {
+			fn sub_assign(&mut self, rhs: &Self) {
+				self.0
+					.iter_mut()
+					.zip(rhs.0.iter())
+					.for_each(|(x, y)| *x -= *y);
+			}
+		}
+
+		impl<const N: usize> num_traits::Zero for R7681<N> {
+			fn zero() -> Self {
+				Self([Z7681(0); N])
+			}
+			fn is_zero(&self) -> bool {
+				self.0.iter().all(num_traits::Zero::is_zero)
+			}
+		}
+
+		/*impl num_traits::One for R7681<N> {
+			fn one() -> Self {
+				let mut x = Self([$Zq(0); $N]);
+				x.0[0].0 = 1;
+				x
+			}
+			fn is_one(&self) -> bool {
+				self.0[0].0 == 1 && self.0.iter().skip(1).all(|x| x.0 == 0)
+			}
+		}*/
+
+		impl<const N: usize> crate::ntt::Ntt for R7681<N> {
+			type Output = crate::ntt::NttDomain<R7681<N>>;
+			fn ntt(mut self) -> crate::ntt::NttDomain<R7681<N>> {
+				let mut m = 0;
+				let mut len = N / 2;
+				while len >= 1 {
+					for start in (0..N).step_by(2 * len) {
+						m += 1;
+						let z = Z7681(R7681::<N>::zetas()[m]);
+						for j in start..start + len {
+							let t = z * self.0[j + len];
+							self.0[j + len] = self.0[j] - t;
+							self.0[j] += t;
+						}
+					}
+					len /= 2;
+				}
+				crate::ntt::NttDomain(self)
+			}
+		}
+
+		impl<const N: usize> crate::ntt::NttInv for crate::ntt::NttDomain<R7681<N>> {
+			type Output = R7681<N>;
+			fn ntt_inv(mut self) -> R7681<N> {
+				let mut m = N;
+				let mut len = 1;
+				while len < N {
+					for start in (0..N).step_by(2 * len) {
+						m -= 1;
+						let z = Z7681(7681 - R7681::<N>::zetas()[m]);
+						for j in start..start + len {
+							let t = self.0.0[j];
+							self.0.0[j] += self.0.0[j + len];
+							self.0.0[j + len] = z * (t - self.0.0[j + len]);
+						}
+					}
+					len *= 2;
+				}
+				self.0.0.iter_mut().for_each(|wj| *wj *= Z7681(N as u16).inverse());
+				self.0
+			}
+		}
+
+		impl<const N: usize> std::ops::Mul for crate::ntt::NttDomain<R7681<N>> {
+			type Output = Self;
+			fn mul(mut self, rhs: Self) -> Self {
+				self.0
+					.0
+					.iter_mut()
+					.zip(rhs.0.0.iter())
+					.for_each(|(x, y)| *x *= *y);
+				self
+			}
+		}
+
+		impl<const N: usize> std::ops::Mul<&Self> for crate::ntt::NttDomain<R7681<N>> {
+			type Output = Self;
+			fn mul(mut self, rhs: &Self) -> Self {
+				self.0
+					.0
+					.iter_mut()
+					.zip(rhs.0.0.iter())
+					.for_each(|(x, y)| *x *= *y);
+				self
+			}
+		}
+
+		impl<const N: usize> std::ops::MulAssign for crate::ntt::NttDomain<R7681<N>> {
+			fn mul_assign(&mut self, rhs: Self) {
+				self.0
+					.0
+					.iter_mut()
+					.zip(rhs.0.0.iter())
+					.for_each(|(x, y)| *x *= *y);
+			}
+		}
+
+		impl<const N: usize> std::ops::MulAssign<&Self> for crate::ntt::NttDomain<R7681<N>> {
+			fn mul_assign(&mut self, rhs: &Self) {
+				self.0
+					.0
+					.iter_mut()
+					.zip(rhs.0.0.iter())
+					.for_each(|(x, y)| *x *= *y);
+			}
+		}
+
+		impl<const N: usize> num_traits::Inv for crate::ntt::NttDomain<R7681<N>> {
+			type Output = Self;
+			fn inv(mut self) -> Self {
+				self.0.0.iter_mut().for_each(|x| *x = x.inverse());
+				self
+			}
+		}

src/ring.rs

@@ -0,0 +1,133 @@
+pub trait Ring {
+	type T;
+}
+
+#[macro_export]
+macro_rules! ring {
+	($Zq:ident, $T:ident, $TT:ident, $TI:ident, $q:expr) => {
+		/// Element of Z/qZ
+		#[derive(Clone, Copy, Debug, Eq, Hash, PartialEq)]
+		pub struct $Zq(pub $T);
+
+		impl $crate::ring::Ring for $Zq {
+			type T = $T;
+		}
+
+		impl $Zq {
+			pub const fn modulus() -> $T {
+				$q
+			}
+			pub const fn new(v: $T) -> Self {
+				Self(v % $q)
+			}
+			const fn small_reduce(v: $T) -> $T {
+				if v < $q { v } else { v - $q }
+			}
+			const fn inverse(self) -> Self {
+				let (mut ro, mut r) = ($q, self.0);
+				let (mut to, mut t) = (0, 1);
+				while r != 0 {
+					let quo = ro / r;
+					(ro, r) = (r, ro - quo * r);
+					(to, t) = (t, to - quo as $TI * t);
+				}
+				if to < 0 {
+					to += $q;
+				}
+				Self(to as $T)
+			}
+		}
+
+		impl std::ops::Add for $Zq {
+			type Output = Self;
+			fn add(self, rhs: Self) -> Self {
+				Self(Self::small_reduce(self.0 + rhs.0))
+			}
+		}
+
+		impl std::ops::AddAssign for $Zq {
+			fn add_assign(&mut self, rhs: Self) {
+				*self = *self + rhs;
+			}
+		}
+
+		impl std::ops::Sub for $Zq {
+			type Output = Self;
+			fn sub(self, rhs: Self) -> Self {
+				Self(Self::small_reduce($q + self.0 - rhs.0))
+			}
+		}
+
+		impl std::ops::SubAssign for $Zq {
+			fn sub_assign(&mut self, rhs: Self) {
+				*self = *self - rhs;
+			}
+		}
+
+		impl std::ops::Neg for $Zq {
+			type Output = Self;
+			fn neg(self) -> Self {
+				Self(Self::small_reduce($q - self.0))
+			}
+		}
+
+		impl num_traits::Zero for $Zq {
+			fn zero() -> Self {
+				Self(0)
+			}
+			fn is_zero(&self) -> bool {
+				self.0 == 0
+			}
+		}
+
+		impl std::ops::Mul for $Zq {
+			type Output = Self;
+			fn mul(self, rhs: Self) -> Self {
+				let product = self.0 as $TT * rhs.0 as $TT;
+				Self((product % $q as $TT) as $T)
+			}
+		}
+
+		impl std::ops::MulAssign for $Zq {
+			fn mul_assign(&mut self, rhs: Self) {
+				*self = *self * rhs;
+			}
+		}
+
+		impl num_traits::One for $Zq {
+			fn one() -> Self {
+				Self(1)
+			}
+			fn is_one(&self) -> bool {
+				self.0 == 1
+			}
+		}
+
+		impl num_traits::Inv for $Zq {
+			type Output = Self;
+			fn inv(self) -> Self {
+				self.inverse()
+			}
+		}
+	};
+}
+
+crate::ring!(Z2, u8, u8, i8, 2);
+
+#[cfg(test)]
+mod test {
+	use super::*;
+
+	use num_traits::*;
+
+	#[test]
+	fn test_basic() {
+		ring!(Z50177, u32, u64, i64, 50177);
+
+		assert_eq!(Z50177::modulus(), 50177);
+
+		assert_eq!(Z50177::new(1), Z50177::new(50178));
+
+		assert_eq!(Z50177::new(123) * Z50177::new(123).inv(), Z50177::one());
+	}
+}

src/tuple.rs

@@ -0,0 +1,26 @@
+use std::ops::{Add, Mul};
+use num_traits::Inv;
+
+#[derive(Clone, Debug, Eq, Hash, PartialEq)]
+pub struct Tuple<A, B>(pub A, pub B);
+
+impl<A, B> Add for Tuple<A, B> where A: Add<A, Output=A>, B: Add<B, Output=B> {
+	type Output = Self;
+	fn add(self, rhs: Self) -> Self {
+		Self(self.0+rhs.0, self.1+rhs.1)
+	}
+}
+
+impl<A, B> Mul for Tuple<A, B> where A: Mul<A, Output=A>, B: Mul<B, Output=B> {
+	type Output = Self;
+	fn mul(self, rhs: Self) -> Self {
+		Self(self.0*rhs.0, self.1*rhs.1)
+	}
+}
+
+impl<A,B> Inv for Tuple<A, B> where A: Inv<Output=A>, B: Inv<Output=B> {
+	type Output = Self;
+	fn inv(self) -> Self {
+		Self(self.0.inv(), self.1.inv())
+	}
+}

src/vector.rs

@@ -0,0 +1,188 @@
+use std::{
+	mem::MaybeUninit,
+	ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign},
+};
+
+use num_traits::{One, Zero};
+
+#[derive(Clone, Debug, Eq, Hash, PartialEq)]
+pub struct Vector<E, const N: usize>(pub [E; N]);
+
+/*impl<const N: usize> Vector<Poly, N> {
+	pub fn norm_inf(&self) -> Zq {
+		Zq(self.0.iter().map(|i| i.norm_inf().0).max().unwrap())
+	}
+}*/
+
+impl<E, const N: usize> Zero for Vector<E, N>
+where
+	E: AddAssign + Zero,
+{
+	fn zero() -> Self {
+		let mut z = MaybeUninit::<[E; N]>::uninit();
+		unsafe {
+			for i in z.assume_init_mut() {
+				*i = E::zero();
+			}
+		}
+		Self(unsafe { z.assume_init() })
+	}
+	fn is_zero(&self) -> bool {
+		self.0.iter().all(Zero::is_zero)
+	}
+}
+
+impl<E, const N: usize> Add<Self> for Vector<E, N>
+where
+	E: AddAssign<E>,
+{
+	type Output = Self;
+	fn add(mut self, rhs: Self) -> Self {
+		self.0.iter_mut().zip(rhs.0).for_each(|(x, y)| *x += y);
+		self
+	}
+}
+
+impl<'a, E, const N: usize> Add<&'a Self> for Vector<E, N>
+where
+	E: 'a + AddAssign<&'a E>,
+{
+	type Output = Self;
+	fn add(mut self, rhs: &'a Self) -> Self {
+		self.0
+			.iter_mut()
+			.zip(rhs.0.iter())
+			.for_each(|(x, y)| *x += y);
+		self
+	}
+}
+
+impl<E, const N: usize> Sub<Self> for Vector<E, N>
+where
+	E: SubAssign<E>,
+{
+	type Output = Self;
+	fn sub(mut self, rhs: Self) -> Self {
+		self.0.iter_mut().zip(rhs.0).for_each(|(x, y)| *x -= y);
+		self
+	}
+}
+
+impl<'a, E, const N: usize> Sub<&'a Self> for Vector<E, N>
+where
+	E: 'a + SubAssign<&'a E>,
+{
+	type Output = Self;
+	fn sub(mut self, rhs: &'a Self) -> Self {
+		self.0
+			.iter_mut()
+			.zip(rhs.0.iter())
+			.for_each(|(x, y)| *x -= y);
+		self
+	}
+}
+
+/// Vector-scalar product
+impl<'a, E, const N: usize> Mul<&'a E> for Vector<E, N>
+where
+	E: 'a + Clone + std::fmt::Debug + Zero + MulAssign<&'a E>,
+{
+	type Output = Vector<E, N>;
+	fn mul(mut self, rhs: &'a E) -> Vector<E, N> {
+		self.0.iter_mut().for_each(|s| *s *= rhs);
+		self
+	}
+}
+
+/// Inner product
+// TODO optimize (reduce and copy less often)
+impl<'a, E, const N: usize> Mul<&'a Self> for Vector<E, N>
+where
+	E: 'a + Clone + Zero + MulAssign<&'a E> + Add<E>,
+{
+	type Output = E;
+	fn mul(self, rhs: &'a Self) -> E {
+		self.0
+			.into_iter()
+			.zip(rhs.0.iter())
+			.fold(E::zero(), |e, (mut x, y)| {
+				x *= y;
+				e + x
+			})
+	}
+}
+
+/// Inner product
+// TODO optimize (reduce and copy less often)
+impl<'a, E, const N: usize> Mul<Self> for &'a Vector<E, N>
+where
+	E: 'a + Clone + Zero + Mul<&'a E, Output = E> + Add<E>,
+{
+	type Output = E;
+	fn mul(self, rhs: Self) -> E {
+		self.0
+			.iter()
+			.zip(rhs.0.iter())
+			.fold(E::zero(), |e, (x, y)| e + x.clone() * y)
+	}
+}
+
+/// Matrix-vector product
+impl<'a, E, const M: usize, const N: usize> Mul<&'a Vector<E, M>> for Vector<Vector<E, M>, N>
+where
+	E: 'a + Clone + std::fmt::Debug + Zero + MulAssign<&'a E> + Add<E> + AddAssign<E>,
+{
+	type Output = Vector<E, N>;
+	fn mul(self, rhs: &'a Vector<E, M>) -> Vector<E, N> {
+		let mut res = Self::Output::zero();
+		res.0.iter_mut().zip(self.0).for_each(|(r, s)| *r = s * rhs);
+		res
+	}
+}
+
+/// Matrix-vector product
+impl<'a, E, const M: usize, const N: usize> Mul<&'a Vector<E, M>> for &'a Vector<Vector<E, M>, N>
+where
+	E: 'a + Clone + std::fmt::Debug + Zero + Mul<&'a E, Output = E> + Add<E> + AddAssign<E>,
+{
+	type Output = Vector<E, N>;
+	fn mul(self, rhs: &'a Vector<E, M>) -> Vector<E, N> {
+		let mut res = Self::Output::zero();
+		res.0
+			.iter_mut()
+			.zip(self.0.iter())
+			.for_each(|(r, s)| *r = s * rhs);
+		res
+	}
+}
+
+impl<E, const N: usize> Vector<Vector<E, N>, N>
+where
+	E: Zero + One,
+	Vector<E, N>: AddAssign + Zero,
+{
+	pub fn id() -> Self {
+		let mut id = Self::zero();
+		for (i, a) in id.0.iter_mut().enumerate() {
+			a.0[i] = E::one();
+		}
+		id
+	}
+}
+
+#[cfg(test)]
+mod test {
+	use super::*;
+
+	use crate::{poly, ring};
+
+	use num_traits::Inv;
+
+	crate::ring!(Z50177, u32, u64, i64, 50177);
+	poly!(R50177, 256, Z50177, u32, u64, 50177, 66);
+
+	#[test]
+	fn test_basic() {
+		Vector::<R50177, 4>::zero();
+	}
+}

test.sage

@@ -0,0 +1,39 @@
+q = 50177
+
+Fq.<yq> = PolynomialRing(ZZ.quotient(q))
+Rq = Fq.quotient(yq**256+1, "xq")
+xq = Rq.gen()
+
+F2q.<y2q> = PolynomialRing(ZZ.quotient(q*2))
+R2q = F2q.quotient(y2q**256+1, "x2q")
+x2q = R2q.gen()
+
+F2.<y2> = PolynomialRing(ZZ.quotient(2))
+R2 = F2.quotient(y2**256+1, "x2")
+x2 = R2.gen()
+
+def center(x, m):
+	x = x % m
+	if x < m/2:
+		return x
+	else:
+		return x - m
+
+def f(x):
+	return (R2([int(i)%2 for i in x.list()]), Rq([int(i)%q for i in x.list()]))
+
+def g(a, b):
+	al = a.list()
+	bl = b.list()
+	return R2q([int(al[i])*q-int(bl[i])*2*(q//2) for i in range(len(al))])
+
+def add(ab, cd):
+	return (ab[0]+cd[0], ab[1]+cd[1])
+
+def mul(ab, cd):
+	return (ab[0]*cd[0], ab[1]*cd[1])
+
+a = -x2q**2 + 3*x2q + 2
+b = 4*x2q + x2q**4
+print(g(*mul(f(a), f(b))))
+print(a*b)