test(core): add normality test based on Shapiro-Francia

tfhe/Cargo.toml

@@ -18,7 +18,6 @@ rust-version = "1.67"
 [dev-dependencies]
 rand = "0.8.5"
 rand_distr = "0.4.3"
-kolmogorov_smirnov = "1.1.0"
 paste = "1.0.7"
 lazy_static = { version = "1.4.0" }
 criterion = "0.4.0"
@@ -29,6 +28,8 @@ bincode = "1.3.3"
 fs2 = { version = "0.4.3" }
 itertools = "0.10.5"
 num_cpus = "1.15"
+# For erf and normality test
+libm = "0.2.6"
 
 [build-dependencies]
 cbindgen = { version = "0.24.3", optional = true }

tfhe/src/core_crypto/commons/generators/encryption.rs

@@ -562,7 +562,7 @@ mod test {
         PolynomialSize,
     };
     use crate::core_crypto::commons::test_tools::{
-        new_encryption_random_generator, new_secret_random_generator,
+        new_encryption_random_generator, new_secret_random_generator, normality_test_f64,
     };
     use crate::core_crypto::commons::traits::UnsignedTorus;
 
@@ -728,6 +728,119 @@ mod test {
         noise_gen_slice_native::<u128>();
     }
 
+    fn test_normal_random_encryption_native<Scalar: UnsignedTorus>() {
+        const RUNS: usize = 10000;
+        const SAMPLES_PER_RUN: usize = 1000;
+        let mut rng = new_encryption_random_generator();
+        let failures: f64 = (0..RUNS)
+            .map(|_| {
+                let mut samples = vec![Scalar::ZERO; SAMPLES_PER_RUN];
+
+                rng.fill_slice_with_random_noise(&mut samples, StandardDev(f64::powi(2., -20)));
+
+                let samples: Vec<f64> = samples
+                    .iter()
+                    .copied()
+                    .map(|x| {
+                        let torus = x.into_torus();
+                        // The upper half of the torus corresponds to the negative domain when
+                        // mapping unsigned integer back to float (MSB or
+                        // sign bit is set)
+                        if torus > 0.5 {
+                            torus - 1.0
+                        } else {
+                            torus
+                        }
+                    })
+                    .collect();
+
+                if normality_test_f64(&samples, 0.05).null_hypothesis_is_valid {
+                    // If we are normal return 0, it's not a failure
+                    0.0
+                } else {
+                    1.0
+                }
+            })
+            .sum::<f64>();
+        let failure_rate = failures / (RUNS as f64);
+        println!("failure_rate: {failure_rate}");
+        // The expected failure rate even on proper gaussian is 5%, so we take a small safety margin
+        assert!(failure_rate <= 0.065);
+    }
+
+    #[test]
+    fn test_normal_random_encryption_native_u32() {
+        test_normal_random_encryption_native::<u32>();
+    }
+
+    #[test]
+    fn test_normal_random_encryption_native_u64() {
+        test_normal_random_encryption_native::<u64>();
+    }
+
+    #[test]
+    fn test_normal_random_encryption_native_u128() {
+        test_normal_random_encryption_native::<u128>();
+    }
+
+    fn test_normal_random_encryption_add_assign_native<Scalar: UnsignedTorus>() {
+        const RUNS: usize = 10000;
+        const SAMPLES_PER_RUN: usize = 1000;
+        let mut rng = new_encryption_random_generator();
+        let failures: f64 = (0..RUNS)
+            .map(|_| {
+                let mut samples = vec![Scalar::ZERO; SAMPLES_PER_RUN];
+
+                rng.unsigned_torus_slice_wrapping_add_random_noise_assign(
+                    &mut samples,
+                    StandardDev(f64::powi(2., -20)),
+                );
+
+                let samples: Vec<f64> = samples
+                    .iter()
+                    .copied()
+                    .map(|x| {
+                        let torus = x.into_torus();
+                        // The upper half of the torus corresponds to the negative domain when
+                        // mapping unsigned integer back to float (MSB or
+                        // sign bit is set)
+                        if torus > 0.5 {
+                            torus - 1.0
+                        } else {
+                            torus
+                        }
+                    })
+                    .collect();
+
+                if normality_test_f64(&samples, 0.05).null_hypothesis_is_valid {
+                    // If we are normal return 0, it's not a failure
+                    0.0
+                } else {
+                    1.0
+                }
+            })
+            .sum::<f64>();
+        let failure_rate = failures / (RUNS as f64);
+        println!("failure_rate: {failure_rate}");
+        // The expected failure rate even on proper gaussian is 5%, so we take a small safety margin
+        assert!(failure_rate <= 0.065);
+    }
+
+    #[test]
+    fn test_normal_random_encryption_add_assign_native_u32() {
+        test_normal_random_encryption_add_assign_native::<u32>();
+    }
+
+    #[test]
+    fn test_normal_random_encryption_add_assign_native_u64() {
+        test_normal_random_encryption_add_assign_native::<u64>();
+    }
+
+    #[test]
+    fn test_normal_random_encryption_add_assign_native_u128() {
+        test_normal_random_encryption_add_assign_native::<u128>();
+    }
+
     fn noise_gen_slice_custom_mod<Scalar: UnsignedTorus>(
         ciphertext_modulus: CiphertextModulus<Scalar>,
     ) {
@@ -783,6 +896,170 @@ mod test {
         noise_gen_slice_custom_mod::<u128>(CiphertextModulus::new_native());
     }
 
+    fn test_normal_random_encryption_custom_mod<Scalar: UnsignedTorus>(
+        ciphertext_modulus: CiphertextModulus<Scalar>,
+    ) {
+        const RUNS: usize = 10000;
+        const SAMPLES_PER_RUN: usize = 1000;
+        let mut rng = new_encryption_random_generator();
+        let failures: f64 = (0..RUNS)
+            .map(|_| {
+                let mut samples = vec![Scalar::ZERO; SAMPLES_PER_RUN];
+
+                rng.fill_slice_with_random_noise_custom_mod(
+                    &mut samples,
+                    StandardDev(f64::powi(2., -20)),
+                    ciphertext_modulus,
+                );
+
+                let samples: Vec<f64> = samples
+                    .iter()
+                    .copied()
+                    .map(|x| {
+                        let torus = x.into_torus();
+                        // The upper half of the torus corresponds to the negative domain when
+                        // mapping unsigned integer back to float (MSB or
+                        // sign bit is set)
+                        if torus > 0.5 {
+                            torus - 1.0
+                        } else {
+                            torus
+                        }
+                    })
+                    .collect();
+
+                if normality_test_f64(&samples, 0.05).null_hypothesis_is_valid {
+                    // If we are normal return 0, it's not a failure
+                    0.0
+                } else {
+                    1.0
+                }
+            })
+            .sum::<f64>();
+        let failure_rate = failures / (RUNS as f64);
+        println!("failure_rate: {failure_rate}");
+        // The expected failure rate even on proper gaussian is 5%, so we take a small safety margin
+        assert!(failure_rate <= 0.065);
+    }
+
+    #[test]
+    fn test_normal_random_encryption_custom_mod_u32() {
+        test_normal_random_encryption_custom_mod::<u32>(
+            CiphertextModulus::try_new_power_of_2(31).unwrap(),
+        );
+    }
+
+    #[test]
+    fn test_normal_random_encryption_custom_mod_u64() {
+        test_normal_random_encryption_custom_mod::<u64>(
+            CiphertextModulus::try_new_power_of_2(63).unwrap(),
+        );
+    }
+
+    #[test]
+    fn test_normal_random_encryption_custom_mod_u128() {
+        test_normal_random_encryption_custom_mod::<u128>(
+            CiphertextModulus::try_new_power_of_2(127).unwrap(),
+        );
+    }
+
+    #[test]
+    fn test_normal_random_encryption_native_custom_mod_u32() {
+        test_normal_random_encryption_custom_mod::<u32>(CiphertextModulus::new_native());
+    }
+
+    #[test]
+    fn test_normal_random_encryption_native_custom_mod_u64() {
+        test_normal_random_encryption_custom_mod::<u64>(CiphertextModulus::new_native());
+    }
+
+    #[test]
+    fn test_normal_random_encryption_native_custom_mod_u128() {
+        test_normal_random_encryption_custom_mod::<u128>(CiphertextModulus::new_native());
+    }
+
+    fn test_normal_random_encryption_add_assign_custom_mod<Scalar: UnsignedTorus>(
+        ciphertext_modulus: CiphertextModulus<Scalar>,
+    ) {
+        const RUNS: usize = 10000;
+        const SAMPLES_PER_RUN: usize = 1000;
+        let mut rng = new_encryption_random_generator();
+        let failures: f64 = (0..RUNS)
+            .map(|_| {
+                let mut samples = vec![Scalar::ZERO; SAMPLES_PER_RUN];
+
+                rng.unsigned_torus_slice_wrapping_add_random_noise_custom_mod_assign(
+                    &mut samples,
+                    StandardDev(f64::powi(2., -20)),
+                    ciphertext_modulus,
+                );
+
+                let samples: Vec<f64> = samples
+                    .iter()
+                    .copied()
+                    .map(|x| {
+                        let torus = x.into_torus();
+                        // The upper half of the torus corresponds to the negative domain when
+                        // mapping unsigned integer back to float (MSB or
+                        // sign bit is set)
+                        if torus > 0.5 {
+                            torus - 1.0
+                        } else {
+                            torus
+                        }
+                    })
+                    .collect();
+
+                if normality_test_f64(&samples, 0.05).null_hypothesis_is_valid {
+                    // If we are normal return 0, it's not a failure
+                    0.0
+                } else {
+                    1.0
+                }
+            })
+            .sum::<f64>();
+        let failure_rate = failures / (RUNS as f64);
+        println!("failure_rate: {failure_rate}");
+        // The expected failure rate even on proper gaussian is 5%, so we take a small safety margin
+        assert!(failure_rate <= 0.065);
+    }
+
+    #[test]
+    fn test_normal_random_encryption_add_assign_custom_mod_u32() {
+        test_normal_random_encryption_add_assign_custom_mod::<u32>(
+            CiphertextModulus::try_new_power_of_2(31).unwrap(),
+        );
+    }
+
+    #[test]
+    fn test_normal_random_encryption_add_assign_custom_mod_u64() {
+        test_normal_random_encryption_add_assign_custom_mod::<u64>(
+            CiphertextModulus::try_new_power_of_2(63).unwrap(),
+        );
+    }
+
+    #[test]
+    fn test_normal_random_encryption_add_assign_custom_mod_u128() {
+        test_normal_random_encryption_add_assign_custom_mod::<u128>(
+            CiphertextModulus::try_new_power_of_2(127).unwrap(),
+        );
+    }
+
+    #[test]
+    fn test_normal_random_encryption_add_assign_native_custom_mod_u32() {
+        test_normal_random_encryption_add_assign_custom_mod::<u32>(CiphertextModulus::new_native());
+    }
+
+    #[test]
+    fn test_normal_random_encryption_add_assign_native_custom_mod_u64() {
+        test_normal_random_encryption_add_assign_custom_mod::<u64>(CiphertextModulus::new_native());
+    }
+
+    #[test]
+    fn test_normal_random_encryption_add_assign_native_custom_mod_u128() {
+        test_normal_random_encryption_add_assign_custom_mod::<u128>(CiphertextModulus::new_native());
+    }
+
     fn mask_gen_slice_native<Scalar: UnsignedTorus>() {
         let mut gen = new_encryption_random_generator();
 

tfhe/src/core_crypto/commons/math/random/tests.rs

@@ -1,8 +1,8 @@
-use crate::core_crypto::commons::dispersion::LogStandardDev;
+use crate::core_crypto::commons::ciphertext_modulus::CiphertextModulus;
 use crate::core_crypto::commons::math::torus::UnsignedTorus;
 use crate::core_crypto::commons::test_tools::*;
 
-fn test_normal_random<T: UnsignedTorus>() {
+fn test_normal_random_three_sigma<T: UnsignedTorus>() {
     //! test if the normal random generation with std_dev is below 3*std_dev (99.7%)
 
     // settings
@@ -22,8 +22,10 @@ fn test_normal_random<T: UnsignedTorus>() {
         .zip(samples_int.iter())
         .for_each(|(out, &elt)| *out = elt.into_torus());
     for x in samples_float.iter_mut() {
+        // The upper half of the torus corresponds to the negative domain when mapping unsigned
+        // integer back to float (MSB or sign bit is set)
         if *x > 0.5 {
-            *x = 1. - *x;
+            *x -= 1.;
         }
     }
 
@@ -48,39 +50,306 @@ fn test_normal_random<T: UnsignedTorus>() {
 }
 
 #[test]
-fn test_normal_random_u32() {
-    test_normal_random::<u32>();
+fn test_normal_random_three_sigma_u32() {
+    test_normal_random_three_sigma::<u32>();
 }
 
 #[test]
-fn test_normal_random_u64() {
-    test_normal_random::<u64>();
+fn test_normal_random_three_sigma_u64() {
+    test_normal_random_three_sigma::<u64>();
 }
 
-fn test_distribution<T: UnsignedTorus>() {
-    //! tests gaussianity against the rand crate generation
-    // settings
-    let std_dev: f64 = f64::powi(2., -5);
-    let mean: f64 = 0.;
-    let k = 10_000_000;
-    let mut generator = new_random_generator();
+#[test]
+fn test_normal_random_f64() {
+    const RUNS: usize = 10000;
+    const SAMPLES_PER_RUN: usize = 1000;
+    let mut rng = new_random_generator();
+    let failures: f64 = (0..RUNS)
+        .map(|_| {
+            let mut samples = vec![0.0f64; SAMPLES_PER_RUN];
 
-    // generate normal random
-    let first = vec![T::ZERO; k];
-    let mut second = vec![T::ZERO; k];
-    generator.fill_slice_with_random_gaussian(&mut second, mean, std_dev);
+            rng.fill_slice_with_random_gaussian(&mut samples, 0.0, 1.0);
 
-    assert_noise_distribution(&first, &second, LogStandardDev(-5.));
+            if normality_test_f64(&samples, 0.05).null_hypothesis_is_valid {
+                // If we are normal return 0, it's not a failure
+                0.0
+            } else {
+                1.0
+            }
+        })
+        .sum::<f64>();
+    let failure_rate = failures / (RUNS as f64);
+    println!("failure_rate: {failure_rate}");
+    // The expected failure rate even on proper gaussian is 5%, so we take a small safety margin
+    assert!(failure_rate <= 0.065);
 }
 
-// // These tests are notoriously flaky
+fn test_normal_random_native<Scalar: UnsignedTorus>() {
+    const RUNS: usize = 10000;
+    const SAMPLES_PER_RUN: usize = 1000;
+    let mut rng = new_random_generator();
+    let failures: f64 = (0..RUNS)
+        .map(|_| {
+            let mut samples = vec![Scalar::ZERO; SAMPLES_PER_RUN];
 
-// #[test]
-// fn test_distribution_u32() {
-//     test_distribution::<u32>();
-// }
+            rng.fill_slice_with_random_gaussian(&mut samples, 0.0, f64::powi(2., -20));
 
-// #[test]
-// fn test_distribution_u64() {
-//     test_distribution::<u64>();
-// }
+            let samples: Vec<f64> = samples
+                .iter()
+                .copied()
+                .map(|x| {
+                    let torus = x.into_torus();
+                    // The upper half of the torus corresponds to the negative domain when mapping
+                    // unsigned integer back to float (MSB or sign bit is set)
+                    if torus > 0.5 {
+                        torus - 1.0
+                    } else {
+                        torus
+                    }
+                })
+                .collect();
+
+            if normality_test_f64(&samples, 0.05).null_hypothesis_is_valid {
+                // If we are normal return 0, it's not a failure
+                0.0
+            } else {
+                1.0
+            }
+        })
+        .sum::<f64>();
+    let failure_rate = failures / (RUNS as f64);
+    println!("failure_rate: {failure_rate}");
+    // The expected failure rate even on proper gaussian is 5%, so we take a small safety margin
+    assert!(failure_rate <= 0.065);
+}
+
+#[test]
+fn test_normal_random_native_u32() {
+    test_normal_random_native::<u32>();
+}
+
+#[test]
+fn test_normal_random_native_u64() {
+    test_normal_random_native::<u64>();
+}
+
+#[test]
+fn test_normal_random_native_u128() {
+    test_normal_random_native::<u128>();
+}
+
+fn test_normal_random_custom_mod<Scalar: UnsignedTorus>(
+    ciphertext_modulus: CiphertextModulus<Scalar>,
+) {
+    const RUNS: usize = 10000;
+    const SAMPLES_PER_RUN: usize = 1000;
+    let mut rng = new_random_generator();
+    let failures: f64 = (0..RUNS)
+        .map(|_| {
+            let mut samples = vec![Scalar::ZERO; SAMPLES_PER_RUN];
+
+            rng.fill_slice_with_random_gaussian_custom_mod(
+                &mut samples,
+                0.0,
+                f64::powi(2., -20),
+                ciphertext_modulus,
+            );
+
+            let samples: Vec<f64> = samples
+                .iter()
+                .copied()
+                .map(|x| {
+                    let torus = x.into_torus();
+                    // The upper half of the torus corresponds to the negative domain when mapping
+                    // unsigned integer back to float (MSB or sign bit is set)
+                    if torus > 0.5 {
+                        torus - 1.0
+                    } else {
+                        torus
+                    }
+                })
+                .collect();
+
+            if normality_test_f64(&samples, 0.05).null_hypothesis_is_valid {
+                // If we are normal return 0, it's not a failure
+                0.0
+            } else {
+                1.0
+            }
+        })
+        .sum::<f64>();
+    let failure_rate = failures / (RUNS as f64);
+    println!("failure_rate: {failure_rate}");
+    // The expected failure rate even on proper gaussian is 5%, so we take a small safety margin
+    assert!(failure_rate <= 0.065);
+}
+
+#[test]
+fn test_normal_random_custom_mod_u32() {
+    test_normal_random_custom_mod::<u32>(CiphertextModulus::try_new_power_of_2(31).unwrap());
+}
+
+#[test]
+fn test_normal_random_custom_mod_u64() {
+    test_normal_random_custom_mod::<u64>(CiphertextModulus::try_new_power_of_2(63).unwrap());
+}
+
+#[test]
+fn test_normal_random_custom_mod_u128() {
+    test_normal_random_custom_mod::<u128>(CiphertextModulus::try_new_power_of_2(127).unwrap());
+}
+
+#[test]
+fn test_normal_random_native_custom_mod_u32() {
+    test_normal_random_custom_mod::<u32>(CiphertextModulus::new_native());
+}
+
+#[test]
+fn test_normal_random_native_custom_mod_u64() {
+    test_normal_random_custom_mod::<u64>(CiphertextModulus::new_native());
+}
+
+#[test]
+fn test_normal_random_native_custom_mod_u128() {
+    test_normal_random_custom_mod::<u128>(CiphertextModulus::new_native());
+}
+
+fn test_normal_random_add_assign_native<Scalar: UnsignedTorus>() {
+    const RUNS: usize = 10000;
+    const SAMPLES_PER_RUN: usize = 1000;
+    let mut rng = new_random_generator();
+    let failures: f64 = (0..RUNS)
+        .map(|_| {
+            let mut samples = vec![Scalar::ZERO; SAMPLES_PER_RUN];
+
+            rng.unsigned_torus_slice_wrapping_add_random_gaussian_assign(
+                &mut samples,
+                0.0,
+                f64::powi(2., -20),
+            );
+
+            let samples: Vec<f64> = samples
+                .iter()
+                .copied()
+                .map(|x| {
+                    let torus = x.into_torus();
+                    // The upper half of the torus corresponds to the negative domain when mapping
+                    // unsigned integer back to float (MSB or sign bit is set)
+                    if torus > 0.5 {
+                        torus - 1.0
+                    } else {
+                        torus
+                    }
+                })
+                .collect();
+
+            if normality_test_f64(&samples, 0.05).null_hypothesis_is_valid {
+                // If we are normal return 0, it's not a failure
+                0.0
+            } else {
+                1.0
+            }
+        })
+        .sum::<f64>();
+    let failure_rate = failures / (RUNS as f64);
+    println!("failure_rate: {failure_rate}");
+    // The expected failure rate even on proper gaussian is 5%, so we take a small safety margin
+    assert!(failure_rate <= 0.065);
+}
+
+#[test]
+fn test_normal_random_add_assign_native_u32() {
+    test_normal_random_add_assign_native::<u32>();
+}
+
+#[test]
+fn test_normal_random_add_assign_native_u64() {
+    test_normal_random_add_assign_native::<u64>();
+}
+
+#[test]
+fn test_normal_random_add_assign_native_u128() {
+    test_normal_random_add_assign_native::<u128>();
+}
+
+fn test_normal_random_add_assign_custom_mod<Scalar: UnsignedTorus>(
+    ciphertext_modulus: CiphertextModulus<Scalar>,
+) {
+    const RUNS: usize = 10000;
+    const SAMPLES_PER_RUN: usize = 1000;
+    let mut rng = new_random_generator();
+    let failures: f64 = (0..RUNS)
+        .map(|_| {
+            let mut samples = vec![Scalar::ZERO; SAMPLES_PER_RUN];
+
+            rng.unsigned_torus_slice_wrapping_add_random_gaussian_custom_mod_assign(
+                &mut samples,
+                0.0,
+                f64::powi(2., -20),
+                ciphertext_modulus,
+            );
+
+            let samples: Vec<f64> = samples
+                .iter()
+                .copied()
+                .map(|x| {
+                    let torus = x.into_torus();
+                    // The upper half of the torus corresponds to the negative domain when mapping
+                    // unsigned integer back to float (MSB or sign bit is set)
+                    if torus > 0.5 {
+                        torus - 1.0
+                    } else {
+                        torus
+                    }
+                })
+                .collect();
+
+            if normality_test_f64(&samples, 0.05).null_hypothesis_is_valid {
+                // If we are normal return 0, it's not a failure
+                0.0
+            } else {
+                1.0
+            }
+        })
+        .sum::<f64>();
+    let failure_rate = failures / (RUNS as f64);
+    println!("failure_rate: {failure_rate}");
+    // The expected failure rate even on proper gaussian is 5%, so we take a small safety margin
+    assert!(failure_rate <= 0.065);
+}
+
+#[test]
+fn test_normal_random_add_assign_custom_mod_u32() {
+    test_normal_random_add_assign_custom_mod::<u32>(
+        CiphertextModulus::try_new_power_of_2(31).unwrap(),
+    );
+}
+
+#[test]
+fn test_normal_random_add_assign_custom_mod_u64() {
+    test_normal_random_add_assign_custom_mod::<u64>(
+        CiphertextModulus::try_new_power_of_2(63).unwrap(),
+    );
+}
+
+#[test]
+fn test_normal_random_add_assign_custom_mod_u128() {
+    test_normal_random_add_assign_custom_mod::<u128>(
+        CiphertextModulus::try_new_power_of_2(127).unwrap(),
+    );
+}
+
+#[test]
+fn test_normal_random_add_assign_native_custom_mod_u32() {
+    test_normal_random_add_assign_custom_mod::<u32>(CiphertextModulus::new_native());
+}
+
+#[test]
+fn test_normal_random_add_assign_native_custom_mod_u64() {
+    test_normal_random_add_assign_custom_mod::<u64>(CiphertextModulus::new_native());
+}
+
+#[test]
+fn test_normal_random_add_assign_native_custom_mod_u128() {
+    test_normal_random_add_assign_custom_mod::<u128>(CiphertextModulus::new_native());
+}

tfhe/src/core_crypto/commons/mod.rs

@@ -67,13 +67,14 @@ pub mod test_tools {
     use crate::core_crypto::commons::generators::{
         EncryptionRandomGenerator, SecretRandomGenerator,
     };
-    use crate::core_crypto::commons::math::random::{RandomGenerable, RandomGenerator, Uniform};
+    use crate::core_crypto::commons::math::random::{
+        ActivatedRandomGenerator, RandomGenerable, RandomGenerator, Uniform,
+    };
     use crate::core_crypto::commons::parameters::{
         CiphertextCount, DecompositionBaseLog, DecompositionLevelCount, GlweDimension,
         LweDimension, PlaintextCount, PolynomialSize,
     };
     use crate::core_crypto::commons::traits::*;
-    use concrete_csprng::generators::SoftwareRandomGenerator;
     use concrete_csprng::seeders::{Seed, Seeder};
 
     fn modular_distance<T: UnsignedInteger>(first: T, other: T) -> T {
@@ -94,15 +95,16 @@ pub mod test_tools {
         }
     }
 
-    pub fn new_random_generator() -> RandomGenerator<SoftwareRandomGenerator> {
+    pub fn new_random_generator() -> RandomGenerator<ActivatedRandomGenerator> {
         RandomGenerator::new(random_seed())
     }
 
-    pub fn new_secret_random_generator() -> SecretRandomGenerator<SoftwareRandomGenerator> {
+    pub fn new_secret_random_generator() -> SecretRandomGenerator<ActivatedRandomGenerator> {
         SecretRandomGenerator::new(random_seed())
     }
 
-    pub fn new_encryption_random_generator() -> EncryptionRandomGenerator<SoftwareRandomGenerator> {
+    pub fn new_encryption_random_generator() -> EncryptionRandomGenerator<ActivatedRandomGenerator>
+    {
         EncryptionRandomGenerator::new(random_seed(), &mut UnsafeRandSeeder)
     }
 
@@ -143,49 +145,85 @@ pub mod test_tools {
         }
     }
 
-    pub fn assert_noise_distribution<First, Second, Element>(
-        first: &First,
-        second: &Second,
-        dist: impl DispersionParameter,
-    ) where
-        First: Container<Element = Element>,
-        Second: Container<Element = Element>,
-        Element: UnsignedTorus,
-    {
-        use rand_distr::Distribution;
+    pub struct NormalityTestResult {
+        pub w_prime: f64,
+        pub p_value: f64,
+        pub null_hypothesis_is_valid: bool,
+    }
 
-        let std_dev = dist.get_standard_dev();
-        let confidence = 0.95;
-        let n_slots = first.container_len();
+    /// Based on Shapiro-Francia normality test
+    pub fn normality_test_f64(samples: &[f64], alpha: f64) -> NormalityTestResult {
+        assert!(
+            samples.len() <= 5000,
+            "normality_test_f64 produces a relevant pvalue for less than 5000 samples"
+        );
 
-        // allocate 2 slices: one for the error samples obtained, the second for fresh samples
-        // according to the std_dev computed
-        let mut sdk_samples = vec![0.0_f64; n_slots];
+        // From "A handy approximation for the error function and its inverse" by Sergei Winitzki
+        fn erf_inv(x: f64) -> f64 {
+            let sign = if x < 0.0 { -1.0 } else { 1.0 };
+            // 1 - x**2
+            let one_minus_x_2 = (1.0 - x) * (1.0 + x);
+            // ln(1 - x**2)
+            let log_term = f64::ln(one_minus_x_2);
+            let a = 0.147;
+            let term_1 = 2.0 / (std::f64::consts::PI * a) + 0.5 * log_term;
+            let term_2 = 1.0 / a * log_term;
 
-        // recover the errors from each ciphertexts
-        sdk_samples
-            .iter_mut()
-            .zip(first.as_ref().iter().zip(second.as_ref().iter()))
-            .for_each(|(out, (&lhs, &rhs))| *out = torus_modular_distance(lhs, rhs));
-
-        // fill the theoretical sample vector according to std_dev using the rand crate
-        let mut theoretical_samples: Vec<f64> = Vec::with_capacity(n_slots);
-        let normal = rand_distr::Normal::new(0.0, std_dev).unwrap();
-        for _i in 0..n_slots {
-            theoretical_samples.push(normal.sample(&mut rand::thread_rng()));
+            sign * f64::sqrt(-term_1 + f64::sqrt(term_1 * term_1 - term_2))
         }
 
-        // compute the kolmogorov smirnov test
-        let result = kolmogorov_smirnov::test_f64(
-            sdk_samples.as_slice(),
-            theoretical_samples.as_slice(),
-            confidence,
-        );
-        assert!(
-            !result.is_rejected,
-            "Not the same distribution with a probability of {}",
-            result.reject_probability
-        );
+        // Normal law CDF
+        fn phi(x: f64) -> f64 {
+            0.5 * (1.0 + libm::erf(x / f64::sqrt(2.0)))
+        }
+
+        fn phi_inv(x: f64) -> f64 {
+            f64::sqrt(2.0) * erf_inv(2.0 * x - 1.0)
+        }
+
+        let n = samples.len();
+        let n_f64 = n as f64;
+        // Sort the input
+        let mut samples: Vec<_> = samples.to_vec();
+        samples.sort_by(|x, y| x.partial_cmp(y).unwrap());
+        let samples = samples;
+        // Compute the mean
+        let mean = samples.iter().copied().sum::<f64>() / n_f64;
+        let frac_three_eight = 3. / 8.;
+        let frac_one_four = 1. / 4.;
+        // Compute Blom scores
+        let m_tilde: Vec<_> = (1..=n)
+            .map(|i| phi_inv((i as f64 - frac_three_eight) / (n_f64 + frac_one_four)))
+            .collect();
+        // Blom scores norm2
+        let m_norm = f64::sqrt(m_tilde.iter().fold(0.0, |acc, x| acc + x * x));
+        // Coefficients
+        let mut coeffs = m_tilde;
+        coeffs.iter_mut().for_each(|x| *x /= m_norm);
+        // Test statistic
+        let denominator = samples.iter().fold(0.0, |acc, x| acc + (x - mean).powi(2));
+        let numerator = samples
+            .iter()
+            .zip(coeffs.iter())
+            .fold(0.0, |acc, (&sample, &coeff)| acc + sample * coeff)
+            .powi(2);
+        let w_prime = numerator / denominator;
+
+        let g_w_prime = f64::ln(1.0 - w_prime);
+        let log_n = n_f64.ln();
+        let log_log_n = log_n.ln();
+        let u = log_log_n - log_n;
+        let mu = 1.0521 * u - 1.2725;
+        let v = log_log_n + 2.0 / log_n;
+        let sigma = -0.26758 * v + 1.0308;
+        let z = (g_w_prime - mu) / sigma;
+        let p_value = 1.0 - phi(z);
+
+        NormalityTestResult {
+            w_prime,
+            p_value,
+            null_hypothesis_is_valid: p_value > alpha,
+        }
     }
 
     /// Return a random plaintext count in [1;max].
@@ -260,4 +298,52 @@ pub mod test_tools {
         let mut generator = new_random_generator();
         generator.random_uniform()
     }
+
+    #[test]
+    pub fn test_normality_tool() {
+        use rand_distr::{Distribution, Normal};
+        const RUNS: usize = 10000;
+        const SAMPLES_PER_RUN: usize = 1000;
+        let mut rng = rand::thread_rng();
+        let normal = Normal::new(0.0, 1.0).unwrap();
+        let failures: f64 = (0..RUNS)
+            .map(|_| {
+                let mut samples = vec![0.0f64; SAMPLES_PER_RUN];
+                samples
+                    .iter_mut()
+                    .for_each(|x| *x = normal.sample(&mut rng));
+                if normality_test_f64(&samples, 0.05).null_hypothesis_is_valid {
+                    // If we are normal return 0, it's not a failure
+                    0.0
+                } else {
+                    1.0
+                }
+            })
+            .sum::<f64>();
+        let failure_rate = failures / (RUNS as f64);
+        println!("failure_rate: {failure_rate}");
+        // The expected failure rate even on proper gaussian is 5%, so we take a small safety margin
+        assert!(failure_rate <= 0.065);
+    }
+
+    #[test]
+    pub fn test_normality_tool_fail_uniform() {
+        const RUNS: usize = 10000;
+        const SAMPLES_PER_RUN: usize = 1000;
+        let mut rng = rand::thread_rng();
+        let failures: f64 = (0..RUNS)
+            .map(|_| {
+                let mut samples = vec![0.0f64; SAMPLES_PER_RUN];
+                samples.iter_mut().for_each(|x| *x = rng.gen());
+                if normality_test_f64(&samples, 0.05).null_hypothesis_is_valid {
+                    // If we are normal return 0, it's not a failure
+                    0.0
+                } else {
+                    1.0
+                }
+            })
+            .sum::<f64>();
+        let failure_rate = failures / (RUNS as f64);
+        assert!(failure_rate == 1.0);
+    }
 }