From 00d68d5165ff8e4ae6c57124c017e76b5d89b653 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pascal=20Eng=C3=A9libert?= Date: Mon, 21 Jul 2025 19:34:51 +0200 Subject: [PATCH] Initial commit --- .gitignore | 2 + Cargo.lock | 25 ++++ Cargo.toml | 7 + rustfmt.toml | 9 ++ src/lib.rs | 7 + src/matrix.rs | 1 + src/ntt.rs | 103 +++++++++++++++ src/poly.rs | 357 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/r7681.rs | 232 ++++++++++++++++++++++++++++++++ src/ring.rs | 133 +++++++++++++++++++ src/tuple.rs | 26 ++++ src/vector.rs | 188 ++++++++++++++++++++++++++ test.sage | 39 ++++++ 13 files changed, 1129 insertions(+) create mode 100644 .gitignore create mode 100644 Cargo.lock create mode 100644 Cargo.toml create mode 100644 rustfmt.toml create mode 100644 src/lib.rs create mode 100644 src/matrix.rs create mode 100644 src/ntt.rs create mode 100644 src/poly.rs create mode 100644 src/r7681.rs create mode 100644 src/ring.rs create mode 100644 src/tuple.rs create mode 100644 src/vector.rs create mode 100644 test.sage diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..0e5e24c --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +/target +*.sage.py diff --git a/Cargo.lock b/Cargo.lock new file mode 100644 index 0000000..3e5edc6 --- /dev/null +++ b/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", +] diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..5cf8237 --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,7 @@ +[package] +name = "gwrizienn" +version = "0.1.0" +edition = "2024" + +[dependencies] +num-traits = "0.2.19" diff --git a/rustfmt.toml b/rustfmt.toml new file mode 100644 index 0000000..26b2841 --- /dev/null +++ b/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 diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..1093f00 --- /dev/null +++ b/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; diff --git a/src/matrix.rs b/src/matrix.rs new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/src/matrix.rs @@ -0,0 +1 @@ + diff --git a/src/ntt.rs b/src/ntt.rs new file mode 100644 index 0000000..9f1952a --- /dev/null +++ b/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(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 Add> for NttDomain +where + E: Add, +{ + type Output = NttDomain; + fn add(self, rhs: NttDomain) -> Self::Output { + NttDomain(self.0 + rhs.0) + } +} + +impl<'a, E> Add<&'a NttDomain> for NttDomain +where + E: 'a + Add<&'a E, Output = E>, +{ + type Output = NttDomain; + fn add(self, rhs: &'a NttDomain) -> Self::Output { + NttDomain(self.0 + &rhs.0) + } +} + +impl Sub> for NttDomain +where + E: Sub, +{ + type Output = NttDomain; + fn sub(self, rhs: NttDomain) -> Self::Output { + NttDomain(self.0 - rhs.0) + } +} + +impl<'a, E> Sub<&'a NttDomain> for NttDomain +where + E: 'a + Sub<&'a E, Output = E>, +{ + type Output = NttDomain; + fn sub(self, rhs: &'a NttDomain) -> Self::Output { + NttDomain(self.0 - &rhs.0) + } +} + +impl Zero for NttDomain +where + E: Zero, +{ + fn zero() -> Self { + NttDomain(E::zero()) + } + fn is_zero(&self) -> bool { + self.0.is_zero() + } +} + +impl AddAssign for NttDomain +where + E: AddAssign, +{ + fn add_assign(&mut self, rhs: Self) { + self.0 += rhs.0; + } +} + +impl<'a, E> AddAssign<&'a Self> for NttDomain +where + E: 'a + AddAssign<&'a E>, +{ + fn add_assign(&mut self, rhs: &'a Self) { + self.0 += &rhs.0; + } +} + +impl SubAssign for NttDomain +where + E: SubAssign, +{ + fn sub_assign(&mut self, rhs: Self) { + self.0 -= rhs.0; + } +} + +impl<'a, E> SubAssign<&'a Self> for NttDomain +where + E: 'a + SubAssign<&'a E>, +{ + fn sub_assign(&mut self, rhs: &'a Self) { + self.0 -= &rhs.0; + } +} diff --git a/src/poly.rs b/src/poly.rs new file mode 100644 index 0000000..320df18 --- /dev/null +++ b/src/poly.rs @@ -0,0 +1,357 @@ +use crate::{ntt::*, ring::Ring}; + +use num_traits::Zero; + +pub trait Poly { + 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); + } +} diff --git a/src/r7681.rs b/src/r7681.rs new file mode 100644 index 0000000..fb975b8 --- /dev/null +++ b/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(pub [Z7681; N]); + + impl crate::poly::Poly for R7681 { + type Element = Z7681; + fn modulus() -> u16 { + 7681 + } + fn new(v: [Z7681; N]) -> Self { + Self(v) + } + } + + impl R7681 { + 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 std::ops::Add for R7681 { + 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 R7681 { + 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 R7681 { + 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 R7681 { + 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 R7681 { + 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 R7681 { + 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 R7681 { + 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 R7681 { + fn sub_assign(&mut self, rhs: &Self) { + self.0 + .iter_mut() + .zip(rhs.0.iter()) + .for_each(|(x, y)| *x -= *y); + } + } + + impl num_traits::Zero for R7681 { + 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 { + 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 R7681 { + type Output = crate::ntt::NttDomain>; + fn ntt(mut self) -> crate::ntt::NttDomain> { + 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::::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> { + type Output = R7681; + fn ntt_inv(mut self) -> R7681 { + 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::::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 std::ops::Mul for crate::ntt::NttDomain> { + 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> { + 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> { + 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> { + 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> { + type Output = Self; + fn inv(mut self) -> Self { + self.0.0.iter_mut().for_each(|x| *x = x.inverse()); + self + } + } diff --git a/src/ring.rs b/src/ring.rs new file mode 100644 index 0000000..dc2b5cb --- /dev/null +++ b/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()); + } +} diff --git a/src/tuple.rs b/src/tuple.rs new file mode 100644 index 0000000..27482a8 --- /dev/null +++ b/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(pub A, pub B); + +impl Add for Tuple where A: Add, B: Add { + type Output = Self; + fn add(self, rhs: Self) -> Self { + Self(self.0+rhs.0, self.1+rhs.1) + } +} + +impl Mul for Tuple where A: Mul, B: Mul { + type Output = Self; + fn mul(self, rhs: Self) -> Self { + Self(self.0*rhs.0, self.1*rhs.1) + } +} + +impl Inv for Tuple where A: Inv, B: Inv { + type Output = Self; + fn inv(self) -> Self { + Self(self.0.inv(), self.1.inv()) + } +} diff --git a/src/vector.rs b/src/vector.rs new file mode 100644 index 0000000..4b0ce28 --- /dev/null +++ b/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(pub [E; N]); + +/*impl Vector { + pub fn norm_inf(&self) -> Zq { + Zq(self.0.iter().map(|i| i.norm_inf().0).max().unwrap()) + } +}*/ + +impl Zero for Vector +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 Add for Vector +where + E: AddAssign, +{ + 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 +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 Sub for Vector +where + E: SubAssign, +{ + 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 +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 +where + E: 'a + Clone + std::fmt::Debug + Zero + MulAssign<&'a E>, +{ + type Output = Vector; + fn mul(mut self, rhs: &'a E) -> Vector { + 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 +where + E: 'a + Clone + Zero + MulAssign<&'a E> + Add, +{ + 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 for &'a Vector +where + E: 'a + Clone + Zero + Mul<&'a E, Output = E> + Add, +{ + 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> for Vector, N> +where + E: 'a + Clone + std::fmt::Debug + Zero + MulAssign<&'a E> + Add + AddAssign, +{ + type Output = Vector; + fn mul(self, rhs: &'a Vector) -> Vector { + 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> for &'a Vector, N> +where + E: 'a + Clone + std::fmt::Debug + Zero + Mul<&'a E, Output = E> + Add + AddAssign, +{ + type Output = Vector; + fn mul(self, rhs: &'a Vector) -> Vector { + let mut res = Self::Output::zero(); + res.0 + .iter_mut() + .zip(self.0.iter()) + .for_each(|(r, s)| *r = s * rhs); + res + } +} + +impl Vector, N> +where + E: Zero + One, + Vector: 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::::zero(); + } +} diff --git a/test.sage b/test.sage new file mode 100644 index 0000000..fe47d9a --- /dev/null +++ b/test.sage @@ -0,0 +1,39 @@ +q = 50177 + +Fq. = PolynomialRing(ZZ.quotient(q)) +Rq = Fq.quotient(yq**256+1, "xq") +xq = Rq.gen() + +F2q. = PolynomialRing(ZZ.quotient(q*2)) +R2q = F2q.quotient(y2q**256+1, "x2q") +x2q = R2q.gen() + +F2. = 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)