diff --git a/Cargo.toml b/Cargo.toml
index 4ac6218..fb6e8e1 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -7,12 +7,13 @@ autobenches = false
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
 [dependencies]
-qfall-math = { git = "https://github.com/qfall/math", rev="5f50c9cd31c869462d959774fb4b51fcd1727dbe" }
+qfall-math = { git = "https://github.com/qfall/math", branch = "dev" }
+flint-sys = "0.7.3"
 sha2 = "0.10.6"
 serde = {version="1.0", features=["derive"]}
 serde_json = "1.0"
 typetag = "0.2"
-criterion = { version = "0.7", features = ["html_reports"] }
+criterion = { version = "0.8", features = ["html_reports"] }
 
 [profile.bench]
 debug = true
diff --git a/src/primitive/psf/gpv.rs b/src/primitive/psf/gpv.rs
index 74fba71..2f22162 100644
--- a/src/primitive/psf/gpv.rs
+++ b/src/primitive/psf/gpv.rs
@@ -112,7 +112,7 @@ impl PSF for PSFGPV {
     /// ```
     fn samp_d(&self) -> MatZ {
         let m = &self.gp.n * &self.gp.k + &self.gp.m_bar;
-        MatZ::sample_d_common(&m, &self.gp.n, &self.s).unwrap()
+        MatZ::sample_discrete_gauss(&m, 1, 0, &self.s).unwrap()
     }
 
     /// Samples an `e` in the domain using SampleD with a short basis that is generated
@@ -157,14 +157,7 @@ impl PSF for PSFGPV {
 
         let center = MatQ::from(&(-1 * &sol));
 
-        sol + MatZ::sample_d_precomputed_gso(
-            short_base,
-            short_base_gso,
-            &self.gp.n,
-            &center,
-            &self.s,
-        )
-        .unwrap()
+        sol + MatZ::sample_d_precomputed_gso(short_base, short_base_gso, &center, &self.s).unwrap()
     }
 
     /// Implements the efficiently computable function `f_a` which here corresponds to
diff --git a/src/primitive/psf/gpv_ring.rs b/src/primitive/psf/gpv_ring.rs
index 435ada4..60a77bc 100644
--- a/src/primitive/psf/gpv_ring.rs
+++ b/src/primitive/psf/gpv_ring.rs
@@ -117,7 +117,7 @@ impl PSF for PSFGPVRing {
     /// ```
     fn samp_d(&self) -> MatPolyOverZ {
         let dimension = self.gp.modulus.get_degree() * (&self.gp.k + 2);
-        let sample = MatZ::sample_d_common(dimension, &self.gp.n, &self.s).unwrap();
+        let sample = MatZ::sample_discrete_gauss(dimension, 1, 0, &self.s).unwrap();
         MatPolyOverZ::from_coefficient_embedding((&sample, self.gp.modulus.get_degree() - 1))
     }
 
@@ -205,7 +205,6 @@ impl PSF for PSFGPVRing {
             + MatPolyOverZ::sample_d(
                 &short_basis,
                 self.gp.modulus.get_degree(),
-                &self.gp.n,
                 &center_embedded,
                 &self.s,
             )
diff --git a/src/primitive/psf/mp_perturbation.rs b/src/primitive/psf/mp_perturbation.rs
index b902309..c5c65e4 100644
--- a/src/primitive/psf/mp_perturbation.rs
+++ b/src/primitive/psf/mp_perturbation.rs
@@ -186,14 +186,8 @@ pub(crate) fn randomized_nearest_plane_gadget(
 
     // just as PSFGPV::samp_p
     long_solution
-        + MatZ::sample_d_precomputed_gso(
-            short_basis_gadget,
-            short_basis_gadget_gso,
-            &psf.gp.n,
-            &center,
-            &s,
-        )
-        .unwrap()
+        + MatZ::sample_d_precomputed_gso(short_basis_gadget, short_basis_gadget_gso, &center, &s)
+            .unwrap()
 }
 
 impl PSF for PSFPerturbation {
@@ -269,7 +263,7 @@ impl PSF for PSFPerturbation {
     /// ```
     fn samp_d(&self) -> MatZ {
         let m = &self.gp.n * &self.gp.k + &self.gp.m_bar;
-        MatZ::sample_d_common(&m, &self.gp.n, &self.s * &self.r).unwrap()
+        MatZ::sample_discrete_gauss(m, 1, 0, &self.s * &self.r).unwrap()
     }
 
     /// Samples an `e` in the domain using SampleD that is generated
@@ -318,8 +312,7 @@ impl PSF for PSFPerturbation {
         vec_u: &MatZq,
     ) -> MatZ {
         // Sample perturbation p <- D_{ZZ^m, r * √Σ_2}
-        let vec_p =
-            MatZ::sample_d_common_non_spherical(&self.gp.n, mat_sqrt_sigma_2, &self.r).unwrap();
+        let vec_p = MatZ::sample_d_common_non_spherical(mat_sqrt_sigma_2, &self.r).unwrap();
 
         // v = u - A * p
         let vec_v = vec_u - mat_a * &vec_p;
diff --git a/src/sample/g_trapdoor/trapdoor_distribution.rs b/src/sample/g_trapdoor/trapdoor_distribution.rs
index 039022b..4c73539 100644
--- a/src/sample/g_trapdoor/trapdoor_distribution.rs
+++ b/src/sample/g_trapdoor/trapdoor_distribution.rs
@@ -114,7 +114,7 @@ impl TrapdoorDistributionRing for SampleZ {
         let nr_cols = i64::try_from(nr_cols).unwrap();
         let mut out_mat = MatPolyOverZ::new(1, nr_cols);
         for j in 0..nr_cols {
-            let sample = PolyOverZ::sample_discrete_gauss(n - 1, n, 0, s).unwrap();
+            let sample = PolyOverZ::sample_discrete_gauss(n - 1, 0, s).unwrap();
             out_mat.set_entry(0, j, &sample).unwrap();
         }
 
diff --git a/src/utils.rs b/src/utils.rs
index 8df1a87..4177430 100644
--- a/src/utils.rs
+++ b/src/utils.rs
@@ -12,4 +12,5 @@
 
 pub mod common_encodings;
 pub mod common_moduli;
+pub mod lossy_compression;
 pub mod rotation_matrix;
diff --git a/src/utils/lossy_compression.rs b/src/utils/lossy_compression.rs
new file mode 100644
index 0000000..f393ecc
--- /dev/null
+++ b/src/utils/lossy_compression.rs
@@ -0,0 +1,392 @@
+// Copyright © 2025 Niklas Siemer
+//
+// This file is part of qFALL-tools.
+//
+// qFALL-tools is free software: you can redistribute it and/or modify it under
+// the terms of the Mozilla Public License Version 2.0 as published by the
+// Mozilla Foundation. See <https://mozilla.org/en-US/MPL/2.0/>.
+
+//! Implements lossy (de-)compression as specified in ML-KEM.
+//!
+//! Reference:
+//! - \[1\] National Institute of Standards and Technology (2024).
+//!   Module-Lattice-Based Key-Encapsulation Mechanism Standard.
+//!   Federal Information Processing Standards Publication (FIPS 203).
+//!   <https://doi.org/10.6028/NIST.FIPS.203>
+
+use flint_sys::fmpz_poly::fmpz_poly_set_coeff_fmpz;
+use qfall_math::{
+    integer::{MatPolyOverZ, PolyOverZ, Z},
+    integer_mod_q::{MatPolynomialRingZq, ModulusPolynomialRingZq, PolynomialRingZq},
+    traits::{GetCoefficient, MatrixDimensions, MatrixGetEntry, MatrixSetEntry, Pow},
+};
+
+/// This trait is implemented by data-structures, which may use lossy compression by dropping lower order bits
+/// as specified in [\[1\]](<index.html#:~:text=[1]>).
+pub trait LossyCompressionFIPS203 {
+    /// Defines the datatype that the compressed value will have.
+    type CompressedType;
+    /// Defines the type of the modulus object.
+    type ModulusType;
+
+    /// Compresses by keeping only `d` higher-order bits.
+    /// This function modifies the value of `self` directly.
+    ///
+    /// The function is specified by `Compress_d(x) := ⌈(2^d / q) * x⌋ mod 2^d`.
+    ///
+    /// Parameters:
+    /// - `d`: specifies the number of bits that is kept to represent values
+    ///
+    /// Returns a new instance of type [`Self::CompressedType`] containing the compressed coefficients with a loss-factor
+    /// defined by `q` and `d`.
+    ///
+    /// # Panics ...
+    /// - if `d` is smaller than `1`.
+    fn lossy_compress(&self, d: impl Into<Z>) -> Self::CompressedType;
+
+    /// Decompresses a previously compressed value by mapping it to the closest recoverable value over the ring `Z_q`.
+    /// This function modifies the value of `self` directly.
+    ///
+    /// The function is specified by `Decompress_d(y) := ⌈(q / 2^d) * y⌋`.
+    ///
+    /// Parameters:
+    /// - `d`: specifies the number of bits that was kept during compression
+    ///
+    /// Returns a new instance of type [`Self`] with decompressed values according to the loss-factor
+    /// defined by `q` and `d`.
+    ///
+    /// # Panics ...
+    /// - if `d` is smaller than `1`.
+    fn lossy_decompress(
+        compressed: &Self::CompressedType,
+        d: impl Into<Z>,
+        modulus: &Self::ModulusType,
+    ) -> Self;
+}
+
+impl LossyCompressionFIPS203 for PolynomialRingZq {
+    type CompressedType = PolyOverZ;
+    type ModulusType = ModulusPolynomialRingZq;
+
+    /// Compresses by keeping only `d` higher-order bits.
+    /// This function modifies the value of `self` directly.
+    ///
+    /// The function is specified by `Compress_d(x) := ⌈(2^d / q) * x⌋ mod 2^d`.
+    ///
+    /// Parameters:
+    /// - `d`: specifies the number of bits that is kept to represent each value
+    ///
+    /// Returns a [`PolyOverZ`] containing the compressed coefficients with a loss-factor
+    /// defined by `q` and `d`.
+    ///
+    /// # Examples
+    /// ```
+    /// use qfall_math::integer_mod_q::PolynomialRingZq;
+    /// use qfall_tools::utils::{common_moduli::new_anticyclic, lossy_compression::LossyCompressionFIPS203};
+    ///
+    /// let modulus = new_anticyclic(16, 257).unwrap();
+    /// let mut poly = PolynomialRingZq::sample_uniform(&modulus);
+    ///
+    /// let compressed = poly.lossy_compress(4);
+    /// ```
+    ///
+    /// # Panics ...
+    /// - if `d` is smaller than `1`.
+    fn lossy_compress(&self, d: impl Into<Z>) -> Self::CompressedType {
+        let d = d.into();
+        assert!(d >= Z::ONE, "Performing this function with d < 1 implies reducing mod 1, leaving no information to recover. Choose a larger parameter d.");
+        let two_pow_d = Z::from(2).pow(d).unwrap();
+        let q = self.get_mod().get_q();
+        let q_div_2 = q.div_floor(2);
+
+        let mut out = PolyOverZ::default();
+
+        for coeff_i in 0..=self.get_degree() {
+            let mut coeff: Z = unsafe { self.get_coeff_unchecked(coeff_i) };
+
+            coeff *= &two_pow_d;
+            coeff += &q_div_2;
+            let mut res = coeff.div_floor(&q) % &two_pow_d;
+
+            unsafe {
+                fmpz_poly_set_coeff_fmpz(out.get_fmpz_poly_struct(), coeff_i, res.get_fmpz());
+            };
+        }
+
+        out
+    }
+
+    /// Decompresses a previously compressed value by mapping it to the closest recoverable value over the ring `Z_q`.
+    /// This function modifies the value of `self` directly.
+    ///
+    /// The function is specified by `Decompress_d(y) := ⌈(q / 2^d) * y⌋`.
+    ///
+    /// Parameters:
+    /// - `compressed`: specifies the compressed value
+    /// - `d`: specifies the number of bits that was kept during compression
+    /// - `modulus`: specifies the modulus of the returned value
+    ///
+    /// Returns a [`PolynomialRingZq`] with decompressed values according to the loss-factor
+    /// defined by `q` and `d`.
+    ///
+    /// # Examples
+    /// ```
+    /// use qfall_math::integer_mod_q::PolynomialRingZq;
+    /// use qfall_tools::utils::{common_moduli::new_anticyclic, lossy_compression::LossyCompressionFIPS203};
+    ///
+    /// let modulus = new_anticyclic(16, 257).unwrap();
+    /// let mut poly = PolynomialRingZq::sample_uniform(&modulus);
+    ///
+    /// let compressed = poly.lossy_compress(4);
+    /// let decompressed = PolynomialRingZq::lossy_decompress(&compressed, 4, &modulus);
+    /// ```
+    ///
+    /// # Panics ...
+    /// - if `d` is smaller than `1`.
+    fn lossy_decompress(
+        compressed: &Self::CompressedType,
+        d: impl Into<Z>,
+        modulus: &Self::ModulusType,
+    ) -> Self {
+        let d = d.into();
+        assert!(d >= Z::ONE, "Performing this function with d < 1 implies reducing mod 1, leaving no information to recover. Choose a larger parameter d.");
+        let two_pow_d_minus_1 = Z::from(2).pow(d - 1).unwrap();
+        let two_pow_d = &two_pow_d_minus_1 * 2;
+        let q = modulus.get_q();
+
+        let mut out = Self::from(modulus);
+
+        for coeff_i in 0..=compressed.get_degree() {
+            let mut coeff: Z = unsafe { compressed.get_coeff_unchecked(coeff_i) };
+
+            coeff *= &q;
+            coeff += &two_pow_d_minus_1;
+            let mut res = coeff.div_floor(&two_pow_d);
+
+            unsafe {
+                fmpz_poly_set_coeff_fmpz(out.get_fmpz_poly_struct(), coeff_i, res.get_fmpz());
+            };
+        }
+
+        out
+    }
+}
+
+impl LossyCompressionFIPS203 for MatPolynomialRingZq {
+    type CompressedType = MatPolyOverZ;
+    type ModulusType = ModulusPolynomialRingZq;
+
+    /// Compresses by keeping only `d` higher-order bits.
+    /// This function modifies the value of `self` directly.
+    ///
+    /// The function is specified by `Compress_d(x) := ⌈(2^d / q) * x⌋ mod 2^d`.
+    ///
+    /// Parameters:
+    /// - `d`: specifies the number of bits that is kept to represent each value
+    ///
+    /// Returns a [`MatPolyOverZ`] containing the compressed coefficients with a loss-factor
+    /// defined by `q` and `d`.
+    ///
+    /// # Examples
+    /// ```
+    /// use qfall_math::integer_mod_q::MatPolynomialRingZq;
+    /// use qfall_tools::utils::{common_moduli::new_anticyclic, lossy_compression::LossyCompressionFIPS203};
+    ///
+    /// let modulus = new_anticyclic(16, 257).unwrap();
+    /// let mut poly = MatPolynomialRingZq::sample_uniform(2, 3, &modulus);
+    ///
+    /// let compressed = poly.lossy_compress(4);
+    /// ```
+    ///
+    /// # Panics ...
+    /// - if `d` is smaller than `1`.
+    fn lossy_compress(&self, d: impl Into<Z>) -> Self::CompressedType {
+        let d = d.into();
+
+        let mut out = MatPolyOverZ::new(self.get_num_rows(), self.get_num_columns());
+
+        for row in 0..self.get_num_rows() {
+            for col in 0..self.get_num_columns() {
+                let entry: PolynomialRingZq = unsafe { self.get_entry_unchecked(row, col) };
+                let res = entry.lossy_compress(&d);
+                unsafe { out.set_entry_unchecked(row, col, res) };
+            }
+        }
+
+        out
+    }
+
+    /// Decompresses a previously compressed value by mapping it to the closest recoverable value over the ring `Z_q`.
+    /// This function modifies the value of `self` directly.
+    ///
+    /// The function is specified by `Decompress_d(y) := ⌈(q / 2^d) * y⌋`.
+    ///
+    /// Parameters:
+    /// - `compressed`: specifies the compressed matrix
+    /// - `d`: specifies the number of bits that was kept during compression
+    /// - `modulus`: specifies the modulus of the returned value
+    ///
+    /// Returns a [`MatPolynomialRingZq`] with decompressed values according to the loss-factor
+    /// defined by `q` and `d`.
+    ///
+    /// # Examples
+    /// ```
+    /// use qfall_math::integer_mod_q::MatPolynomialRingZq;
+    /// use qfall_tools::utils::{common_moduli::new_anticyclic, lossy_compression::LossyCompressionFIPS203};
+    ///
+    /// let modulus = new_anticyclic(16, 257).unwrap();
+    /// let mut poly = MatPolynomialRingZq::sample_uniform(2, 3, &modulus);
+    ///
+    /// let compressed = poly.lossy_compress(4);
+    /// let decompressed = MatPolynomialRingZq::lossy_decompress(&compressed, 4, &modulus);
+    /// ```
+    ///
+    /// # Panics ...
+    /// - if `d` is smaller than `1`.
+    fn lossy_decompress(
+        compressed: &Self::CompressedType,
+        d: impl Into<Z>,
+        modulus: &Self::ModulusType,
+    ) -> Self {
+        let d = d.into();
+
+        let mut out = MatPolynomialRingZq::new(
+            compressed.get_num_rows(),
+            compressed.get_num_columns(),
+            modulus,
+        );
+
+        for row in 0..compressed.get_num_rows() {
+            for col in 0..compressed.get_num_columns() {
+                let entry: PolyOverZ = unsafe { compressed.get_entry_unchecked(row, col) };
+                let res = PolynomialRingZq::lossy_decompress(&entry, &d, modulus);
+                unsafe { out.set_entry_unchecked(row, col, res) };
+            }
+        }
+
+        out
+    }
+}
+
+#[cfg(test)]
+mod test_compression_poly {
+    use crate::utils::{common_moduli::new_anticyclic, lossy_compression::LossyCompressionFIPS203};
+    use qfall_math::{
+        integer::Z,
+        integer_mod_q::PolynomialRingZq,
+        traits::{Distance, GetCoefficient, Pow},
+    };
+
+    /// Ensures that decompressing compressed values results in values close to the original for small d.
+    #[test]
+    fn round_trip_small_d() {
+        let d = 4;
+        let q = Z::from(257);
+        let modulus = new_anticyclic(16, &q).unwrap();
+        let poly = PolynomialRingZq::sample_uniform(&modulus);
+
+        let compressed = poly.lossy_compress(d);
+        let decompressed = PolynomialRingZq::lossy_decompress(&compressed, d, &modulus);
+
+        for i in 0..modulus.get_degree() {
+            let orig_coeff: Z = poly.get_coeff(i).unwrap();
+            let coeff: Z = decompressed.get_coeff(i).unwrap();
+
+            let mut distance = orig_coeff.distance(coeff);
+            if distance > &q / 2 {
+                distance = &q - distance;
+            }
+
+            assert!(distance <= Z::from(2).pow(q.log_ceil(2).unwrap() - d - 1).unwrap());
+        }
+    }
+
+    /// Ensures that decompressing compressed values results in values close to the original.
+    #[test]
+    fn round_trip() {
+        let d = 11;
+        let q = Z::from(3329);
+        let modulus = new_anticyclic(16, &q).unwrap();
+        let poly = PolynomialRingZq::sample_uniform(&modulus);
+
+        let compressed = poly.lossy_compress(d);
+        let decompressed = PolynomialRingZq::lossy_decompress(&compressed, d, &modulus);
+
+        for i in 0..modulus.get_degree() {
+            let orig_coeff: Z = poly.get_coeff(i).unwrap();
+            let coeff: Z = decompressed.get_coeff(i).unwrap();
+
+            let mut distance = orig_coeff.distance(coeff);
+            if distance > &q / 2 {
+                distance = &q - distance;
+            }
+
+            assert!(distance <= Z::from(2).pow(q.log_ceil(2).unwrap() - d - 1).unwrap());
+        }
+    }
+
+    /// Ensures that the function panics if `d = 0` or smaller.
+    #[test]
+    #[should_panic]
+    fn too_small_d() {
+        let d = 0;
+        let q = Z::from(3329);
+        let modulus = new_anticyclic(16, &q).unwrap();
+        let poly = PolynomialRingZq::sample_uniform(&modulus);
+
+        poly.lossy_compress(d);
+    }
+}
+
+#[cfg(test)]
+mod test_compression_matrix {
+    use crate::utils::{common_moduli::new_anticyclic, lossy_compression::LossyCompressionFIPS203};
+    use qfall_math::{
+        integer::Z,
+        integer_mod_q::{MatPolynomialRingZq, PolynomialRingZq},
+        traits::{Distance, GetCoefficient, MatrixDimensions, MatrixGetEntry, Pow},
+    };
+
+    /// Ensures that decompressing compressed values results in values close to the original.
+    #[test]
+    fn round_trip() {
+        let d = 11;
+        let q = Z::from(3329);
+        let modulus = new_anticyclic(16, &q).unwrap();
+        let matrix = MatPolynomialRingZq::sample_uniform(2, 2, &modulus);
+
+        let compressed = matrix.lossy_compress(d);
+        let decompressed = MatPolynomialRingZq::lossy_decompress(&compressed, d, &modulus);
+
+        for row in 0..matrix.get_num_rows() {
+            for col in 0..matrix.get_num_columns() {
+                let orig_entry: PolynomialRingZq = matrix.get_entry(row, col).unwrap();
+                let entry: PolynomialRingZq = decompressed.get_entry(row, col).unwrap();
+
+                for i in 0..modulus.get_degree() {
+                    let orig_coeff: Z = orig_entry.get_coeff(i).unwrap();
+                    let coeff: Z = entry.get_coeff(i).unwrap();
+
+                    let mut distance = orig_coeff.distance(coeff);
+                    if distance > &q / 2 {
+                        distance = &q - distance;
+                    }
+
+                    assert!(distance <= Z::from(2).pow(q.log_ceil(2).unwrap() - d - 1).unwrap());
+                }
+            }
+        }
+    }
+
+    /// Ensures that the function panics if `d = 0` or smaller.
+    #[test]
+    #[should_panic]
+    fn too_small_d() {
+        let d = 0;
+        let q = Z::from(3329);
+        let modulus = new_anticyclic(16, &q).unwrap();
+        let poly = MatPolynomialRingZq::sample_uniform(2, 3, &modulus);
+
+        poly.lossy_compress(d);
+    }
+}