diff --git a/Cargo.lock b/Cargo.lock index 98e140e..49827b6 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -60,12 +60,6 @@ dependencies = [ "miniz_oxide", ] -[[package]] -name = "fuchsia-cprng" -version = "0.1.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a06f77d526c1a601b7c4cdd98f54b5eaabffc14d5f2f0296febdc7f357c6d3ba" - [[package]] name = "getrandom" version = "0.2.15" @@ -89,15 +83,6 @@ version = "2.0.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b248f5224d1d606005e02c97f5aa4e88eeb230488bcc03bc9ca4d7991399f2b5" -[[package]] -name = "kolmogorov_smirnov" -version = "1.1.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "96ceadb630ffd3ceb5c4d739232d8a9d51d23c383b14f3b64ffc67368cbc59b5" -dependencies = [ - "rand 0.3.23", -] - [[package]] name = "libc" version = "0.2.169" @@ -149,7 +134,7 @@ dependencies = [ "num-complex", "num-rational", "num-traits", - "rand 0.8.5", + "rand", "rand_distr", "simba", "typenum", @@ -215,7 +200,6 @@ name = "onod3000" version = "0.1.0" dependencies = [ "flate2", - "kolmogorov_smirnov", "pyo3", "statrs", ] @@ -322,29 +306,6 @@ dependencies = [ "proc-macro2", ] -[[package]] -name = "rand" -version = "0.3.23" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "64ac302d8f83c0c1974bf758f6b041c6c8ada916fbb44a609158ca8b064cc76c" -dependencies = [ - "libc", - "rand 0.4.6", -] - -[[package]] -name = "rand" -version = "0.4.6" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "552840b97013b1a26992c11eac34bdd778e464601a4c2054b5f0bff7c6761293" -dependencies = [ - "fuchsia-cprng", - "libc", - "rand_core 0.3.1", - "rdrand", - "winapi", -] - [[package]] name = "rand" version = "0.8.5" @@ -353,7 +314,7 @@ checksum = "34af8d1a0e25924bc5b7c43c079c942339d8f0a8b57c39049bef581b46327404" dependencies = [ "libc", "rand_chacha", - "rand_core 0.6.4", + "rand_core", ] [[package]] @@ -363,24 +324,9 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88" dependencies = [ "ppv-lite86", - "rand_core 0.6.4", + "rand_core", ] -[[package]] -name = "rand_core" -version = "0.3.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7a6fdeb83b075e8266dcc8762c22776f6877a63111121f5f8c7411e5be7eed4b" -dependencies = [ - "rand_core 0.4.2", -] - -[[package]] -name = "rand_core" -version = "0.4.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9c33a3c44ca05fa6f1807d8e6743f3824e8509beca625669633be0acbdf509dc" - [[package]] name = "rand_core" version = "0.6.4" @@ -397,7 +343,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "32cb0b9bc82b0a0876c2dd994a7e7a2683d3e7390ca40e6886785ef0c7e3ee31" dependencies = [ "num-traits", - "rand 0.8.5", + "rand", ] [[package]] @@ -406,15 +352,6 @@ version = "0.2.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" -[[package]] -name = "rdrand" -version = "0.4.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "678054eb77286b51581ba43620cc911abf02758c91f93f479767aed0f90458b2" -dependencies = [ - "rand_core 0.3.1", -] - [[package]] name = "safe_arch" version = "0.7.4" @@ -446,7 +383,7 @@ dependencies = [ "approx", "nalgebra", "num-traits", - "rand 0.8.5", + "rand", ] [[package]] @@ -500,28 +437,6 @@ dependencies = [ "safe_arch", ] -[[package]] -name = "winapi" -version = "0.3.9" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5c839a674fcd7a98952e593242ea400abe93992746761e38641405d28b00f419" -dependencies = [ - "winapi-i686-pc-windows-gnu", - "winapi-x86_64-pc-windows-gnu", -] - -[[package]] -name = "winapi-i686-pc-windows-gnu" -version = "0.4.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ac3b87c63620426dd9b991e5ce0329eff545bccbbb34f3be09ff6fb6ab51b7b6" - -[[package]] -name = "winapi-x86_64-pc-windows-gnu" -version = "0.4.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" - [[package]] name = "zerocopy" version = "0.7.35" diff --git a/Cargo.toml b/Cargo.toml index b625ad1..55410fb 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -9,7 +9,6 @@ crate-type = ["cdylib", "staticlib", "lib"] [dependencies] flate2 = "1.0.35" -kolmogorov_smirnov = "1.1.0" pyo3 = {version = "0.23.4", optional = true} statrs = "0.18.0" diff --git a/readme.md b/readme.md index d90181f..ba83174 100644 --- a/readme.md +++ b/readme.md @@ -103,7 +103,7 @@ ChiBit | 0.9299 | 0.9299 | | ChiByte | 0.9819 | 0.9819 | | MeanByte | 0.7964 | 0.7964 | | Compression | N/A | 0.0146 | | -KS | 0.5313 | 0.2476 | Mismatch | +KS | 0.5313 | 0.9998 | Different method | Pi | 0.3784 | 0.3784 | | Shells | 0.7559 | 0.7559 | | Gaps | 0.0000 | 0.0000 | | @@ -125,7 +125,7 @@ ChiBit | 0.5590 | 0.5590 | | ChiByte | 0.3748 | 0.3748 | | MeanByte | 0.0183 | 0.0183 | | Compression | N/A | 0.9308 | | -KS | 0.2900 | 0.0178 | Mismatch | +KS | 0.2900 | 0.1158 | Different method | Pi | 0.9765 | 0.9765 | | Shells | 0.7235 | 0.7235 | | Gaps | 0.6425 | 0.6425 | | @@ -145,7 +145,7 @@ ChiBit | 0.8428 | 0.8428 | | ChiByte | 0.5785 | 0.5785 | | MeanByte | 0.9601 | 0.9601 | | Compression | N/A | 0.9308 | | -KS | 0.7395 | 0.0441 | Mismatch | +KS | 0.7395 | 0.7777 | Different method | Pi | 0.2806 | 0.2806 | | Shells | 0.7711 | 0.7711 | | Gaps | 0.1937 | 0.1937 | | @@ -157,14 +157,9 @@ UnCorrelation | 0.4674 | 0.4674 | | ### Conclusions -Although the implementation of the randomness tests closely follows the original Java logic, differences in the results arise in the Kolmogorov-Smirnov test. The external libraries used for running the test probably differ in their implementation. Also the uniform distribution to test against is provided by the Apache Commons Math library. Since it would be way beyond the scope of this porting project to try to fully match the functionality of the dependencies used we just accept the minor difference. Another issue might be that Java emphasizes predictability and portability, enforcing strict IEEE 754 behavior across platforms. Rust on the other hand prioritizes performance and flexibility, allowing platform-specific optimizations that may deviate slightly from strict IEEE semantics. +Although the implementation of the randomness tests closely follows the original Java logic, the Kolmogorov-Smirnov test deliberately differs. The Java version compares the data against a pseudo-random uniform sample, which makes the p-value change from run to run for the same input. This port instead runs a one-sample KS test of the byte values against the exact discrete uniform distribution over 0..255, with p-values from the asymptotic Kolmogorov distribution (Stephens' small-sample correction). The result is deterministic and calibrated: on `/dev/urandom` it fails about 1% of the time at alpha = 0.01, slightly less on small inputs because the discrete null makes the test a little conservative. -These differences are generally negligible for practical purposes and do not affect the overall functionality or -statistical significance of the test. - -### Bonus - -We implemented a Well Equidistributed Long-period Linear pseudo-random number generator that is used with a random seed derived from an epoch timestamp that is being used as uniform distribution for the KS test. +All other differences are generally negligible for practical purposes and do not affect the overall functionality or statistical significance of the tests. ## License diff --git a/src/uniformity/ks.rs b/src/uniformity/ks.rs index 9c76e9a..7b3a74f 100644 --- a/src/uniformity/ks.rs +++ b/src/uniformity/ks.rs @@ -23,79 +23,104 @@ // SOFTWARE. -use kolmogorov_smirnov::test_f64; -use std::time::{SystemTime, UNIX_EPOCH}; -use statrs::distribution::ContinuousCDF; - use crate::Onod; -use crate::well19937c::Well19937c; impl Onod { /// KS randomness test - /// Performs the Kolmogorov-Smirnov test to evaluate the uniformity of data distribution - /// and returns the test statistic (D-statistic), z-score, and p-value. + /// Performs a one-sample Kolmogorov-Smirnov test of the byte values against the + /// exact discrete uniform distribution over 0..=255 and returns the test statistic + /// (D-statistic), the scaled statistic sqrt(n)*D, and the p-value. pub fn ks(samples: &[u8]) -> (f64, f64, f64) { if samples.is_empty() { return (-1.0, 0.0, 1.0); // Invalid input } - - // Normalize the input samples to [0, 1) range - let mut normalized_samples: Vec = samples.iter().map(|&x| x as f64 / 255.0).collect(); - - // Generate a uniform distribution for comparison - let mut uniform_distribution: Vec = Self::generate_uniform_distribution(normalized_samples.len(), Self::get_timestamp_seed()); - - // Sort both distributions - normalized_samples.sort_by(|a, b| a.partial_cmp(b).unwrap()); - uniform_distribution.sort_by(|a, b| a.partial_cmp(b).unwrap()); - - let debug = false; - // Debugging: Optional print sorted values - if debug { - println!("Sorted Normalized Samples: {:?}", normalized_samples); - println!("Sorted Uniform Distribution: {:?}", uniform_distribution); + + let mut counts = [0u64; 256]; + for &b in samples { + counts[b as usize] += 1; } - - // Perform the Kolmogorov-Smirnov test - let confidence = 0.05; // Significance level - let result = test_f64(&normalized_samples, &uniform_distribution, confidence); - - // Extract the KS statistic (D-statistic) - let ks_statistic = result.statistic; - - // Debugging: Print ECDF differences and max difference (D-statistic) - if debug { - for (i, (&sample, &uniform)) in normalized_samples.iter().zip(&uniform_distribution).enumerate() { - let diff = (sample - uniform).abs(); - println!( - "Index: {}, Sample: {:.6}, Uniform: {:.6}, Difference: {:.6}", - i, sample, uniform, diff - ); - } - println!("D-Statistic: {:.6}", ks_statistic); + + // Both CDFs are step functions that jump only at the byte values, so the + // supremum of their difference is reached at one of those 256 points. + let n = samples.len() as f64; + let mut cumulative = 0u64; + let mut ks_statistic: f64 = 0.0; + for (k, &count) in counts.iter().enumerate() { + cumulative += count; + let empirical = cumulative as f64 / n; + let expected = (k + 1) as f64 / 256.0; + ks_statistic = ks_statistic.max((empirical - expected).abs()); } - - // Calculate the z-score - let sample_size = normalized_samples.len() as f64; - let z_score = ks_statistic * sample_size.sqrt(); - - // Calculate the p-value - let p_value = 2.0 * (1.0 - statrs::distribution::Normal::new(0.0, 1.0).unwrap().cdf(z_score.abs())); - + + let sqrt_n = n.sqrt(); + let z_score = ks_statistic * sqrt_n; + + // Asymptotic Kolmogorov distribution with Stephens' small-sample correction. + // For a discrete null this is slightly conservative. + let p_value = Self::kolmogorov_survival((sqrt_n + 0.12 + 0.11 / sqrt_n) * ks_statistic); + (ks_statistic, z_score, p_value) } - /// Generates a uniform distribution of the same length as the input data. - fn generate_uniform_distribution(len: usize, seed: u32) -> Vec { - let mut rng = Well19937c::new(seed); - (0..len).map(|_| rng.next_f64()).collect() + /// P(K > lambda) for the Kolmogorov distribution. + fn kolmogorov_survival(lambda: f64) -> f64 { + if lambda <= 0.0 { + return 1.0; + } + if lambda < 1.18 { + // P(K <= lambda) = sqrt(2*pi)/lambda * sum_k exp(-(2k-1)^2 * pi^2 / (8*lambda^2)) + let w = -std::f64::consts::PI.powi(2) / (8.0 * lambda * lambda); + let cdf: f64 = (1..=8) + .map(|k| (((2 * k - 1) as f64).powi(2) * w).exp()) + .sum::() + * (2.0 * std::f64::consts::PI).sqrt() + / lambda; + (1.0 - cdf).clamp(0.0, 1.0) + } else { + // P(K > lambda) = 2 * sum_k (-1)^(k-1) * exp(-2 * k^2 * lambda^2) + let survival: f64 = (1..=100) + .map(|k| { + let sign = if k % 2 == 1 { 1.0 } else { -1.0 }; + sign * (-2.0 * (k * k) as f64 * lambda * lambda).exp() + }) + .sum::() + * 2.0; + survival.clamp(0.0, 1.0) + } + } +} + +#[cfg(test)] +mod tests { + use crate::Onod; + + #[test] + fn perfectly_uniform_bytes_pass() { + let samples: Vec = (0..1024).map(|i| i as u8).collect(); + let (d, _, p) = Onod::ks(&samples); + assert_eq!(d, 0.0); + assert_eq!(p, 1.0); } - fn get_timestamp_seed() -> u32 { - let duration = SystemTime::now() - .duration_since(UNIX_EPOCH) - .expect("Time went backwards"); - // Use seconds or nanoseconds as the seed - (duration.as_secs() as u32) ^ (duration.subsec_nanos()) + #[test] + fn skewed_bytes_fail() { + let samples: Vec = (0..4096).map(|i| (i % 128) as u8).collect(); + let (_, _, p) = Onod::ks(&samples); + assert!(p < 1e-6); } -} \ No newline at end of file + + #[test] + fn deterministic() { + let samples: Vec = (0..10_000u32).map(|i| (i.wrapping_mul(2654435761) >> 24) as u8).collect(); + assert_eq!(Onod::ks(&samples), Onod::ks(&samples)); + } + + #[test] + fn kolmogorov_reference_values() { + // Critical values of the Kolmogorov distribution. + assert!((Onod::kolmogorov_survival(1.3581) - 0.05).abs() < 1e-4); + assert!((Onod::kolmogorov_survival(1.6276) - 0.01).abs() < 1e-4); + // Both series agree at the switch point. + assert!((Onod::kolmogorov_survival(1.18 - 1e-12) - Onod::kolmogorov_survival(1.18)).abs() < 1e-9); + } +}