diff --git a/.gitignore b/.gitignore index 0e5e24c..ea8c4bf 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1 @@ /target -*.sage.py diff --git a/Cargo.lock b/Cargo.lock index 3e5edc6..8c86fea 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -8,13 +8,71 @@ version = "1.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "c08606f8c3cbf4ce6ec8e28fb0014a2c086708fe954eaa885384a6165172e7e8" +[[package]] +name = "cfg-if" +version = "1.0.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2fd1289c04a9ea8cb22300a459a72a385d7c73d3259e2ed7dcb2af674838cfa9" + +[[package]] +name = "crossbeam-deque" +version = "0.8.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9dd111b7b7f7d55b72c0a6ae361660ee5853c9af73f70c3c2ef6858b950e2e51" +dependencies = [ + "crossbeam-epoch", + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-epoch" +version = "0.9.18" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e" +dependencies = [ + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-utils" +version = "0.8.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d0a5c400df2834b80a4c3327b3aad3a4c4cd4de0629063962b03235697506a28" + +[[package]] +name = "either" +version = "1.15.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719" + +[[package]] +name = "getrandom" +version = "0.2.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "335ff9f135e4384c8150d6f27c6daed433577f86b4750418338c01a1a2528592" +dependencies = [ + "cfg-if", + "libc", + "wasi", +] + [[package]] name = "gwrizienn" version = "0.1.0" dependencies = [ "num-traits", + "rand", + "rand_core", + "rayon", + "zeroize", ] +[[package]] +name = "libc" +version = "0.2.175" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6a82ae493e598baaea5209805c49bbf2ea7de956d50d7da0da1164f9c6d28543" + [[package]] name = "num-traits" version = "0.2.19" @@ -23,3 +81,129 @@ checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" dependencies = [ "autocfg", ] + +[[package]] +name = "ppv-lite86" +version = "0.2.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85eae3c4ed2f50dcfe72643da4befc30deadb458a9b590d720cde2f2b1e97da9" +dependencies = [ + "zerocopy", +] + +[[package]] +name = "proc-macro2" +version = "1.0.101" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "89ae43fd86e4158d6db51ad8e2b80f313af9cc74f5c0e03ccb87de09998732de" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "quote" +version = "1.0.40" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1885c039570dc00dcb4ff087a89e185fd56bae234ddc7f056a945bf36467248d" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "rand" +version = "0.8.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "34af8d1a0e25924bc5b7c43c079c942339d8f0a8b57c39049bef581b46327404" +dependencies = [ + "libc", + "rand_chacha", + "rand_core", +] + +[[package]] +name = "rand_chacha" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88" +dependencies = [ + "ppv-lite86", + "rand_core", +] + +[[package]] +name = "rand_core" +version = "0.6.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c" +dependencies = [ + "getrandom", +] + +[[package]] +name = "rayon" +version = "1.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "368f01d005bf8fd9b1206fb6fa653e6c4a81ceb1466406b81792d87c5677a58f" +dependencies = [ + "either", + "rayon-core", +] + +[[package]] +name = "rayon-core" +version = "1.13.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "22e18b0f0062d30d4230b2e85ff77fdfe4326feb054b9783a3460d8435c8ab91" +dependencies = [ + "crossbeam-deque", + "crossbeam-utils", +] + +[[package]] +name = "syn" +version = "2.0.106" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ede7c438028d4436d71104916910f5bb611972c5cfd7f89b8300a8186e6fada6" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "unicode-ident" +version = "1.0.18" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5a5f39404a5da50712a4c1eecf25e90dd62b613502b7e925fd4e4d19b5c96512" + +[[package]] +name = "wasi" +version = "0.11.1+wasi-snapshot-preview1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ccf3ec651a847eb01de73ccad15eb7d99f80485de043efb2f370cd654f4ea44b" + +[[package]] +name = "zerocopy" +version = "0.8.26" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1039dd0d3c310cf05de012d8a39ff557cb0d23087fd44cad61df08fc31907a2f" +dependencies = [ + "zerocopy-derive", +] + +[[package]] +name = "zerocopy-derive" +version = "0.8.26" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9ecf5b4cc5364572d7f4c329661bcc82724222973f2cab6f050a4e5c22f75181" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "zeroize" +version = "1.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ced3678a2879b30306d323f4542626697a464a97c0a07c9aebf7ebca65cd4dde" diff --git a/Cargo.toml b/Cargo.toml index 5cf8237..7570ae1 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -4,4 +4,16 @@ version = "0.1.0" edition = "2024" [dependencies] -num-traits = "0.2.19" +num-traits = "0.2" +rand = { version = "0.8", optional = true } +rand_core = { version = "0.6", optional = true } +zeroize = { version = "1", optional = true, default-features = false } + +[dev-dependencies] +rayon = "1.8" + +[features] +default = ["rand", "zeroize"] + +rand = ["dep:rand", "dep:rand_core"] +zeroize = ["dep:zeroize"] diff --git a/README.md b/README.md index 0395d6d..df07f47 100644 --- a/README.md +++ b/README.md @@ -12,10 +12,17 @@ Goals: Non-goals: * Generic (it's generic but only for primitive types) * Complete (it's simple because it's not a complete implementation of modern algebra) +* Dynamic (vector dimensions and moduli are strongly typed) + +Supported: +* ring Zq +* ring Zq/(x^N+1) with additive operations +* ring Zq/(x^N+1) with multiplicative operations if q=p or q=2p with p prime and 2N divides p-1 +* vectors and matrices of the above rings ## Name -It's Breton for "root", because we use roots of unity to compute the NTT for faster polynomial multiplication. +Pronounce _grizienn_. It's Breton for "root", because we use roots of unity to compute the NTT for faster O(N log N) polynomial multiplication. ## License diff --git a/examples/dilithium.rs b/examples/dilithium.rs new file mode 100644 index 0000000..b87aefa --- /dev/null +++ b/examples/dilithium.rs @@ -0,0 +1,77 @@ +//! Example of simplified Dilithium with q=8380417 + +use gwrizienn::{ + matrix::Matrix, + ntt::{Ntt, NttInv}, + ring::Ring, + vector::Vector, + *, +}; + +// Implement Zq +ring!(Zq, u32, u64, i64, 8380417); + +// Implement Rq = Zq/(x^256+1) +// zeta=1753 is the first 512-th root of unity mod q +poly!(Rq, 256, Zq, u32, u64, u32, u64, 8380417, 1753); + +fn high_bits(mut v: Vector) -> Vector { + for vi in v.0.iter_mut() { + for vij in vi.0.iter_mut() { + vij.0 -= vij.0 % 190464; + } + } + v +} + +fn main() { + let mut rng = rand::thread_rng(); + let uniform = Zq::uniform(); + let ball_c = Zq::uniform_ball(1); + let ball_s = Zq::uniform_ball(2); + let ball_y = Zq::uniform_ball(131071); + + // generate secret key + let a = Matrix::::random(uniform, &mut rng); + let s1 = Vector::::random(ball_s, &mut rng); + let s2 = Vector::::random(ball_s, &mut rng); + // random value for signing + let y = Vector::::random(ball_y, &mut rng); + // challenge + let c = Rq::random(ball_c, &mut rng); + + // use NTT for fast multiplication + let a = a.ntt(); + let s1 = s1.ntt(); + let s2 = s2.ntt(); + let y = y.ntt(); + let c = c.ntt(); + + // generate public key + let t = &a * &s1 + s2; + // commitment + let w = &a * &y; + // proof + let z = y + s1 * &c; + + // verify + assert_eq!( + high_bits((&a * &z - t * &c).ntt_inv()), + high_bits(w.ntt_inv()) + ); + + // let uniform = Zq::uniform(); + // let ball_c = Zq::uniform_ball(1); + // let ball_s = Zq::uniform_ball(2); + // let ball_y = Zq::uniform_ball(131071); + // + // let a = Matrix::::random(uniform, &mut rng).ntt(); + // let s1 = Vector::::random(ball_s, &mut rng).ntt(); + // let s2 = Vector::::random(ball_s, &mut rng).ntt(); + // let y = Vector::::random(ball_y, &mut rng).ntt(); + // let c = Rq::random(ball_c, &mut rng).ntt(); + // + // let t = &a * &s1 + s2; + // let w = &a * &y; + // let z = y + s1 * &c; +} diff --git a/examples/ntwe.rs b/examples/ntwe.rs new file mode 100644 index 0000000..e4f1a37 --- /dev/null +++ b/examples/ntwe.rs @@ -0,0 +1,101 @@ +//! Example of simplified NTWE scheme (Gärtner 2024) + +use gwrizienn::{ + matrix::Matrix, + ntt::{Ntt, NttInv}, + ring::{Lift, Ring}, + vector::{Vector, VectorRef}, + *, +}; +use num_traits::{Inv, One, Zero}; + +// Implement Zq +ring!(Zq, u32, u32, i32, 50177); + +// Implement Rq = Zq/(x^256+1) +// zeta=66 is the first 512-th root of unity mod q +poly!(Rq, 256, Zq, u32, u32, u32, u64, 50177, 66); + +// Implement Zb and Rb with a big modulus for NTT in Z2q +// chosen because > N*q^2 +ring!(Zb, u64, u128, i128, 644539222529); +poly!( + Rb, + 256, + Zb, + u64, + u128, + u64, + u128, + 644539222529, + 483489047161 +); + +//Implement Z2q +ring!(Z2q, u32, u64, i64, 100354); + +// Implement R2q +// A different macro must be used because 2q is not prime. +poly2!( + Rq, + Rb, + R2q, + 256, + Z2q, + u32, + u64, + u32, + u64, + i64, + i64, + 50177, + 644539222529 +); + +const L: usize = 3; +const M: usize = 2; + +fn main() { + let mut rng = rand::thread_rng(); + let uniform = Zq::uniform(); + let ball_c = Z2q::uniform_positive_semiball(1); + let ball_s = Zq::uniform_ball(1); + let ball_y = Zb::uniform_ball(55); + + // generate secret key s = [f s0 e] + let mut s = Vector::::random(ball_s, &mut rng); + let f0 = &mut s.0[0]; + *f0 *= Zq(2); + *f0 += Zq::one(); + let s_ntt = s.clone().ntt(); + let f = &s_ntt.0[0]; + let s0: VectorRef<_, L> = s_ntt.get_sub(1); + let e: VectorRef<_, M> = s_ntt.get_sub(1 + L); + // generate public key + let a0 = Matrix::::random(uniform, &mut rng); + let b = (&a0.clone().ntt() * s0 + e) * &f.clone().inv(); + let mut a = Matrix::::zero(); + a.set_column(0, (b.ntt_inv().lift() * Z2q(100352)).get_ref()); + a[0][0] += Z2q(50177); + a.set_columns(1, &(a0.lift() * Z2q(2))); + for i in 0..M { + a[i][1 + L + i] += Z2q(2); + } + let a: Matrix = a.lift(); + let a = a.ntt(); + // random value for signing + let y = Vector::::random(ball_y, &mut rng); + let y = y.ntt(); + // commitment + let w: Vector = (&a * &y).ntt_inv().lift(); + // challenge + let c = R2q::random(ball_c, &mut rng); + let cb: Rb = c.lift(); + // proof + let s: Vector = s.lift().lift(); + let z = y + s.ntt() * &cb.ntt(); + // verify + let mut w2: Vector = (&a * &z).ntt_inv().lift(); + w2[0] -= c * Z2q(50177); + assert_eq!(w, w2); +} diff --git a/src/lib.rs b/src/lib.rs index 1093f00..da33889 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,7 +1,33 @@ +#![no_std] +#![warn(missing_docs)] +#![deny(non_ascii_idents)] +#![deny(unnameable_types)] +#![deny(unreachable_pub)] +#![deny(unstable_features)] +#![warn(unused_qualifications)] +#![allow(clippy::tabs_in_doc_comments)] + +//! Modular and polynomial arithmetic. + pub mod matrix; pub mod ntt; pub mod poly; -//pub mod r7681; pub mod ring; -pub mod tuple; pub mod vector; + +/// Something that can be sampled using a distribution D over T +#[cfg(feature = "rand")] +pub trait Random, T> { + /// Sample an element of T from distribution `distr` + fn random(distr: D, rng: &mut R) -> Self; +} + +/// Type having a unity, or multiplicative neutral element +/// +/// This trait is useful because `num_traits::One` requires `Mul`, which is not implemented for polynomials. +pub trait One { + /// the unity + fn one() -> Self; + /// is it the unity? + fn is_one(&self) -> bool; +} diff --git a/src/matrix.rs b/src/matrix.rs index 8b13789..3bca0f7 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -1 +1,401 @@ +//! Matrices of ring elements or polynomials. +use crate::vector::{Dot, Vector, VectorRef}; + +use core::{ + mem::MaybeUninit, + ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign}, +}; +use num_traits::Zero; + +/// Matrix of E of M rows and N columns +#[derive(Clone, Debug, Eq, Hash, PartialEq)] +pub struct Matrix(pub [Vector; M]); + +/// Referenced matrix of E of M rows and N columns +#[derive(Debug, Eq, Hash, PartialEq)] +pub struct MatrixRef<'a, E, const M: usize, const N: usize>(pub &'a [Vector; M]); + +/// Mutable referenced matrix of E of M rows and N columns +#[derive(Debug, Eq, Hash, PartialEq)] +pub struct MatrixMut<'a, E, const M: usize, const N: usize>(pub &'a mut [Vector; M]); + +impl<'a, E, const M: usize, const N: usize> Clone for MatrixRef<'a, E, M, N> { + fn clone(&self) -> Self { + *self + } +} + +impl<'a, E, const M: usize, const N: usize> Copy for MatrixRef<'a, E, M, N> {} + +impl Matrix { + /// Set column `idx` + pub fn set_column(&mut self, idx: usize, v: VectorRef) { + for (row, x) in self.0.iter_mut().zip(v.0.iter()) { + row.0[idx] = x.clone(); + } + } + /// Set `S` columns starting at column `idx` to values from the matrix `v` + pub fn set_columns(&mut self, idx: usize, v: &Matrix) { + for (s_row, v_row) in self.0.iter_mut().zip(v.0.iter()) { + for (s_x, v_x) in s_row.0[idx..idx + S].iter_mut().zip(v_row.0.iter()) { + *s_x = v_x.clone(); + } + } + } +} + +impl Zero for Matrix +where + Vector: Zero, + Self: Add, +{ + fn zero() -> Self { + let mut z = MaybeUninit::<[Vector; M]>::uninit(); + unsafe { + for i in z.assume_init_mut() { + *i = Vector::zero(); + } + } + Self(unsafe { z.assume_init() }) + } + fn is_zero(&self) -> bool { + self.0.iter().all(Zero::is_zero) + } +} + +impl crate::One for Matrix +where + E: crate::One + Zero, + Self: Zero, +{ + fn one() -> Self { + let mut id = Self::zero(); + for (i, a) in id.0.iter_mut().enumerate() { + a.0[i] = E::one(); + } + id + } + fn is_one(&self) -> bool { + self.0.iter().enumerate().all(|(i, row)| { + row.0 + .iter() + .enumerate() + .all(|(j, v)| if i == j { v.is_one() } else { v.is_zero() }) + }) + } +} + +impl Add for Matrix +where + Vector: 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 M: usize, const N: usize> Add<&'a Self> for Matrix +where + Vector: 'a + AddAssign<&'a Vector>, +{ + 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<'a, E, const M: usize, const N: usize> Add> for Matrix +where + Vector: 'a + AddAssign<&'a Vector>, +{ + type Output = Self; + fn add(mut self, rhs: MatrixRef<'a, E, M, N>) -> Self { + self.0 + .iter_mut() + .zip(rhs.0.iter()) + .for_each(|(x, y)| *x += y); + self + } +} + +impl AddAssign for Matrix +where + Vector: AddAssign>, +{ + fn add_assign(&mut self, rhs: Self) { + self.0.iter_mut().zip(rhs.0).for_each(|(x, y)| *x += y); + } +} + +impl<'a, E, const M: usize, const N: usize> AddAssign<&'a Self> for Matrix +where + Vector: 'a + AddAssign<&'a Vector>, +{ + fn add_assign(&mut self, rhs: &'a Self) { + self.0 + .iter_mut() + .zip(rhs.0.iter()) + .for_each(|(x, y)| *x += y); + } +} + +impl<'a, E, const M: usize, const N: usize> AddAssign> for Matrix +where + Vector: 'a + AddAssign<&'a Vector>, +{ + fn add_assign(&mut self, rhs: MatrixRef<'a, E, M, N>) { + self.0 + .iter_mut() + .zip(rhs.0.iter()) + .for_each(|(x, y)| *x += y); + } +} + +impl Sub for Matrix +where + Vector: 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 SubAssign for Matrix +where + Vector: SubAssign>, +{ + fn sub_assign(&mut self, rhs: Self) { + self.0.iter_mut().zip(rhs.0).for_each(|(x, y)| *x -= y); + } +} + +impl<'a, E, const M: usize, const N: usize> Sub<&'a Self> for Matrix +where + Vector: 'a + SubAssign<&'a Vector>, +{ + 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 + } +} + +impl<'a, E, const M: usize, const N: usize> SubAssign<&'a Self> for Matrix +where + Vector: 'a + SubAssign<&'a Vector>, +{ + fn sub_assign(&mut self, rhs: &'a Self) { + self.0 + .iter_mut() + .zip(rhs.0.iter()) + .for_each(|(x, y)| *x -= y); + } +} + +impl<'a, E, const M: usize, const N: usize> Sub> for Matrix +where + Vector: 'a + SubAssign<&'a Vector>, +{ + type Output = Self; + fn sub(mut self, rhs: MatrixRef<'a, E, M, N>) -> Self { + self.0 + .iter_mut() + .zip(rhs.0.iter()) + .for_each(|(x, y)| *x -= y); + self + } +} + +impl<'a, E, const M: usize, const N: usize> SubAssign> for Matrix +where + Vector: 'a + SubAssign<&'a Vector>, +{ + fn sub_assign(&mut self, rhs: MatrixRef<'a, E, M, N>) { + self.0 + .iter_mut() + .zip(rhs.0.iter()) + .for_each(|(x, y)| *x -= y); + } +} + +impl crate::ntt::Ntt for Matrix +where + E: crate::ntt::Ntt, +{ + type Output = Matrix; + fn ntt(self) -> Matrix { + let mut z = MaybeUninit::<[Vector; M]>::uninit(); + unsafe { + for (i, j) in z.assume_init_mut().iter_mut().zip(self.0) { + *i = j.ntt(); + } + } + Matrix(unsafe { z.assume_init() }) + } +} + +impl crate::ntt::NttInv for Matrix +where + E: crate::ntt::NttInv, +{ + type Output = Matrix; + fn ntt_inv(self) -> Matrix { + let mut z = MaybeUninit::<[Vector; M]>::uninit(); + unsafe { + for (i, j) in z.assume_init_mut().iter_mut().zip(self.0) { + *i = j.ntt_inv(); + } + } + Matrix(unsafe { z.assume_init() }) + } +} + +// TODO only for safe in-place liftable elements +impl crate::ring::Lift, T, U> + for Matrix +where + Fr: crate::ring::Lift, +{ + fn lift(self) -> Matrix { + let mut z = MaybeUninit::<[Vector; M]>::uninit(); + unsafe { + for (i, j) in z.assume_init_mut().iter_mut().zip(self.0) { + *i = j.lift(); + } + } + Matrix(unsafe { z.assume_init() }) + } +} + +impl Mul for Matrix +where + E: MulAssign, + F: Copy, +{ + type Output = Self; + fn mul(mut self, rhs: F) -> Self { + for row in self.0.iter_mut() { + for v in row.0.iter_mut() { + *v *= rhs; + } + } + self + } +} + +impl<'a, E, const M: usize, const N: usize> Mul<&'a Vector> for &'a Matrix +where + &'a Vector: Dot<&'a Vector, Output = E>, + Vector: Zero, +{ + 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.dot(rhs)); + res + } +} + +impl<'a, E, const M: usize, const N: usize> Mul> for &'a Matrix +where + &'a Vector: Dot, Output = E>, + Vector: Zero, +{ + type Output = Vector; + fn mul(self, rhs: VectorRef<'a, E, N>) -> Vector { + let mut res = Self::Output::zero(); + res.0 + .iter_mut() + .zip(self.0.iter()) + .for_each(|(r, s)| *r = s.dot(rhs)); + res + } +} + +#[cfg(feature = "rand")] +impl crate::Random for Matrix +where + D: Clone + rand::distributions::Distribution, + E: crate::Random, +{ + fn random(distr: D, rng: &mut R) -> Self { + let mut z = MaybeUninit::<[Vector; M]>::uninit(); + unsafe { + for i in z.assume_init_mut() { + *i = Vector::random(distr.clone(), rng); + } + } + Self(unsafe { z.assume_init() }) + } +} + +impl core::ops::Index for Matrix { + type Output = Vector; + fn index(&self, idx: usize) -> &Vector { + &self.0[idx] + } +} + +impl core::ops::IndexMut for Matrix { + fn index_mut(&mut self, idx: usize) -> &mut Vector { + &mut self.0[idx] + } +} + +#[cfg(feature = "zeroize")] +impl zeroize::Zeroize for Matrix { + fn zeroize(&mut self) { + self.0.zeroize(); + } +} + +impl core::fmt::Display for Matrix { + fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result { + write!(f, "[")?; + for i in self.0.iter() { + write!(f, "{},", i)?; + } + write!(f, "]") + } +} + +#[cfg(test)] +#[allow(unreachable_pub)] +#[allow(unnameable_types)] +mod test { + use super::*; + + use crate::{ntt::Ntt, poly, ring}; + + ring!(Z50177, u32, u32, i64, 50177); + poly!(R50177, 256, Z50177, u32, u32, u32, u64, 50177, 66); + + #[test] + fn test_basic() { + Matrix::::zero(); + } + + #[test] + fn test_mul() { + let a = Matrix::::zero(); + let b = Vector::::zero(); + let x = R50177::one(); + let _ = a.clone().ntt() * &x.ntt(); + let _ = &a.clone().ntt() * &b.ntt(); + } +} diff --git a/src/ntt.rs b/src/ntt.rs index 9f1952a..91d17ba 100644 --- a/src/ntt.rs +++ b/src/ntt.rs @@ -1,72 +1,84 @@ -use num_traits::Zero; -use std::ops::{Add, AddAssign, Sub, SubAssign}; +//! Number-Theoretic Transform (NTT) is kind of a discrete Fourier Transform that transforms a polynomial to the NTT domain, where multiplication can be done efficiently in O(N) instead of O(N^2). Addition can also be performed efficiently in the NTT domain. The `NttDomain` struct ensures mathematical soundness and domains are not used together. +use core::{ + marker::PhantomData, + ops::{Add, AddAssign, Index, IndexMut, Mul, MulAssign, Sub, SubAssign}, +}; +use num_traits::{Inv, Zero}; + +/// Element E in the NTT domain #[derive(Debug, Clone, Eq, PartialEq)] -pub struct NttDomain(pub E); +pub struct NttDomain(pub E, pub PhantomData<[EE; N]>); +/// An element that can be transformed by the NTT pub trait Ntt: Sized { + /// Type in the NTT domain type Output: Sized; + /// Apply the NTT fn ntt(self) -> Self::Output; } +/// An element in the NTT domain that can be transformed back by the inverse NTT pub trait NttInv: Sized { + /// Type out of the NTT domain type Output: Sized; + /// Apply the inverse NTT fn ntt_inv(self) -> Self::Output; } -impl Add> for NttDomain +impl Add> for NttDomain where E: Add, { - type Output = NttDomain; - fn add(self, rhs: NttDomain) -> Self::Output { - NttDomain(self.0 + rhs.0) + type Output = NttDomain; + fn add(self, rhs: NttDomain) -> Self::Output { + NttDomain(self.0 + rhs.0, Default::default()) } } -impl<'a, E> Add<&'a NttDomain> for NttDomain +impl<'a, E, EE, const N: usize> 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) + type Output = NttDomain; + fn add(self, rhs: &'a NttDomain) -> Self::Output { + NttDomain(self.0 + &rhs.0, Default::default()) } } -impl Sub> for NttDomain +impl Sub> for NttDomain where E: Sub, { - type Output = NttDomain; - fn sub(self, rhs: NttDomain) -> Self::Output { - NttDomain(self.0 - rhs.0) + type Output = NttDomain; + fn sub(self, rhs: NttDomain) -> Self::Output { + NttDomain(self.0 - rhs.0, Default::default()) } } -impl<'a, E> Sub<&'a NttDomain> for NttDomain +impl<'a, E, EE, const N: usize> 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) + type Output = NttDomain; + fn sub(self, rhs: &'a NttDomain) -> Self::Output { + NttDomain(self.0 - &rhs.0, Default::default()) } } -impl Zero for NttDomain +impl Zero for NttDomain where E: Zero, { fn zero() -> Self { - NttDomain(E::zero()) + NttDomain(E::zero(), Default::default()) } fn is_zero(&self) -> bool { self.0.is_zero() } } -impl AddAssign for NttDomain +impl AddAssign for NttDomain where E: AddAssign, { @@ -75,7 +87,7 @@ where } } -impl<'a, E> AddAssign<&'a Self> for NttDomain +impl<'a, E, EE, const N: usize> AddAssign<&'a Self> for NttDomain where E: 'a + AddAssign<&'a E>, { @@ -84,7 +96,7 @@ where } } -impl SubAssign for NttDomain +impl SubAssign for NttDomain where E: SubAssign, { @@ -93,7 +105,7 @@ where } } -impl<'a, E> SubAssign<&'a Self> for NttDomain +impl<'a, E, EE, const N: usize> SubAssign<&'a Self> for NttDomain where E: 'a + SubAssign<&'a E>, { @@ -101,3 +113,207 @@ where self.0 -= &rhs.0; } } + +impl Mul for NttDomain +where + E: crate::poly::Poly, + EE: Copy + MulAssign, +{ + type Output = Self; + fn mul(mut self, rhs: Self) -> Self { + self.0 + .coefficients_mut() + .iter_mut() + .zip(rhs.0.coefficients().iter()) + .for_each(|(x, y)| *x *= *y); + self + } +} + +impl Mul<&Self> for NttDomain +where + E: crate::poly::Poly, + EE: Copy + MulAssign, +{ + type Output = Self; + fn mul(mut self, rhs: &Self) -> Self { + self.0 + .coefficients_mut() + .iter_mut() + .zip(rhs.0.coefficients().iter()) + .for_each(|(x, y)| *x *= *y); + self + } +} + +impl MulAssign for NttDomain +where + E: crate::poly::Poly, + EE: Copy + MulAssign, +{ + fn mul_assign(&mut self, rhs: Self) { + self.0 + .coefficients_mut() + .iter_mut() + .zip(rhs.0.coefficients().iter()) + .for_each(|(x, y)| *x *= *y); + } +} + +impl MulAssign<&Self> for NttDomain +where + E: crate::poly::Poly, + EE: Copy + MulAssign, +{ + fn mul_assign(&mut self, rhs: &Self) { + self.0 + .coefficients_mut() + .iter_mut() + .zip(rhs.0.coefficients().iter()) + .for_each(|(x, y)| *x *= *y); + } +} + +impl Mul for NttDomain +where + E: crate::poly::Poly, + EE: Copy + MulAssign, +{ + type Output = Self; + fn mul(mut self, rhs: EE) -> Self { + self.0.coefficients_mut().iter_mut().for_each(|x| *x *= rhs); + self + } +} + +impl MulAssign for NttDomain +where + E: crate::poly::Poly, + EE: Copy + MulAssign, +{ + fn mul_assign(&mut self, rhs: EE) { + self.0.coefficients_mut().iter_mut().for_each(|x| *x *= rhs); + } +} + +impl Inv for NttDomain +where + E: crate::poly::Poly, + EE: Copy + Inv, +{ + type Output = Self; + fn inv(mut self) -> Self { + self.0 + .coefficients_mut() + .iter_mut() + .for_each(|x| *x = x.inv()); + self + } +} + +impl NttInv for NttDomain +where + E: crate::poly::InternalNttInv, +{ + type Output = E; + fn ntt_inv(self) -> E { + E::__ntt_inv(self.0) + } +} + +impl crate::poly::Poly for NttDomain +where + E: crate::poly::Poly, + EE: crate::ring::Ring<>::T> + Copy, + >::T: Clone, +{ + type T = >::T; + type Element = EE; + + fn coefficients(&self) -> &[EE; N] { + self.0.coefficients() + } + + fn coefficients_mut(&mut self) -> &mut [EE; N] { + self.0.coefficients_mut() + } + + fn from_scalar(e: EE) -> Self { + Self(E::from_coefficients([e; N]), Default::default()) + } + + fn to_coefficients(self) -> [EE; N] { + self.0.to_coefficients() + } + + fn from_coefficients(v: [EE; N]) -> Self { + Self(E::from_coefficients(v), Default::default()) + } +} + +#[cfg(feature = "rand")] +impl crate::Random, T> + for NttDomain +where + E: crate::Random, T>, + crate::ring::UniformZq: rand::distributions::Distribution, +{ + fn random(distr: crate::ring::UniformZq, rng: &mut R) -> Self { + Self(E::random(distr, rng), Default::default()) + } +} + +impl Index for NttDomain +where + E: Index, +{ + type Output = EE; + fn index(&self, idx: usize) -> &EE { + &self.0[idx] + } +} + +impl IndexMut for NttDomain +where + E: Index + IndexMut, +{ + fn index_mut(&mut self, idx: usize) -> &mut EE { + &mut self.0[idx] + } +} + +#[cfg(feature = "zeroize")] +impl zeroize::Zeroize for NttDomain { + fn zeroize(&mut self) { + self.0.zeroize(); + } +} + +impl core::fmt::Display for NttDomain { + fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result { + self.0.fmt(f) + } +} + +/// Compute self + a * b in a way that may be faster +pub trait Fma { + /// Compute self + a * b in a way that may be faster + /// + /// You are encouraged to benchmark this against the naive operations. + /// Compiler optimizations may render this function uninteresting. + fn fma(&mut self, a: A, b: B); +} + +impl Fma<&Self, &Self> for NttDomain +where + E: crate::poly::Poly, + EE: Copy + Mul + AddAssign, +{ + fn fma(&mut self, a: &Self, b: &Self) { + self.0 + .coefficients_mut() + .iter_mut() + .zip(a.0.coefficients().iter().zip(b.0.coefficients().iter())) + .for_each(|(si, (ai, bi))| *si += *ai * *bi); + } +} diff --git a/src/poly.rs b/src/poly.rs index 320df18..ae6d106 100644 --- a/src/poly.rs +++ b/src/poly.rs @@ -1,36 +1,106 @@ -use crate::{ntt::*, ring::Ring}; +//! Macros and traits for polynomial quotient ring over modular rings, (Z/qZ)\[x\]/(x^N+1). -use num_traits::Zero; +use core::ops::{Index, IndexMut}; -pub trait Poly { - type Element: Ring; +/// Polynomial of degree at most N-1 +pub trait Poly: Index + IndexMut { + /// Type of the representatives of the coefficients + type T; + /// Type of the coefficients + type Element: crate::ring::Ring; - fn modulus() -> T; - fn new(v: [Self::Element; N]) -> Self; + /// Get the list of coefficients (least degree first) + fn coefficients(&self) -> &[Self::Element; N]; + + /// Get the list of coefficients (least degree first) + fn coefficients_mut(&mut self) -> &mut [Self::Element; N]; + + /// New polynomial from a scalar (safely reducing) + fn from_scalar(e: Self::Element) -> Self; + + /// Get the list of coefficients (least degree first) + fn to_coefficients(self) -> [Self::Element; N]; + + /// New polynomial from the list of coefficients (least degree first) + fn from_coefficients(v: [Self::Element; N]) -> Self; } -pub(crate) const fn bitrev(x: usize, n: usize) -> usize { +/// Internal helper trait for inverse NTT (needed for macros) +#[doc(hidden)] +pub trait InternalNttInv { + type Input; + type Output; + fn __ntt_inv(input: Self::Input) -> Self::Output; +} + +/// Reverse the first n bits of x +#[doc(hidden)] +pub const fn bitrev(x: usize, n: usize) -> usize { x.reverse_bits() >> (usize::BITS - n.ilog2()) } +/// Implement basic traits for a modular polynomial ring +/// +/// - `$Rq`: the struct that will be created, that represents a polynomial in (Z/qZ)/(x^N+1) +/// - `$N`: degree N +/// - `$Zq`: (defined) underlying modular ring Z/qZ +/// - `$T`: (defined) underlying scalar type +/// - `$TT`: (defined) scalar type big enough to represent the square of the biggest element of Zq +/// - `$TN`: (defined) scalar type big enough to represent N times half the biggest element of Zq +/// - `$TTN`: (defined) scalar type big enough to represent N times the square of half the biggest element of Zq +/// - `$q`: modulus #[macro_export] +#[doc(hidden)] 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) + ($Rq:ident, $N:literal, $Zq:ident, $T:ident, $TT:ident, $TN:ident, $TTN:ident, $q:expr) => { + /// Element of $Z_q[x]/(x^N+1)$ #[derive(Clone, Copy, Debug, Eq, Hash, PartialEq)] + #[repr(transparent)] pub struct $Rq(pub [$Zq; $N]); - impl $crate::poly::Poly<$T, $N> for $Rq { + impl $crate::poly::Poly<$N> for $Rq { + type T = $T; type Element = $Zq; - fn modulus() -> $T { - $q + + fn coefficients(&self) -> &[$Zq; $N] { + &self.0 } - fn new(v: [$Zq; $N]) -> Self { + + fn coefficients_mut(&mut self) -> &mut [$Zq; $N] { + &mut self.0 + } + + fn from_scalar(e: $Zq) -> Self { + let mut v = [$Zq(0); $N]; + v[0] = e; + Self(v) + } + + fn to_coefficients(self) -> [$Zq; $N] { + self.0 + } + + fn from_coefficients(v: [$Zq; $N]) -> Self { Self(v) } } - impl std::ops::Add for $Rq { + //#[cfg(feature = "rand")] + impl $crate::Random for $Rq + where + D: Clone + rand::distributions::Distribution, + $Zq: $crate::Random, + { + fn random(distr: D, rng: &mut R) -> Self { + let mut p = [$Zq(0); $N]; + for i in p.iter_mut() { + *i = $Zq::random(distr.clone(), rng); + } + Self(p) + } + } + + impl core::ops::Add for $Rq { type Output = Self; fn add(mut self, rhs: Self) -> Self { self.0 @@ -41,7 +111,7 @@ macro_rules! impl_poly_ring { } } - impl std::ops::Add<&Self> for $Rq { + impl core::ops::Add<&Self> for $Rq { type Output = Self; fn add(mut self, rhs: &Self) -> Self { self.0 @@ -52,7 +122,7 @@ macro_rules! impl_poly_ring { } } - impl std::ops::AddAssign for $Rq { + impl core::ops::AddAssign for $Rq { fn add_assign(&mut self, rhs: Self) { self.0 .iter_mut() @@ -61,7 +131,7 @@ macro_rules! impl_poly_ring { } } - impl std::ops::AddAssign<&Self> for $Rq { + impl core::ops::AddAssign<&Self> for $Rq { fn add_assign(&mut self, rhs: &Self) { self.0 .iter_mut() @@ -70,7 +140,7 @@ macro_rules! impl_poly_ring { } } - impl std::ops::Sub for $Rq { + impl core::ops::Sub for $Rq { type Output = Self; fn sub(mut self, rhs: Self) -> Self { self.0 @@ -81,7 +151,7 @@ macro_rules! impl_poly_ring { } } - impl std::ops::Sub<&Self> for $Rq { + impl core::ops::Sub<&Self> for $Rq { type Output = Self; fn sub(mut self, rhs: &Self) -> Self { self.0 @@ -92,7 +162,7 @@ macro_rules! impl_poly_ring { } } - impl std::ops::SubAssign for $Rq { + impl core::ops::SubAssign for $Rq { fn sub_assign(&mut self, rhs: Self) { self.0 .iter_mut() @@ -101,7 +171,7 @@ macro_rules! impl_poly_ring { } } - impl std::ops::SubAssign<&Self> for $Rq { + impl core::ops::SubAssign<&Self> for $Rq { fn sub_assign(&mut self, rhs: &Self) { self.0 .iter_mut() @@ -109,14 +179,225 @@ macro_rules! impl_poly_ring { .for_each(|(x, y)| *x -= *y); } } + + impl core::ops::Add<$Zq> for $Rq { + type Output = Self; + fn add(mut self, rhs: $Zq) -> Self { + self.0[0] += rhs; + self + } + } + + impl core::ops::AddAssign<$Zq> for $Rq { + fn add_assign(&mut self, rhs: $Zq) { + self.0[0] += rhs; + } + } + + impl core::ops::Sub<$Zq> for $Rq { + type Output = Self; + fn sub(mut self, rhs: $Zq) -> Self { + self.0[0] -= rhs; + self + } + } + + impl core::ops::SubAssign<$Zq> for $Rq { + fn sub_assign(&mut self, rhs: $Zq) { + self.0[0] -= rhs; + } + } + + impl core::ops::Mul<$Zq> for $Rq { + type Output = Self; + fn mul(mut self, rhs: $Zq) -> Self { + //self.0.iter_mut().for_each(|x| *x *= rhs); + + // Barrett with re-used b + const R: $TT = ($q as $TT).next_power_of_two() as _; + let b = rhs.0 as $TT; + let b2 = (b * R / $q as $TT); + self.0.iter_mut().for_each(|a| { + a.0 = (a.0 as $TT * b - a.0 as $TT * b2 / R * $q as $TT) as $T; + if a.0 >= $q { + a.0 -= $q; + } + }); + self + } + } + + impl core::ops::Mul<&$Zq> for $Rq { + type Output = Self; + fn mul(mut self, rhs: &$Zq) -> Self { + //self.0.iter_mut().for_each(|x| *x *= *rhs); + + // Barrett with re-used b + const R: $TT = ($q as $TT).next_power_of_two() as _; + let b = rhs.0 as $TT; + let b2 = (b * R / $q as $TT); + self.0.iter_mut().for_each(|a| { + a.0 = (a.0 as $TT * b - a.0 as $TT * b2 / R * $q as $TT) as $T; + if a.0 >= $q { + a.0 -= $q; + } + }); + self + } + } + + impl core::ops::MulAssign<$Zq> for $Rq { + fn mul_assign(&mut self, rhs: $Zq) { + //self.0.iter_mut().for_each(|x| *x *= rhs); + + // Barrett with re-used b + const R: $TT = ($q as $TT).next_power_of_two() as _; + let b = rhs.0 as $TT; + let b2 = (b * R / $q as $TT); + self.0.iter_mut().for_each(|a| { + a.0 = (a.0 as $TT * b - a.0 as $TT * b2 / R * $q as $TT) as $T; + if a.0 >= $q { + a.0 -= $q; + } + }); + } + } + + impl core::ops::MulAssign<&$Zq> for $Rq { + fn mul_assign(&mut self, rhs: &$Zq) { + //self.0.iter_mut().for_each(|x| *x *= *rhs); + + // Barrett with re-used b + const R: $TT = ($q as $TT).next_power_of_two() as _; + let b = rhs.0 as $TT; + let b2 = (b * R / $q as $TT); + self.0.iter_mut().for_each(|a| { + a.0 = (a.0 as $TT * b - a.0 as $TT * b2 / R * $q as $TT) as $T; + if a.0 >= $q { + a.0 -= $q; + } + }); + } + } + + impl core::ops::Neg for $Rq { + type Output = Self; + fn neg(mut self) -> Self { + self.0.iter_mut().for_each(|x| *x = -*x); + self + } + } + + 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 $crate::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 core::ops::Index for $Rq { + type Output = $Zq; + fn index(&self, idx: usize) -> &$Zq { + &self.0[idx] + } + } + + impl core::ops::IndexMut for $Rq { + fn index_mut(&mut self, idx: usize) -> &mut $Zq { + &mut self.0[idx] + } + } + + impl $crate::ring::Norm2 for $Rq { + type Output = $TN; + type OutputSquared = $TTN; + fn norm2(&self) -> $TN { + todo!() + } + fn norm2_squared(&self) -> $TTN { + let mut ret = 0; + for i in self.0 { + ret += i.norm2_squared() as $TTN; + } + ret + } + } + + impl $crate::ring::NormInf for $Rq { + type Output = $T; + fn norm_inf(&self) -> $T { + let mut ret = 0; + for i in self.0 { + let v = i.norm_inf(); + if v > ret { + ret = v; + } + } + ret + } + } + + impl From<[$T; $N]> for $Rq { + fn from(mut p: [$T; $N]) -> Self { + for i in p.iter_mut() { + *i %= $q; + } + Self(unsafe { core::mem::transmute(p) }) + } + } + + impl core::fmt::Display for $Rq { + fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result { + write!(f, "[")?; + for i in self.0 { + write!(f, "{},", i.0)?; + } + write!(f, "]") + } + } }; } +/// Implement a modular polynomial ring over a field +/// +/// - `$Rq`: the struct that will be created, that represents a polynomial in (Z/qZ)/(x^N+1) +/// - `$N`: degree N +/// - `$Zq`: (defined) underlying modular ring Z/qZ +/// - `$T`: (defined) underlying primitive scalar type +/// - `$TT`: (defined) primitive scalar type big enough to represent the square of the biggest element of Zq +/// - `$TN`: (defined) primitive scalar type big enough to represent N times half the biggest element of Zq +/// - `$TTN`: (defined) primitive scalar type big enough to represent N times the square of half the biggest element of Zq +/// - `$q`: prime modulus q +/// - `$zeta`: 2N-th root of unity (element z of Zq such that z^m != 1 mod q for 1 <= m < 2N and z^(2N) = 1 mod q) #[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); + ( + $Rq:ident, + $N:literal, + $Zq:ident, + $T:ident, + $TT:ident, + $TN:ident, + $TTN:ident, + $q:expr, + $zeta:expr + ) => { + $crate::impl_poly_ring!($Rq, $N, $Zq, $T, $TT, $TN, $TTN, $q); impl $Rq { + #[allow(dead_code)] pub const fn zetas() -> [$T; $N] { let mut zetas = [0; $N]; let mut i: usize = 1; @@ -128,14 +409,50 @@ macro_rules! poly { } zetas } + #[allow(dead_code)] pub const fn one() -> Self { let mut x = Self([$Zq(0); $N]); x.0[0].0 = 1; x } + #[allow(dead_code)] pub fn is_one(&self) -> bool { self.0[0].0 == 1 && self.0.iter().skip(1).all(|x| x.0 == 0) } + } + + #[doc(hidden)] + impl $crate::poly::InternalNttInv for $Rq { + type Input = Self; + type Output = Self; + #[doc(hidden)] + fn __ntt_inv(mut s: Self) -> Self { + const ZETAS: [$T; $N] = $Rq::zetas(); + const R: $TT = ($q as $TT).next_power_of_two() as _; + let mut m = $N; + let mut len = 1; + while len < $N { + for start in (0..$N).step_by(2 * len) { + m -= 1; + let z = ($q - ZETAS[m]) as $TT; + let z2 = (z * R / $q as $TT); + for j in start..start + len { + let t = s.0[j]; + s.0[j] += s.0[j + len]; + let a = (t - s.0[j + len]).0 as $TT; + let mut res = (a * z - a * z2 / R * $q as $TT) as $T; + if res >= $q { + res -= $q; + } + s.0[j + len] = $Zq(res); + } + } + len *= 2; + } + const N_INV: $Zq = $Zq($N).inverse(); + s.0.iter_mut().for_each(|wj| *wj *= N_INV); + s + } /*/// Slow multiplication, for testing fn naive_mul(&self, rhs: &Self) -> Self { let mut buf = [$Zq(0); 2*$N]; @@ -148,175 +465,287 @@ macro_rules! poly { }*/ } - 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> { + type Output = $crate::ntt::NttDomain<$Rq, $Zq, $N>; + fn ntt(mut self) -> $crate::ntt::NttDomain<$Rq, $Zq, $N> { + const ZETAS: [$T; $N] = $Rq::zetas(); + const R: $TT = ($q as $TT).next_power_of_two() as _; 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]); + let z = ZETAS[m] as $TT; + let z2 = (z * R / $q as $TT); for j in start..start + len { - let t = z * self.0[j + len]; + let a = self.0[j + len].0 as $TT; + let mut res = (a * z - a * z2 / R * $q as $TT) as $T; + if res >= $q { + res -= $q; + } + let t = $Zq(res); 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 + $crate::ntt::NttDomain(self, Default::default()) } } }; } +/// Implement a modular polynomial ring with modulus 2q +/// +/// - `$Rq`: (defined) Rq +/// - `$Rb`: (defined) polynomial ring Rb with b prime such that b > N*q^2 +/// - `$R2q`: the struct that will be created +/// - `$N`: degree +/// - `$Z2q`: (defined) underlying modular ring +/// - `$T`: (defined) underlying scalar type of Z2q +/// - `$TT`: (defined) scalar type big enough to represent the square of the biggest element of Z2q +/// - `$TN`: (defined) scalar type big enough to represent N times half the biggest element of Z2q +/// - `$TTN`: (defined) scalar type big enough to represent N times the square of half the biggest element of Z2q +/// - `$TI`: TI for Z2q +/// - `$TIb`: TI for Zb +/// - `$q`: prime modulus q +/// - `$b`: prime modulus b #[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); + ( + $Rq:ident, + $Rb:ident, + $R2q:ident, + $N:literal, + $Z2q:ident, + $T:ident, + $TT:ident, + $TN:ident, + $TTN:ident, + $TI:ident, + $TIb:ident, + $q:expr, + $b:expr + ) => { + $crate::impl_poly_ring!($R2q, $N, $Z2q, $T, $TT, $TN, $TTN, 2 * $q); + impl $crate::ring::Lift<$Z2q, <$Rq as $crate::poly::Poly<$N>>::T, $T> + for <$Rq as $crate::poly::Poly<$N>>::Element + { + fn lift(self) -> $Z2q { + if self.0 > $q / 2 { + $Z2q(self.0 + $q) + } else { + $Z2q(self.0) + } + } + } + // TODO ensure we can't impl transmuting lift on rings with different T's (safety) + impl $crate::ring::Lift<$R2q, <$Rq as $crate::poly::Poly<$N>>::T, $T> for $Rq { + fn lift(self) -> $R2q { + let mut ret: $R2q = unsafe { core::mem::transmute(self) }; + ret.0.iter_mut().for_each(|x| { + //x.0 %= 2 * $q; + if x.0 > $q / 2 { + x.0 += $q; + } + }); + ret + } + /*unsafe fn lift_unchecked(self) -> $R2q { + let mut ret: $R2q = core::mem::transmute(self); + ret.0.iter_mut().for_each(|x| if x.0 > $q / 2 { + x.0 += $q; + }); + ret + }*/ + } + impl + $crate::ring::Lift< + <$Rb as $crate::poly::Poly<$N>>::Element, + $T, + <$Rb as $crate::poly::Poly<$N>>::T, + > for $Z2q + { + fn lift(self) -> <$Rb as $crate::poly::Poly<$N>>::Element { + use $crate::ring::Ring; + let x = self.center::<$TI>(); + unsafe { + if x < 0 { + <$Rb as $crate::poly::Poly<$N>>::Element::from_scalar_unchecked( + (x as $TIb + + <$Rb as $crate::poly::Poly<$N>>::Element::modulus() as $TIb) + as _, + ) + } else { + <$Rb as $crate::poly::Poly<$N>>::Element::from_scalar_unchecked(x as _) + } + } + } + } + impl $crate::ring::Lift<$Rb, $T, <$Rb as $crate::poly::Poly<$N>>::T> for &$R2q { + fn lift(self) -> $Rb { + use $crate::poly::Poly; + let mut ret: $Rb = num_traits::Zero::zero(); + for (r, s) in ret + .coefficients_mut() + .into_iter() + .zip(self.coefficients().iter()) + { + *r = s.lift(); + } + ret + } + } + impl $crate::ring::Lift<$Rb, $T, <$Rb as $crate::poly::Poly<$N>>::T> for $R2q { + fn lift(self) -> $Rb { + use $crate::poly::Poly; + let mut ret: $Rb = num_traits::Zero::zero(); + for (r, s) in ret + .coefficients_mut() + .into_iter() + .zip(self.coefficients().iter()) + { + *r = s.lift(); + } + ret + } + } + impl $crate::ring::Lift<$Z2q, <$Rb as $crate::poly::Poly<$N>>::T, $T> + for <$Rb as $crate::poly::Poly<$N>>::Element + { + fn lift(self) -> $Z2q { + use $crate::ring::Ring; + $Z2q::from_scalar( + (self.center::<$TIb>() % $Z2q::modulus() as $TIb + $Z2q::modulus() as $TIb) + as _, + ) + } + } + impl $crate::ring::Lift<$R2q, <$Rb as $crate::poly::Poly<$N>>::T, $T> for &$Rb { + fn lift(self) -> $R2q { + use $crate::poly::Poly; + let mut ret: $R2q = num_traits::Zero::zero(); + for (r, s) in ret + .coefficients_mut() + .into_iter() + .zip(self.coefficients().iter()) + { + *r = s.lift(); + } + ret + } + } + impl $crate::ring::Lift<$R2q, <$Rb as $crate::poly::Poly<$N>>::T, $T> for $Rb { + fn lift(self) -> $R2q { + use $crate::poly::Poly; + let mut ret: $R2q = num_traits::Zero::zero(); + for (r, s) in ret + .coefficients_mut() + .into_iter() + .zip(self.coefficients().iter()) + { + *r = s.lift(); + } + ret + } + } 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 - } - } }; } +/*impl + Zero, const N: usize> crate::Random for Rq where >::Element: crate::Random { + type Distribution = <>::Element as crate::Random>::Distribution; + fn random(distr: Self::Distribution, rng: &mut R) -> Self { + use rand::distributions::uniform::UniformSampler; + let mut p = Self::zero(); + for i in p.coefficients_mut().iter_mut() { + *i = distr.sample(rng); + } + p + } +}*/ + #[cfg(test)] +#[allow(unreachable_pub, unnameable_types)] mod test { use super::*; - use crate::ring; + use crate::{ + Random, + ntt::{Ntt, NttInv}, + ring, + ring::{Lift, Norm2, Ring}, + }; - use num_traits::Inv; + use num_traits::{Inv, Zero}; - 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); + ring!(Z644539222529, u64, u128, i128, 644539222529); + poly!( + R644539222529, + 256, + Z644539222529, + u64, + u128, + u64, + u128, + 644539222529, + 483489047161 + ); + type Rb = R644539222529; + + ring!(Z12289, u16, u32, i32, 12289); + poly!(R12289, 256, Z12289, u16, u32, u32, u64, 12289, 3); + ring!(Z24578, u16, u32, i32, 24578); + poly2!( + R12289, + R644539222529, + R24578, + 256, + Z24578, + u16, + u32, + u32, + u64, + i32, + i64, + 12289, + 644539222529 + ); + + ring!(Z50177, u32, u32, i64, 50177); + poly!(R50177, 256, Z50177, u32, u32, u32, u64, 50177, 66); + ring!(Z100354, u32, u64, i64, 100354); + poly2!( + R50177, + R644539222529, + R100354, + 256, + Z100354, + u32, + u64, + u32, + u64, + i64, + i64, + 50177, + 644539222529 + ); + + #[test] + fn test_add() { + let mut rng = rand::thread_rng(); + let a = R100354::random(Z100354::uniform(), &mut rng); + let b = R100354::random(Z100354::uniform(), &mut rng); + let aa: Rb = a.lift(); + let bb: Rb = b.lift(); + let c: R100354 = (aa.ntt() + bb.ntt()).ntt_inv().lift(); + assert_eq!(c, a + b); + } #[test] fn test_bitrev() { @@ -324,34 +753,346 @@ mod test { 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]); + let a = R50177::from_coefficients([Z50177::new(42); 256]); + let b = R50177::from_coefficients([Z50177::new(123); 256]); assert_eq!(a, a.clone().ntt().ntt_inv()); assert_eq!(a + b, (a.clone().ntt() + b.clone().ntt()).ntt_inv()); + assert_eq!( + a.ntt() * R50177::from_scalar(Z50177(2)).ntt(), + a.ntt() * Z50177(2) + ); } #[test] fn test_ntt_inverse() { - let a = R50177::new([Z50177::new(42); 256]).ntt(); + let a = R50177::from_coefficients([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]); + let a: Rb = R100354::from_coefficients([Z100354::new(41); 256]).lift(); + let mut b = R100354::zero(); b.0[0].0 = 2; - let c = R100354::new([Z100354::new(82); 256]); + let b: Rb = b.lift(); + let c = R100354::from_coefficients([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); + let c2: R100354 = (a.ntt() * b.ntt()).ntt_inv().lift(); + assert_eq!(c, c2); + } + + #[test] + fn test_mul_composite() { + let a = R12289::from([ + 8963, 2006, 2392, 337, 6683, 4568, 9821, 3415, 10960, 961, 9366, 9677, 1, 11799, 1932, + 723, 6931, 11767, 2196, 8786, 4833, 6507, 3251, 11263, 2816, 904, 3881, 1375, 9946, + 2597, 9292, 3569, 8407, 8971, 11646, 10336, 3562, 921, 7639, 7178, 2039, 3694, 184, + 6823, 6632, 1273, 1141, 6170, 11174, 4372, 2513, 7132, 4005, 2627, 9186, 8299, 10108, + 2887, 913, 5202, 6033, 9056, 9617, 978, 7177, 7445, 4360, 6008, 9460, 10442, 4488, + 2090, 5871, 1170, 1513, 6115, 9290, 1002, 7251, 6249, 12166, 11111, 102, 2326, 3315, + 1067, 9714, 3312, 11190, 2689, 11148, 9469, 8091, 3244, 4107, 5820, 11748, 8006, 8865, + 1201, 8810, 9925, 7638, 10318, 425, 5659, 9883, 5247, 4835, 10959, 3595, 2736, 6528, + 5004, 8242, 4798, 7218, 163, 947, 1463, 5050, 5407, 10905, 5239, 690, 7457, 7980, 8886, + 598, 6703, 2437, 1313, 11211, 8895, 9111, 10340, 3787, 1262, 10194, 5655, 7307, 5801, + 5652, 1555, 1022, 6855, 7523, 3674, 9394, 9442, 10760, 10617, 12104, 546, 3836, 6769, + 4008, 7433, 2964, 1191, 9785, 1785, 3935, 5427, 9660, 1669, 260, 6624, 4233, 450, 978, + 6346, 1306, 1782, 6631, 1630, 4202, 1037, 10656, 9664, 8468, 5014, 12046, 3470, 5805, + 5041, 5269, 899, 2409, 8550, 6487, 7600, 9571, 2653, 11755, 6642, 7370, 4990, 7921, + 1510, 3850, 10647, 669, 2924, 5495, 11344, 7006, 685, 1537, 4486, 9920, 7974, 9596, + 6775, 2203, 4387, 10214, 3964, 5776, 10414, 7604, 8602, 2199, 1229, 2641, 11250, 274, + 9644, 9343, 9449, 5165, 3220, 3140, 10844, 3356, 11583, 2603, 7571, 2241, 1362, 11873, + 195, 8165, 11254, 11398, 4675, 469, 7555, 186, 6316, 3083, 6240, 2796, 10415, 6736, + 1285, + ]); + let b = R12289::from([ + 1697, 5301, 7439, 9695, 9771, 3282, 8609, 4110, 7111, 8874, 8830, 8654, 2690, 7029, + 4604, 7833, 5373, 1045, 8569, 11228, 7396, 4880, 11678, 10935, 4340, 2624, 3172, 3224, + 10286, 254, 2651, 3755, 11288, 6468, 8645, 3432, 7492, 10186, 2381, 4906, 5904, 526, + 10816, 6145, 1843, 9388, 8456, 2087, 8026, 3093, 6335, 11790, 2834, 1830, 10249, 4951, + 10533, 5850, 776, 6326, 3747, 10339, 8470, 12051, 5853, 3513, 10404, 10495, 4283, 926, + 12042, 1002, 6502, 9696, 2447, 5745, 5494, 8044, 10610, 12262, 7033, 10338, 833, 1360, + 119, 3107, 3432, 1717, 8799, 383, 6302, 6660, 10307, 6068, 5376, 8804, 11643, 6010, + 4136, 8565, 3728, 3487, 63, 809, 5050, 8562, 943, 905, 5157, 271, 3348, 8824, 1135, + 1817, 11973, 4016, 4951, 828, 2303, 8928, 6881, 3366, 10766, 9356, 3607, 6624, 3044, + 11597, 2650, 4976, 11922, 1022, 11312, 3962, 4384, 8729, 5350, 1008, 11357, 6182, 6652, + 10880, 11822, 2552, 5265, 10456, 7233, 9560, 3435, 4362, 9632, 5725, 1607, 12189, + 10306, 5543, 10099, 4453, 9418, 9858, 3638, 8452, 7423, 3961, 3585, 1285, 10245, 6957, + 1023, 1883, 10043, 3012, 8374, 1100, 3335, 5435, 3692, 5038, 10787, 4587, 2567, 6081, + 7839, 11505, 11123, 10752, 2183, 12081, 3233, 1473, 12191, 10491, 3800, 10623, 5104, + 4256, 30, 11135, 11155, 5685, 4129, 2538, 919, 9235, 4786, 567, 2589, 10983, 11902, + 12049, 9038, 5491, 102, 12098, 8225, 6309, 849, 2071, 3242, 5116, 1781, 6524, 4409, + 7252, 8737, 10746, 11159, 1399, 1185, 7744, 2147, 8746, 9189, 3705, 6091, 3275, 9728, + 1022, 2087, 11101, 3602, 2538, 1815, 9572, 3741, 5274, 3564, 10057, 2491, 10035, 2791, + 5175, 5180, 3291, 4961, 10286, + ]); + let c = R12289::from([ + 8410, 8384, 11844, 1355, 6937, 12001, 6441, 11622, 2127, 5658, 10995, 3172, 120, 3420, + 10879, 244, 7655, 12193, 6843, 7422, 1474, 8685, 2119, 2020, 2078, 2676, 11339, 6148, + 12065, 781, 3043, 7275, 9636, 6116, 9888, 11178, 5362, 3134, 4005, 1649, 11116, 8359, + 2326, 4073, 7375, 9169, 5906, 5372, 7055, 10246, 2075, 12266, 9455, 11442, 583, 11932, + 9409, 961, 6594, 7316, 1604, 284, 3819, 6624, 1175, 8504, 860, 9791, 11793, 6498, 4167, + 11722, 11941, 351, 12269, 1215, 11194, 5044, 3535, 10450, 178, 603, 11655, 7290, 10489, + 11391, 3407, 9535, 3780, 7217, 11472, 8981, 9902, 187, 11286, 7711, 526, 3881, 8830, + 10136, 7958, 181, 7711, 12274, 10742, 501, 12209, 9748, 2300, 9874, 9873, 7956, 9478, + 6238, 5038, 8991, 1163, 721, 6916, 3434, 4110, 4687, 8470, 11276, 7141, 2039, 6143, + 5709, 6021, 218, 95, 6159, 10183, 5940, 291, 7373, 4863, 9760, 11147, 3396, 9698, 4983, + 9870, 4613, 4639, 11473, 8713, 7620, 10354, 7917, 9953, 3272, 3732, 3655, 10053, 184, + 7225, 10436, 851, 1097, 5451, 1411, 1622, 5090, 10403, 7382, 11538, 11799, 162, 7826, + 4416, 2085, 5289, 11636, 10936, 6717, 0, 9603, 7622, 639, 6668, 9382, 1898, 11232, + 10174, 5709, 1935, 2267, 11796, 10392, 2015, 1509, 379, 3218, 3754, 6631, 8536, 4547, + 3003, 4550, 3155, 1304, 6725, 6308, 12212, 10261, 9123, 1496, 2447, 5607, 4689, 6876, + 5033, 4743, 3042, 10373, 4014, 11387, 9096, 1943, 1910, 5016, 977, 3195, 2988, 9842, + 9100, 1732, 2338, 1844, 6493, 1533, 3407, 10010, 931, 3823, 11514, 7758, 11187, 2999, + 5230, 5470, 1512, 1851, 5833, 9027, 9552, 9057, 4387, 10206, 10226, 12059, 3967, 8860, + 3967, 5306, + ]); + assert_eq!((a.clone().ntt() * b.clone().ntt()).ntt_inv(), c); + + let a = R24578::from([ + 1477, 19948, 5843, 16637, 21769, 17984, 5379, 18670, 7164, 4754, 3955, 7815, 7243, + 17937, 14287, 22903, 16013, 9976, 23937, 23365, 10970, 13974, 9019, 16106, 18742, + 23541, 23760, 14111, 22822, 21569, 9513, 10424, 11341, 19484, 13789, 4674, 3199, 18257, + 1344, 17648, 9489, 6963, 21275, 18159, 24394, 7612, 3327, 10396, 7670, 146, 22066, + 15392, 22762, 22765, 675, 23926, 3390, 22017, 10652, 23925, 16349, 7913, 16667, 8924, + 19613, 3610, 3840, 6082, 6531, 19260, 19258, 23292, 7438, 19797, 22205, 14800, 1670, + 696, 3414, 14028, 4687, 23237, 19797, 19714, 10321, 15329, 7597, 15849, 10463, 18331, + 13991, 4171, 12864, 825, 7882, 9152, 8450, 15649, 22515, 7954, 21064, 2974, 18588, + 22480, 19131, 8482, 203, 19449, 209, 18942, 18292, 8251, 21451, 18580, 17504, 7536, + 4739, 13464, 7430, 22197, 15064, 22845, 24288, 1522, 12893, 13725, 9375, 692, 6698, + 22999, 1797, 4431, 17843, 4557, 20265, 1230, 19471, 20172, 23388, 4245, 15473, 19228, + 10201, 8743, 5135, 23429, 12923, 17707, 6124, 10184, 7909, 22593, 14561, 18614, 24480, + 9551, 2406, 18235, 23810, 7957, 16240, 18033, 8884, 10845, 23957, 18869, 13528, 9312, + 17078, 8895, 12545, 4775, 13367, 19183, 7016, 11046, 13573, 3776, 21869, 11845, 3247, + 19569, 2085, 18450, 11480, 18248, 19671, 5367, 5970, 13716, 19238, 12657, 12938, 9676, + 13113, 12294, 10613, 4421, 1985, 23475, 16958, 11698, 12627, 18013, 4421, 383, 14897, + 6031, 21806, 10094, 20000, 8344, 16791, 14878, 15050, 21337, 1728, 13555, 12631, 19226, + 7635, 23359, 5542, 20532, 19311, 17608, 8588, 3984, 17978, 14537, 14686, 6483, 19900, + 7513, 20518, 10473, 11312, 10161, 9263, 4283, 24222, 10432, 14265, 5471, 4052, 22753, + 11350, 6981, 17612, 8992, 6001, 10563, 22726, 18472, 8450, 11095, + ]); + let b = R24578::from([ + 715, 2634, 1346, 11801, 16559, 8048, 21303, 14436, 736, 7711, 12396, 14975, 1553, + 15937, 19039, 11001, 16175, 15678, 6262, 19101, 8093, 5016, 7715, 3993, 16328, 2836, + 21890, 11307, 18702, 1943, 8684, 13271, 21075, 6955, 9864, 5243, 17173, 21726, 1629, + 9131, 19015, 14397, 11801, 9329, 15163, 24285, 19916, 16789, 397, 11544, 18282, 6545, + 2615, 4492, 12723, 21314, 6006, 9475, 13821, 22382, 9018, 14263, 3829, 21460, 537, 205, + 24341, 4843, 14898, 22384, 8754, 20872, 14704, 16403, 215, 19995, 2161, 21449, 9236, + 1103, 16728, 103, 10486, 17074, 6231, 18069, 18585, 19617, 24339, 3684, 7943, 10883, + 13297, 19588, 12558, 22971, 15529, 5559, 17839, 23963, 571, 17560, 15751, 3212, 13224, + 17719, 4380, 6676, 23112, 3035, 19061, 10803, 2752, 19612, 17677, 10, 22441, 11368, + 6746, 14395, 6068, 7973, 3379, 15521, 12514, 4521, 10043, 16727, 14490, 8492, 23407, + 6534, 20980, 3932, 20529, 9822, 13548, 17105, 6378, 18341, 8590, 18151, 9836, 23434, + 23296, 23932, 12050, 5204, 13777, 14436, 13098, 11682, 11344, 738, 13553, 1521, 2368, + 23616, 5162, 1763, 8801, 6197, 13124, 15310, 5510, 7843, 13360, 14979, 20960, 19974, + 11913, 13697, 19376, 15259, 7945, 24478, 4311, 8005, 9599, 14779, 18339, 309, 8053, + 14194, 16718, 22447, 10590, 22382, 14239, 1514, 23525, 16386, 11845, 9979, 24219, 7936, + 23668, 5209, 4398, 15746, 19759, 13882, 7206, 9640, 7389, 20702, 2756, 14373, 2559, + 10456, 3388, 4717, 1512, 16716, 4081, 10233, 12905, 14026, 4119, 19021, 2153, 24551, + 11276, 2096, 12859, 777, 3163, 18545, 8449, 21996, 17832, 1949, 21855, 11647, 4003, + 6673, 2363, 9220, 5119, 12052, 18448, 7933, 20422, 24394, 6851, 17020, 8368, 6279, + 15055, 6261, 20943, 17379, 8506, 17057, 18002, 16142, + ]); + let c = R24578::from([ + 19851, 22117, 2203, 12013, 5450, 13847, 23681, 10098, 20401, 18766, 14538, 22302, 4993, + 17981, 18214, 13084, 12073, 1327, 5808, 18879, 10641, 18248, 21134, 20659, 8581, 3868, + 13782, 24187, 4996, 10732, 18502, 16646, 16500, 3994, 6984, 13974, 16795, 18824, 16618, + 21634, 14772, 23495, 8243, 17058, 9571, 4633, 24351, 18443, 1553, 8998, 4718, 12965, + 4490, 17418, 6276, 12912, 7797, 11824, 23974, 16964, 10014, 20121, 6766, 4841, 21515, + 4957, 3557, 1675, 20178, 12160, 2740, 18904, 15505, 13835, 19934, 8795, 23743, 21119, + 2256, 12135, 13244, 15937, 16331, 17009, 14022, 23375, 23420, 24215, 2447, 10465, 5925, + 21029, 4442, 23698, 3105, 5539, 19289, 11571, 17278, 5710, 8948, 22608, 21872, 19745, + 18627, 24259, 2690, 15282, 13564, 4139, 5779, 15209, 13862, 22063, 19152, 16855, 13022, + 19626, 2072, 8648, 6200, 15326, 13128, 16642, 24256, 7286, 21106, 3618, 5034, 12906, + 18581, 20403, 20932, 20674, 11258, 11899, 20413, 23595, 5151, 12387, 23714, 17653, + 1412, 22872, 19646, 14316, 5241, 11587, 5399, 12291, 18059, 8618, 16511, 14193, 1999, + 6022, 19639, 21320, 480, 1638, 9555, 20791, 16339, 12036, 7102, 17143, 22438, 6698, + 11231, 5860, 16882, 14933, 23205, 12317, 22449, 5019, 5178, 5567, 7201, 7496, 20129, + 416, 19736, 7684, 18099, 10348, 90, 14894, 17620, 17120, 438, 6649, 18868, 23446, + 24507, 412, 15338, 2599, 2846, 16114, 3006, 5703, 7771, 13323, 16095, 14926, 18755, + 3148, 16778, 3084, 20285, 4323, 14991, 4743, 17369, 1138, 9137, 5837, 9720, 20233, + 15588, 1392, 21955, 19691, 23005, 13373, 13563, 4971, 1607, 13098, 10897, 4026, 18082, + 16495, 2934, 17758, 19009, 18193, 13152, 5278, 15764, 6469, 3862, 670, 12749, 15749, + 3558, 12866, 17724, 1215, 11419, 20670, 344, 8418, 2514, 2862, + ]); + let aa: Rb = a.lift(); + let bb: Rb = b.lift(); + let c2: R24578 = (aa.ntt() * bb.ntt()).ntt_inv().lift(); + assert_eq!(c, c2); + + let a: Rb = R12289::from_scalar(Z12289::from_scalar(2)).lift().lift(); + let b: Rb = R12289::from_scalar(Z12289::from_scalar(2)).lift().lift(); + let c: Rb = R12289::from_scalar(Z12289::from_scalar(4)).lift().lift(); assert_eq!((a.ntt() * b.ntt()).ntt_inv(), c); + + let a: Rb = R12289::from_scalar(Z12289::from_scalar(2)).lift().lift(); + let b: Rb = R12289::from_scalar(Z12289::from_scalar(3)).lift().lift(); + let c: Rb = R12289::from_scalar(Z12289::from_scalar(6)).lift().lift(); + assert_eq!((a.ntt() * b.ntt()).ntt_inv(), c); + + let a: Rb = R12289::from_scalar(Z12289::from_scalar(2)).lift().lift(); + let b: Rb = R12289::from_scalar(Z12289::from_scalar(12288)) + .lift() + .lift(); + let c: Rb = R12289::from_scalar(Z12289::from_scalar(12287)) + .lift() + .lift(); + assert_eq!((a.ntt() * b.ntt()).ntt_inv(), c); + + let a: Rb = R24578::from([ + 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + ]) + .lift(); + let b: Rb = R24578::from([ + 2, 24576, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, + ]) + .lift(); + let c: Rb = R24578::from([ + 24577, 2, 24575, 2, 24576, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + ]) + .lift(); + assert_eq!(a.ntt() * b.ntt(), c.ntt()); + assert_eq!((a.ntt() * b.ntt()).ntt_inv(), c); + + let a = R24578::from([ + 0, 0, 24577, 0, 1, 1, 0, 1, 24577, 24577, 24577, 0, 1, 0, 0, 1, 1, 24577, 0, 1, 0, 1, + 24577, 0, 24577, 0, 1, 24577, 24577, 24577, 0, 0, 24577, 1, 0, 24577, 0, 1, 24577, 1, + 1, 0, 1, 0, 1, 0, 24577, 24577, 0, 1, 1, 24577, 0, 24577, 0, 24577, 0, 24577, 0, 1, 0, + 0, 24577, 0, 1, 0, 0, 24577, 0, 1, 24577, 24577, 0, 24577, 1, 24577, 24577, 1, 0, 1, + 24577, 1, 0, 1, 1, 24577, 1, 24577, 1, 0, 1, 0, 1, 1, 0, 1, 24577, 1, 1, 1, 24577, + 24577, 24577, 0, 0, 1, 0, 1, 24577, 24577, 1, 1, 1, 24577, 0, 1, 24577, 24577, 24577, + 0, 24577, 1, 24577, 1, 0, 0, 24577, 1, 24577, 0, 0, 0, 0, 1, 0, 24577, 24577, 24577, + 24577, 0, 0, 0, 0, 24577, 1, 0, 1, 1, 0, 0, 24577, 0, 1, 1, 0, 24577, 1, 0, 0, 24577, + 24577, 24577, 1, 1, 0, 0, 24577, 1, 24577, 0, 1, 1, 1, 0, 24577, 0, 1, 0, 24577, 1, 0, + 1, 0, 0, 0, 24577, 1, 0, 24577, 24577, 0, 24577, 1, 24577, 24577, 1, 0, 0, 0, 1, 24577, + 1, 0, 24577, 24577, 24577, 1, 0, 24577, 1, 0, 1, 1, 24577, 0, 0, 24577, 0, 0, 0, 0, 0, + 24577, 0, 0, 0, 24577, 1, 0, 1, 1, 0, 24577, 24577, 1, 1, 24577, 24577, 24577, 0, 0, 1, + 1, 0, 0, 1, 1, 1, 24577, 1, 1, 24577, 24577, 0, 24577, 24577, + ]) + .lift(); + let b = R24578::from([ + 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, + 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, + 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, + 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, + 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, + 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, + 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, + ]) + .lift(); + let c = R24578::from([ + 24574, 3, 2, 2, 24574, 24571, 24572, 24573, 24570, 1, 4, 1, 7, 24569, 24572, 2, 2, 13, + 24575, 0, 2, 6, 24575, 0, 24567, 24574, 2, 24571, 24575, 2, 24577, 2, 24577, 24576, + 24571, 3, 24577, 24574, 24575, 8, 8, 2, 24575, 24573, 5, 3, 2, 6, 4, 2, 6, 2, 6, 24568, + 1, 2, 4, 24572, 10, 7, 1, 24577, 24563, 0, 5, 7, 3, 5, 24576, 24566, 0, 24571, 24573, + 0, 24570, 3, 5, 24568, 24574, 24569, 24576, 4, 24577, 2, 2, 24568, 1, 24572, 1, 3, + 24577, 24573, 24576, 24576, 24577, 24576, 5, 24573, 15, 6, 7, 6, 24570, 0, 24575, 4, 5, + 24576, 1, 24574, 5, 24577, 7, 4, 9, 0, 7, 24575, 6, 3, 24576, 6, 24575, 13, 24572, + 24575, 2, 24569, 3, 6, 1, 4, 4, 5, 24577, 5, 24576, 24577, 9, 24576, 24577, 24576, + 24574, 3, 2, 24571, 6, 24572, 5, 8, 24575, 24577, 24574, 24577, 8, 24577, 24577, 24575, + 8, 3, 24571, 24574, 24570, 3, 24574, 24563, 24577, 24576, 0, 0, 8, 24573, 24577, 24576, + 1, 24576, 4, 6, 10, 4, 24571, 6, 1, 24574, 24574, 1, 7, 4, 0, 24570, 6, 24573, 0, 12, + 0, 24577, 5, 4, 2, 24575, 24575, 2, 24577, 24576, 24574, 4, 3, 24576, 24575, 24577, 4, + 2, 4, 2, 3, 2, 11, 24571, 3, 24576, 24575, 24576, 24575, 24575, 24568, 24576, 3, 4, 7, + 0, 24575, 5, 24570, 24573, 24573, 5, 2, 24577, 24561, 24564, 2, 1, 24573, 0, 3, 3, 5, + 2, 7, 24577, 24574, 14, 24573, 4, 24567, 24573, + ]); + let c2: R24578 = (a.ntt() * b.ntt()).ntt_inv().lift(); + assert_eq!(c, c2); + + let mut rng = rand::thread_rng(); + let a: Rb = R100354::random(Z100354::uniform(), &mut rng).lift(); + let b: Rb = R100354::random(Z100354::uniform(), &mut rng).lift(); + let c: Rb = R100354::random(Z100354::uniform(), &mut rng).lift(); + assert_eq!((a.clone().ntt() * Rb::one().ntt()).ntt_inv(), a); + assert_eq!( + (a.clone().ntt() + b.clone().ntt()) * c.clone().ntt(), + a.clone().ntt() * c.clone().ntt() + b.clone().ntt() * c.clone().ntt() + ); + assert_eq!( + ((a.clone().ntt() + b.clone().ntt()) * c.clone().ntt()).ntt_inv(), + (a.clone().ntt() * c.clone().ntt() + b.clone().ntt() * c.clone().ntt()).ntt_inv() + ); + assert_eq!( + ((a.clone() + b.clone()).ntt() * c.clone().ntt()).ntt_inv(), + (a.clone().ntt() * c.clone().ntt() + b.clone().ntt() * c.clone().ntt()).ntt_inv() + ); + assert_eq!( + ((a.clone().ntt() + b.clone().ntt()) * c.clone().ntt()).ntt_inv(), + (a.clone().ntt() * c.clone().ntt()).ntt_inv() + + (b.clone().ntt() * c.clone().ntt()).ntt_inv() + ); + assert_eq!( + ((a.clone() + &b).ntt() * c.clone().ntt()).ntt_inv(), + (a.clone().ntt() * c.clone().ntt()).ntt_inv() + + (b.clone().ntt() * c.clone().ntt()).ntt_inv() + ); + } + + #[test] + fn test_lift() { + let mut a = R50177::from_coefficients([Z50177::new(1); 256]); + a.0[1].0 = 50176; + let mut b = R100354::from_coefficients([Z100354::new(1); 256]); + b.0[1].0 = 100353; + assert_eq!(a.clone().lift(), b); + //assert_eq!(unsafe {a.lift_unchecked()}, b); + } + + #[test] + fn test_lift_big() { + let mut rng = rand::thread_rng(); + let a = R100354::random(Z100354::uniform(), &mut rng); + let aa: Rb = a.lift(); + let a2: R100354 = aa.lift(); + assert_eq!(a, a2); + } + + #[test] + fn test_norm2() { + assert_eq!( + R50177::from_coefficients([Z50177::new(1); 256]).norm2_squared(), + 256 + ); + assert_eq!( + R50177::from_coefficients([Z50177::new(50176); 256]).norm2_squared(), + 256 + ); + assert_eq!( + R100354::from_coefficients([Z100354::new(1); 256]).norm2_squared(), + 256 + ); + assert_eq!( + R100354::from_coefficients([Z100354::new(100353); 256]).norm2_squared(), + 256 + ); } } diff --git a/src/ring.rs b/src/ring.rs index dc2b5cb..5aa0662 100644 --- a/src/ring.rs +++ b/src/ring.rs @@ -1,28 +1,162 @@ -pub trait Ring { - type T; +//! Macros and traits for modular rings Z/qZ. + +use num_traits::cast::AsPrimitive; + +/// An element of ring Z/qZ represented by type `T` +pub trait Ring: AsRef + AsMut + Sized { + /// q + fn modulus() -> T; + + /// q/2 rounded down + fn modulus_half() -> T; + + /// New element from a scalar, safely reducing modulo q + fn from_scalar(v: T) -> Self; + + /// New element from a scalar, without reduction modulo q + /// + /// # Safety + /// If v >= q, following operations may produce wrong results or panic. + unsafe fn from_scalar_unchecked(v: T) -> Self; + + /// New element from a scalar, assuming 0 <= v < 2q + fn from_scalar_small_reduce(v: T) -> Self; + + /// positive representative in [0..q-1] + fn to_scalar(self) -> T; + + /// Sample a uniform element in the ring + fn sample_uniform(rng: impl rand::Rng) -> Self; + + /// Uniform distribution in the ring + fn uniform() -> UniformZq; + + /// Uniform distribution in a ball of maximum infinity norm `max` + fn uniform_ball(max_norm: T) -> UniformZq; + + /// Uniform distribution in a positive semiball of maximum infinity norm `max` + fn uniform_positive_semiball(max_norm: T) -> UniformZq; + + /// Centered representative in [-q/2, q/2] + fn center(self) -> TI + where + Self: Copy, + T: AsPrimitive + Copy + Ord, + TI: 'static + Copy + core::ops::Sub, + { + if self.to_scalar() <= Self::modulus_half() { + self.to_scalar().as_() + } else { + self.to_scalar().as_() - Self::modulus().as_() + } + } } +/// Implement a modular ring +/// +/// - `$Zq`: the struct that will be created +/// - `$T`: (defined) scalar type big enough to represent 2q +/// - `$TT`: (defined) scalar type big enough to represent (q-1)^2 +/// - `$TI`: (defined) signed scalar type big enough to represent (q-1)^2 +/// - `$q`: prime modulus +/// +/// ``` +/// gwrizienn::ring!(Z50177, u32, u64, i64, 50177); +/// assert_eq!(Z50177(50176) + Z50177(43), Z50177(42)); +/// assert_eq!(Z50177(123) * Z50177(123).inverse(), Z50177(1)); +/// ``` #[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)] + #[repr(transparent)] pub struct $Zq(pub $T); - impl $crate::ring::Ring for $Zq { - type T = $T; + impl AsRef<$T> for $Zq { + fn as_ref(&self) -> &$T { + &self.0 + } + } + + impl AsMut<$T> for $Zq { + fn as_mut(&mut self) -> &mut $T { + &mut self.0 + } + } + + impl $crate::ring::Ring<$T> for $Zq { + fn modulus() -> $T { + $q + } + + fn modulus_half() -> $T { + $q / 2 + } + + fn from_scalar(v: $T) -> Self { + Self(v % $q) + } + + fn from_scalar_small_reduce(v: $T) -> Self { + Self(Self::small_reduce(v)) + } + + fn to_scalar(self) -> $T { + self.0 + } + + unsafe fn from_scalar_unchecked(v: $T) -> Self { + Self(v) + } + + //#[cfg(feature = "rand")] + //#[allow(dead_code)] + fn sample_uniform(mut rng: impl rand::Rng) -> Self { + Self(rng.gen_range(0..$q)) + } + + //#[cfg(feature = "rand")] + //#[allow(dead_code)] + fn uniform() -> $crate::ring::UniformZq<$T, Self> { + use rand::distributions::uniform::UniformSampler; + $crate::ring::UniformZq::new_inclusive(Self(0), Self($q - 1)) + } + + /// Uniform distribution in a ball of maximum infinity norm `max` + //#[cfg(feature = "rand")] + //#[allow(dead_code)] + fn uniform_ball(max_norm: $T) -> $crate::ring::UniformZq<$T, Self> { + use rand::distributions::uniform::UniformSampler; + $crate::ring::UniformZq::new_inclusive(Self($q - max_norm), Self($q + max_norm)) + } + + /// Uniform distribution in a positive semiball of maximum infinity norm `max` + //#[cfg(feature = "rand")] + //#[allow(dead_code)] + fn uniform_positive_semiball(max_norm: $T) -> $crate::ring::UniformZq<$T, Self> { + use rand::distributions::uniform::UniformSampler; + $crate::ring::UniformZq::new_inclusive(Self(0), Self(max_norm)) + } } impl $Zq { - pub const fn modulus() -> $T { - $q - } + /// Create an element + #[allow(dead_code)] pub const fn new(v: $T) -> Self { Self(v % $q) } + + /// Reduce v modulo q, assuming that 0 <= v < 2q + #[allow(dead_code)] const fn small_reduce(v: $T) -> $T { if v < $q { v } else { v - $q } } + + /// Compute the inverse modulo q + /// + /// May output garbage if `self` is not invertible. + #[allow(dead_code)] const fn inverse(self) -> Self { let (mut ro, mut r) = ($q, self.0); let (mut to, mut t) = (0, 1); @@ -38,36 +172,48 @@ macro_rules! ring { } } - impl std::ops::Add for $Zq { + impl num_traits::AsPrimitive<$T> for $Zq { + fn as_(self) -> $T { + self.0 + } + } + + impl num_traits::AsPrimitive<$Zq> for $T { + fn as_(self) -> $Zq { + $Zq(self) + } + } + + impl core::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 { + impl core::ops::AddAssign for $Zq { fn add_assign(&mut self, rhs: Self) { *self = *self + rhs; } } - impl std::ops::Sub for $Zq { + impl core::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 { + impl core::ops::SubAssign for $Zq { fn sub_assign(&mut self, rhs: Self) { *self = *self - rhs; } } - impl std::ops::Neg for $Zq { + impl core::ops::Neg for $Zq { type Output = Self; fn neg(self) -> Self { - Self(Self::small_reduce($q - self.0)) + if self.0 == 0 { self } else { Self($q - self.0) } } } @@ -80,15 +226,26 @@ macro_rules! ring { } } - impl std::ops::Mul for $Zq { + impl core::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) + + // Barrett + // let a = self.0 as $TT; + // let b = rhs.0 as $TT; + // const r: $TT = ($q as $TT).next_power_of_two() as _; + // let res = (a * b - a * (b * r / $q as $TT) / r * $q as $TT) as $T; + // Self(if res >= $q { + // res - $q + // } else { + // res + // }) } } - impl std::ops::MulAssign for $Zq { + impl core::ops::MulAssign for $Zq { fn mul_assign(&mut self, rhs: Self) { *self = *self * rhs; } @@ -109,12 +266,168 @@ macro_rules! ring { self.inverse() } } + + // TODO how to check our feature from where the macro is called? + //#[cfg(feature = "rand")] + impl rand::distributions::uniform::SampleUniform for $Zq { + type Sampler = $crate::ring::UniformZq<$T, $Zq>; + } + + //#[cfg(feature = "rand")] + impl> $crate::Random for $Zq { + fn random(distr: D, rng: &mut R) -> Self { + distr.sample(rng) + } + } + impl> $crate::Random for $Zq { + fn random(distr: D, rng: &mut R) -> Self { + let r = distr.sample(rng); + if r < 0 { + Self(($q + r) as $T) + } else { + Self(r as $T) + } + } + } + + impl $crate::ring::Norm2 for $Zq { + type Output = $T; + type OutputSquared = $TT; + fn norm2(&self) -> $T { + use $crate::ring::Ring; + self.center::<$TI>().unsigned_abs() as $T + } + fn norm2_squared(&self) -> $TT { + use $crate::ring::Ring; + let v: $TI = self.center::<$TI>(); + (v * v) as $TT + } + } + + impl $crate::ring::NormInf for $Zq { + type Output = $T; + fn norm_inf(&self) -> $T { + use $crate::ring::Ring; + self.center::<$TI>().abs() as _ + } + } + + impl core::fmt::Display for $Zq { + fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result { + self.0.fmt(f) + } + } }; } -crate::ring!(Z2, u8, u8, i8, 2); +/// Uniform distribution over a ring +#[cfg(feature = "rand")] +#[derive(Clone, Copy, Debug)] +pub struct UniformZq( + pub(crate) rand::distributions::uniform::UniformInt, + pub(crate) core::marker::PhantomData, +); + +#[cfg(feature = "rand")] +impl> + rand::distributions::uniform::UniformSampler for UniformZq +where + rand::distributions::uniform::UniformInt: + rand::distributions::uniform::UniformSampler, +{ + type X = Zq; + fn new(low: B1, high: B2) -> Self + where + B1: rand::distributions::uniform::SampleBorrow + Sized, + B2: rand::distributions::uniform::SampleBorrow + Sized, + { + UniformZq( + rand::distributions::uniform::UniformInt::::new( + low.borrow().to_scalar(), + high.borrow().to_scalar(), + ), + Default::default(), + ) + } + fn new_inclusive(low: B1, high: B2) -> Self + where + B1: rand::distributions::uniform::SampleBorrow + Sized, + B2: rand::distributions::uniform::SampleBorrow + Sized, + { + UniformZq( + rand::distributions::uniform::UniformInt::::new_inclusive( + low.borrow().to_scalar(), + high.borrow().to_scalar(), + ), + Default::default(), + ) + } + fn sample(&self, rng: &mut R) -> Self::X { + Zq::from_scalar_small_reduce(self.0.sample(rng)) + } +} + +impl rand::distributions::Distribution for UniformZq +where + Self: rand::distributions::uniform::UniformSampler, + Zq: rand::distributions::uniform::SampleUniform, +{ + fn sample(&self, rng: &mut R) -> Zq { + rand::distributions::uniform::UniformSampler::sample(self, rng) + } +} + +//crate::ring!(Z2, u8, u8, i8, 2); + +/// Lifting between types +pub trait Lift { + /// Lift to another modular ring + /// + /// Injective only if the destination ring is bigger. + fn lift(self) -> To; + //unsafe fn lift_unchecked(self) -> To; +} + +/*impl, To: Ring> Lift for Fr { + fn lift(self) -> To { + To::from_scalar(self.to_scalar()) + } + unsafe fn lift_unchecked(self) -> To { + To::from_scalar_unchecked(self.to_scalar()) + } +}*/ + +/*impl, U, Fr: Ring, To: Ring> Lift for Fr { + fn lift(self) -> To { + To::from_scalar(self.to_scalar().into()) + } + /*unsafe fn lift_unchecked(self) -> To { + To::from_scalar_unchecked(self.to_scalar().into()) + }*/ +}*/ + +/// Euclidean norm +pub trait Norm2 { + /// Output type + type Output; + /// Output type + type OutputSquared; + /// Euclidean norm + fn norm2(&self) -> Self::Output; + /// Square of the euclidean norm (may be faster than the norm, depending on implementation) + fn norm2_squared(&self) -> Self::OutputSquared; +} + +/// Infinity norm +pub trait NormInf { + /// Output type + type Output; + /// Infinity norm + fn norm_inf(&self) -> Self::Output; +} #[cfg(test)] +#[allow(unreachable_pub, unnameable_types)] mod test { use super::*; @@ -129,5 +442,10 @@ mod test { assert_eq!(Z50177::new(1), Z50177::new(50178)); assert_eq!(Z50177::new(123) * Z50177::new(123).inv(), Z50177::one()); + + assert_eq!(Z50177::new(1).norm2_squared(), 1); + assert_eq!(Z50177::new(50176).norm2_squared(), 1); + assert_eq!(Z50177::new(2).norm2_squared(), 4); + assert_eq!(Z50177::new(50175).norm2_squared(), 4); } } diff --git a/src/vector.rs b/src/vector.rs index 4b0ce28..dbe6144 100644 --- a/src/vector.rs +++ b/src/vector.rs @@ -1,22 +1,82 @@ -use std::{ +//! Vector types that can contain modular scalars (`Ring`) or polynomials (`Poly`). +//! Matrices must not be implemented as vectors of vectors; use the `matrix` module instead. + +//use crate::ntt::Fma; + +use core::{ mem::MaybeUninit, - ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign}, + ops::{Add, AddAssign, Index, IndexMut, Mul, MulAssign, Sub, SubAssign}, }; +use num_traits::Zero; -use num_traits::{One, Zero}; +/// Vector of E +pub trait VectorTrait: Index + IndexMut {} +impl VectorTrait for Vector {} + +/// Vector of dimension N #[derive(Clone, Debug, Eq, Hash, PartialEq)] pub struct Vector(pub [E; N]); +/// Reference to a vector +#[derive(Debug, Eq, Hash, PartialEq)] +pub struct VectorRef<'a, E, const N: usize>(pub &'a [E; N]); + +/// Mutable reference to a vector +#[derive(Debug, Eq, Hash, PartialEq)] +pub struct VectorMut<'a, E, const N: usize>(pub &'a mut [E; N]); + +impl<'a, E, const N: usize> Clone for VectorRef<'a, E, N> { + fn clone(&self) -> Self { + *self + } +} + +impl<'a, E, const N: usize> Copy for VectorRef<'a, E, N> {} + /*impl Vector { pub fn norm_inf(&self) -> Zq { Zq(self.0.iter().map(|i| i.norm_inf().0).max().unwrap()) } }*/ +impl Vector { + /// Concatenate two vectors + pub fn concat(self, other: Vector) -> Vector { + assert_eq!(M + N, MN); + let mut z = MaybeUninit::<[E; MN]>::uninit(); + unsafe { + for (new, old) in z + .assume_init_mut() + .iter_mut() + .zip(self.0.into_iter().chain(other.0.into_iter())) + { + *new = old; + } + } + Vector(unsafe { z.assume_init() }) + } + /// Get a reference to a subvector + pub fn get_sub<'a, const M: usize>(&'a self, idx: usize) -> VectorRef<'a, E, M> + where + &'a [E]: TryInto<&'a [E; M]>, + { + VectorRef(if let Ok(sub) = &self.0[idx..idx + M].try_into() { + sub + } else { + panic!("Subvector out of bounds") + }) + } + /// Get a reference to the vector + pub fn get_ref(&self) -> VectorRef<'_, E, N> { + VectorRef(&self.0) + } +} + impl Zero for Vector where - E: AddAssign + Zero, + E: Zero, + Self: Add, { fn zero() -> Self { let mut z = MaybeUninit::<[E; N]>::uninit(); @@ -57,6 +117,53 @@ where } } +impl<'a, E, const N: usize> Add> for Vector +where + E: 'a + AddAssign<&'a E>, +{ + type Output = Self; + fn add(mut self, rhs: VectorRef<'a, E, N>) -> Self { + self.0 + .iter_mut() + .zip(rhs.0.iter()) + .for_each(|(x, y)| *x += y); + self + } +} + +impl AddAssign for Vector +where + E: AddAssign, +{ + fn add_assign(&mut self, rhs: Self) { + self.0.iter_mut().zip(rhs.0).for_each(|(x, y)| *x += y); + } +} + +impl<'a, E, const N: usize> AddAssign<&'a Self> for Vector +where + E: 'a + AddAssign<&'a E>, +{ + fn add_assign(&mut self, rhs: &'a Self) { + self.0 + .iter_mut() + .zip(rhs.0.iter()) + .for_each(|(x, y)| *x += y); + } +} + +impl<'a, E, const N: usize> AddAssign> for Vector +where + E: 'a + AddAssign<&'a E>, +{ + fn add_assign(&mut self, rhs: VectorRef<'a, E, N>) { + self.0 + .iter_mut() + .zip(rhs.0.iter()) + .for_each(|(x, y)| *x += y); + } +} + impl Sub for Vector where E: SubAssign, @@ -68,6 +175,15 @@ where } } +impl SubAssign for Vector +where + E: SubAssign, +{ + fn sub_assign(&mut self, rhs: Self) { + self.0.iter_mut().zip(rhs.0).for_each(|(x, y)| *x -= y); + } +} + impl<'a, E, const N: usize> Sub<&'a Self> for Vector where E: 'a + SubAssign<&'a E>, @@ -82,26 +198,98 @@ where } } -/// Vector-scalar product -impl<'a, E, const N: usize> Mul<&'a E> for Vector +impl<'a, E, const N: usize> SubAssign<&'a Self> for Vector where - E: 'a + Clone + std::fmt::Debug + Zero + MulAssign<&'a E>, + E: 'a + SubAssign<&'a E>, +{ + fn sub_assign(&mut self, rhs: &'a Self) { + self.0 + .iter_mut() + .zip(rhs.0.iter()) + .for_each(|(x, y)| *x -= y); + } +} + +impl<'a, E, const N: usize> Sub> for Vector +where + E: 'a + SubAssign<&'a E>, +{ + type Output = Self; + fn sub(mut self, rhs: VectorRef<'a, E, N>) -> Self { + self.0 + .iter_mut() + .zip(rhs.0.iter()) + .for_each(|(x, y)| *x -= y); + self + } +} + +impl<'a, E, const N: usize> SubAssign> for Vector +where + E: 'a + SubAssign<&'a E>, +{ + fn sub_assign(&mut self, rhs: VectorRef<'a, E, N>) { + self.0 + .iter_mut() + .zip(rhs.0.iter()) + .for_each(|(x, y)| *x -= y); + } +} + +/// Vector-scalar product +impl Mul for Vector +where + E: MulAssign, { type Output = Vector; - fn mul(mut self, rhs: &'a E) -> Vector { + fn mul(mut self, rhs: F) -> 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 +/// Vector-scalar product +impl MulAssign for Vector where - E: 'a + Clone + Zero + MulAssign<&'a E> + Add, + E: MulAssign, +{ + fn mul_assign(&mut self, rhs: F) { + self.0.iter_mut().for_each(|s| *s *= rhs); + } +} + +/// Inner product +pub trait Dot { + /// Result type (should be the vector's element type) + type Output; + /// Inner product + fn dot(self, rhs: Rhs) -> Self::Output; +} + +impl Dot for Vector +where + E: Zero + MulAssign + Add, { type Output = E; - fn mul(self, rhs: &'a Self) -> E { + // TODO optimize (reduce and copy less often) + fn dot(self, rhs: Self) -> E { + self.0 + .into_iter() + .zip(rhs.0) + .fold(E::zero(), |e, (mut x, y)| { + x *= y; + e + x + }) + } +} + +impl<'a, E, const N: usize> Dot<&'a Self> for Vector +where + E: 'a + Zero + MulAssign<&'a E> + Add, +{ + type Output = E; + // TODO optimize (reduce and copy less often) + fn dot(self, rhs: &'a Self) -> E { self.0 .into_iter() .zip(rhs.0.iter()) @@ -112,74 +300,222 @@ where } } -/// Inner product -// TODO optimize (reduce and copy less often) -impl<'a, E, const N: usize> Mul for &'a Vector +impl<'a, E, const N: usize> Dot> for Vector where - E: 'a + Clone + Zero + Mul<&'a E, Output = E> + Add, + E: 'a + Zero + MulAssign<&'a E> + Add, { type Output = E; - fn mul(self, rhs: Self) -> E { + // TODO optimize (reduce and copy less often) + fn dot(self, rhs: VectorRef<'a, E, N>) -> E { + self.0 + .into_iter() + .zip(rhs.0.iter()) + .fold(E::zero(), |e, (mut x, y)| { + x *= y; + e + x + }) + } +} + +impl<'a, E, const N: usize> Dot for &'a Vector +where + E: 'a + Clone + Mul<&'a E, Output = E> + Add + Zero, + //E: Fma<&'a E, &'a E> +{ + type Output = E; + fn dot(self, rhs: Self) -> E { self.0 .iter() .zip(rhs.0.iter()) .fold(E::zero(), |e, (x, y)| e + x.clone() * y) + //.fold(E::zero(), |mut e, (x, y)| { + // e.fma(x, y); + // e + //}) } } -/// Matrix-vector product -impl<'a, E, const M: usize, const N: usize> Mul<&'a Vector> for Vector, N> +impl<'a, E, const N: usize> Dot> for &'a Vector where - E: 'a + Clone + std::fmt::Debug + Zero + MulAssign<&'a E> + Add + AddAssign, + E: 'a + Clone + Mul<&'a E, Output = E> + Add + Zero, + //E: Fma<&'a E, &'a E> { - 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 + type Output = E; + fn dot(self, rhs: VectorRef<'a, E, N>) -> E { + self.0 + .iter() + .zip(rhs.0.iter()) + .fold(E::zero(), |e, (x, y)| e + x.clone() * y) + //.fold(E::zero(), |mut e, (x, y)| { + // e.fma(x, y); + // e + //}) } } -/// 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 Index for Vector { + type Output = E; + fn index(&self, idx: usize) -> &E { + &self.0[idx] } } -impl Vector, N> +impl IndexMut for Vector { + fn index_mut(&mut self, idx: usize) -> &mut E { + &mut self.0[idx] + } +} + +#[cfg(feature = "rand")] +impl crate::Random for Vector where - E: Zero + One, - Vector: AddAssign + Zero, + D: Clone + rand::distributions::Distribution, + E: crate::Random, { - pub fn id() -> Self { - let mut id = Self::zero(); - for (i, a) in id.0.iter_mut().enumerate() { - a.0[i] = E::one(); + fn random(distr: D, rng: &mut R) -> Self { + let mut z = MaybeUninit::<[E; N]>::uninit(); + unsafe { + for i in z.assume_init_mut() { + *i = E::random(distr.clone(), rng); + } } - id + Self(unsafe { z.assume_init() }) + } +} + +impl crate::ntt::Ntt for Vector +where + E: crate::ntt::Ntt, +{ + type Output = Vector; + fn ntt(self) -> Vector { + let mut z = MaybeUninit::<[O; N]>::uninit(); + unsafe { + for (i, j) in z.assume_init_mut().iter_mut().zip(self.0) { + *i = j.ntt(); + } + } + Vector(unsafe { z.assume_init() }) + } +} + +impl crate::ntt::NttInv for Vector +where + E: crate::ntt::NttInv, +{ + type Output = Vector; + fn ntt_inv(self) -> Vector { + let mut z = MaybeUninit::<[O; N]>::uninit(); + unsafe { + for (i, j) in z.assume_init_mut().iter_mut().zip(self.0) { + *i = j.ntt_inv(); + } + } + Vector(unsafe { z.assume_init() }) + } +} + +impl crate::ring::Lift, T, U> for Vector +where + Fr: crate::ring::Lift, +{ + fn lift(self) -> Vector { + let mut z = MaybeUninit::<[To; N]>::uninit(); + unsafe { + for (i, j) in z.assume_init_mut().iter_mut().zip(self.0) { + *i = j.lift(); + } + } + Vector(unsafe { z.assume_init() }) + } + /*unsafe fn lift_unchecked(self) -> $R2q { + let mut ret: $R2q = core::mem::transmute(self); + ret.0.iter_mut().for_each(|x| if x.0 > $q / 2 { + x.0 += $q; + }); + ret + }*/ +} + +impl crate::ring::Norm2 for Vector +where + TT: Zero + AddAssign, + E: crate::ring::Norm2, +{ + type Output = T; + type OutputSquared = TT; + fn norm2(&self) -> T { + todo!() + } + fn norm2_squared(&self) -> TT { + let mut ret = TT::zero(); + for i in self.0.iter() { + ret += i.norm2_squared(); + } + ret + } +} + +impl crate::ring::NormInf for Vector +where + T: Ord, + E: crate::ring::NormInf, +{ + type Output = T; + fn norm_inf(&self) -> T { + self.0 + .iter() + .map(crate::ring::NormInf::norm_inf) + .max() + .unwrap() + } +} + +impl core::ops::Neg for Vector +where + Self: Zero + Sub, +{ + type Output = Self; + fn neg(self) -> Self { + Self::zero() - self + } +} + +#[cfg(feature = "zeroize")] +impl zeroize::Zeroize for Vector { + fn zeroize(&mut self) { + self.0.zeroize(); + } +} + +impl core::fmt::Display for Vector { + fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result { + write!(f, "[")?; + for i in self.0.iter() { + write!(f, "{},", i)?; + } + write!(f, "]") + } +} + +impl<'a, E, const N: usize> From> for Vector +where + &'a [E; N]: Into<[E; N]>, +{ + fn from(s: VectorRef<'a, E, N>) -> Vector { + Vector(s.0.into()) } } #[cfg(test)] +#[allow(unreachable_pub, unnameable_types)] 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); + ring!(Z50177, u32, u32, i64, 50177); + poly!(R50177, 256, Z50177, u32, u32, u32, u64, 50177, 66); #[test] fn test_basic() { diff --git a/test.sage b/test.sage deleted file mode 100644 index fe47d9a..0000000 --- a/test.sage +++ /dev/null @@ -1,39 +0,0 @@ -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)