randogoth/onod3000
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

fixes and a PRNG

Browse source
This commit is contained in:
randogoth 2025-01-18 14:22:17 +02:00
parent b4c2e5ea9b
commit 2d6b61b6cf
9 changed files with 336 additions and 203 deletions
Download patch file Download diff file Expand all files Collapse all files

74
readme.md
View file

@ -96,42 +96,44 @@ To test the implementation three files of binary random data have been analyzed
### Testing 1024 bytes from test1.bin (/dev/random) ### Testing 1024 bytes from test1.bin (/dev/random)
Test | Java ent3000 | Rust onod3000 | Note | Test | Java ent3000 | Rust onod3000 | Note |
------------- | ------------ | ------------- | ---- | ------------- | ------------ | ------------- | ------------- |
Shannon | N/A | 0.0000 | | Shannon | N/A | 0.0000 | |
Monobit | 0.8771 | 0.8771 | | Monobit | 0.6746 | 0.6746 | |
ChiBit | 0.0949 | 0.0949 | | ChiBit | 0.9299 | 0.9299 | |
ChiByte | 0.1307 | 0.1307 | | ChiByte | 0.9819 | 0.9819 | |
MeanByte | 0.5378 | 0.5378 | | MeanByte | 0.7964 | 0.7964 | |
Compression | N/A | 0.0146 | | Compression | N/A | 0.0146 | |
KS | 0.5142 | 0.9814 | MISMATCH | KS | 0.5313 | 0.2476 | Mismatch |
Pi | 0.4015 | 0.7293 | MISMATCH | Pi | 0.3784 | 0.3784 | |
Shells | 0.5479 | 0.1470 | MISMATCH | Shells | 0.7559 | 0.7559 | |
Gaps | 0.0000 | 0.0000 | | Gaps | 0.0000 | 0.0000 | |
Avalanche | 0.9798 | 0.9798 | | Avalanche | 0.8794 | 0.8794 | |
Runs | 0.4916 | 0.4916 | | Runs | 0.0850 | 0.0850 | |
RunUps | 0.2573 | 0.2573 | | RunUps | 0.6537 | 0.6537 | |
Prediction | 0.3164 | 0.1255 | MISMATCH | Prediction | 0.6164 | 0.6164 | |
UnCorrelation | 0.4556 | 0.4554 | | UnCorrelation | 0.0921 | 0.0918 | Close enough |
Feel free to provide more results, and I'll continue updating the table!
### Testing 259,200 bytes from test2.bin (/dev/random) ### Testing 259,200 bytes from test2.bin (/dev/random)
Test | Java ent3000 | Rust onod3000 | Note | Test | Java ent3000 | Rust onod3000 | Note |
------------- | ------------ | ------------- | ---- | ------------- | ------------ | ------------- | ------------- |
Shannon | N/A | 0.6912 | | Shannon | N/A | 0.6853 | |
Monobit | 0.0619 | 0.0619 | | Monobit | 0.2103 | 0.2103 | |
ChiBit | 0.5461 | 0.5461 | | ChiBit | 0.5590 | 0.5590 | |
ChiByte | 0.4672 | 0.4672 | | ChiByte | 0.3748 | 0.3748 | |
MeanByte | 0.0500 | 0.0500 | | MeanByte | 0.0183 | 0.0183 | |
Compression | N/A | 0.9308 | | Compression | N/A | 0.9308 | |
KS | 0.9229 | 0.0028 | MISMATCH | KS | 0.2900 | 0.0178 | Mismatch |
Pi | 0.2828 | 0.3953 | MISMATCH | Pi | 0.9765 | 0.9765 | |
Shells | 0.4604 | 0.3163 | MISMATCH | Shells | 0.7235 | 0.7235 | |
Gaps | 0.3854 | 0.3854 | | Gaps | 0.6425 | 0.6425 | |
Avalanche | 0.9974 | 0.9974 | | Avalanche | 0.9995 | 0.9995 | |
Runs | 0.5873 | 0.5873 | | Runs | 0.3310 | 0.3310 | |
RunUps | 0.0402 | 0.0402 | | RunUps | 0.9393 | 0.9393 | |
Prediction | 0.7409 | 0.3569 | MISMATCH | Prediction | 0.1519 | 0.1519 | |
UnCorrelation | 0.0949 | 0.0949 | | UnCorrelation | 0.4497 | 0.4497 | |
### Testing 259,200 bytes from test3.bin (Hardware QRNG) ### Testing 259,200 bytes from test3.bin (Hardware QRNG)
@ -143,23 +145,25 @@ ChiBit | 0.8428 | 0.8428 | |
ChiByte | 0.5785 | 0.5785 | | ChiByte | 0.5785 | 0.5785 | |
MeanByte | 0.9601 | 0.9601 | | MeanByte | 0.9601 | 0.9601 | |
Compression | N/A | 0.9308 | | Compression | N/A | 0.9308 | |
KS | 0.7395 | 0.0247 | MISMATCH | KS | 0.7395 | 0.0441 | Mismatch |
Pi | 0.2806 | 0.0189 | MISMATCH | Pi | 0.2806 | 0.2806 | |
Shells | 0.7711 | 0.3254 | MISMATCH | Shells | 0.7711 | 0.7711 | |
Gaps | 0.1937 | 0.1937 | | Gaps | 0.1937 | 0.1937 | |
Avalanche | 0.9932 | 0.9932 | | Avalanche | 0.9932 | 0.9932 | |
Runs | 0.2500 | 0.2500 | | Runs | 0.2500 | 0.2500 | |
RunUps | 0.9710 | 0.9710 | | RunUps | 0.9710 | 0.9710 | |
Prediction | 0.7173 | 0.7399 | MISMATCH | Prediction | 0.7173 | 0.7173 | |
UnCorrelation | 0.4674 | 0.4674 | | UnCorrelation | 0.4674 | 0.4674 | |
The port is work in progress and effort will be put into rigorously re-evaluating the implementations. Although the implementation of the randomness tests closely follows the original Java logic, minor differences in the results arise in some of the tests. 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.
Java emphasizes predictability and portability, enforcing strict IEEE 754 behavior across platforms. Rust prioritizes performance and flexibility, allowing platform-specific optimizations that may deviate slightly from strict IEEE semantics.
These differences are generally negligible for practical purposes and do not affect the overall functionality or These differences are generally negligible for practical purposes and do not affect the overall functionality or
statistical significance of the test. 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.
## License ## License
This project is licensed under the MIT License. See the [LICENSE](LICENSE) file for details. This project is licensed under the MIT License. See the [LICENSE](LICENSE) file for details.

30
src/chisquaretest.rs Normal file
View file

@ -0,0 +1,30 @@
use statrs::distribution::{ChiSquared, ContinuousCDF};
pub fn chi_square_test(observed: &[u64], expected: &[f64]) -> f64 {
// Preconditions
if observed.len() != expected.len() || observed.len() < 2 {
panic!("Observed and expected arrays must have the same length and length >= 2.");
}
if expected.iter().any(|&e| e <= 0.0) {
panic!("Expected array must contain only strictly positive values.");
}
// Rescale expected array if necessary
let sum_observed: f64 = observed.iter().map(|&o| o as f64).sum();
let sum_expected: f64 = expected.iter().sum();
let rescaled_expected: Vec<f64> = expected.iter().map(|&e| e * sum_observed / sum_expected).collect();
// Calculate chi-squared statistic
let chi_squared_stat: f64 = observed
.iter()
.zip(rescaled_expected.iter())
.map(|(&o, &e)| (o as f64 - e).powi(2) / e)
.sum();
// Perform chi-squared test
let degrees_of_freedom = observed.len() as f64 - 1.0;
let chi_squared_dist = ChiSquared::new(degrees_of_freedom).expect("Failed to create ChiSquared distribution");
let p_value = 1.0 - chi_squared_dist.cdf(chi_squared_stat);
p_value
}

2
src/lib.rs
View file

@ -1,6 +1,8 @@
pub struct Onod; pub struct Onod;
mod uniformity; mod uniformity;
pub mod chisquaretest;
pub mod well19937c;
pub mod ffi; pub mod ffi;
#[cfg(feature = "python")] #[cfg(feature = "python")]
pub mod python; pub mod python;

69
src/uniformity/ks.rs
View file

@ -1,47 +1,76 @@
use kolmogorov_smirnov::test_f64; use kolmogorov_smirnov::test_f64;
use std::time::{SystemTime, UNIX_EPOCH};
use statrs::distribution::ContinuousCDF;
use crate::Onod; use crate::Onod;
use crate::well19937c::Well19937c;
impl Onod { impl Onod {
/// KS randomness test /// KS randomness test
/// Performs the Kolmogorov-Smirnov test to evaluate uniformity of data distribution and returns a p-value. /// Performs the Kolmogorov-Smirnov test to evaluate the uniformity of data distribution
// Note: This implementation of the Kolmogorov-Smirnov (KS) test differs slightly from the /// and returns the test statistic (D-statistic), z-score, and p-value.
// Java implementation due to differences in library behavior and floating-point handling.
// The Java implementation (Apache Commons Math) adds random jitter to handle ties in small
// datasets, uses strict inequality for small sample sizes, and applies specific precision
// rules based on the IEEE 754 standard for `double`. The Rust implementation, using the
// `kolmogorov_smirnov` crate, does not add jitter or handle ties in the same way, and
// adheres to the crate's internal handling of floating-point comparisons. These differences
// may result in slight variations in p-values or KS statistics between the two versions.
pub fn ks(samples: &[u8]) -> (f64, f64, f64) { pub fn ks(samples: &[u8]) -> (f64, f64, f64) {
if samples.is_empty() { if samples.is_empty() {
return (-1.0, 0.0, 1.0); // empty data return (-1.0, 0.0, 1.0); // Invalid input
} }
// Normalize the input samples to [0, 1] range // Normalize the input samples to [0, 1) range
let normalized_samples: Vec<f64> = samples.iter().map(|&x| x as f64 / 255.0).collect(); let mut normalized_samples: Vec<f64> = samples.iter().map(|&x| x as f64 / 255.0).collect();
// Generate a uniform distribution for comparison // Generate a uniform distribution for comparison
let uniform_distribution: Vec<f64> = (0..normalized_samples.len()) let mut uniform_distribution: Vec<f64> = Self::generate_uniform_distribution(normalized_samples.len(), Self::get_timestamp_seed());
.map(|i| i as f64 / (normalized_samples.len() as f64 - 1.0))
.collect(); // 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);
}
// Perform the Kolmogorov-Smirnov test // Perform the Kolmogorov-Smirnov test
let confidence = 0.01; // Significance level let confidence = 0.05; // Significance level
let result = test_f64(&normalized_samples, &uniform_distribution, confidence); let result = test_f64(&normalized_samples, &uniform_distribution, confidence);
// Extract the KS statistic (D-statistic) // Extract the KS statistic (D-statistic)
let ks_statistic = result.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);
}
// Calculate the z-score // Calculate the z-score
let sample_size = normalized_samples.len() as f64; let sample_size = normalized_samples.len() as f64;
let z_score = ks_statistic * sample_size.sqrt(); let z_score = ks_statistic * sample_size.sqrt();
// Extract the p-value // Calculate the p-value
let p_value = 1.0 - result.reject_probability; let p_value = 2.0 * (1.0 - statrs::distribution::Normal::new(0.0, 1.0).unwrap().cdf(z_score.abs()));
(ks_statistic, z_score, p_value) (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<f64> {
let mut rng = Well19937c::new(seed);
(0..len).map(|_| rng.next_f64()).collect()
}
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())
}
} }

70
src/uniformity/pi.rs
View file

@ -1,63 +1,37 @@
/// Pi randomness test
/// Uses a Monte Carlo simulation to estimate randomness by calculating the approximation of Pi.
use statrs::distribution::{Normal, ContinuousCDF}; use statrs::distribution::{Normal, ContinuousCDF};
use crate::Onod; use crate::Onod;
impl Onod { impl Onod {
/// Pi randomness test using nalgebra for vectorized operations.
/// Pi randomness test
/// Uses a Monte Carlo simulation to estimate randomness by calculating the approximation of Pi.
/// This implementation of the Pi randomness test closely follows the logic of the original Java implementation.
/// However, minor differences in the results may arise due to the following reasons:
///
/// 1. **Floating-Point Precision**:
/// Rust and Java both use 64-bit floating-point numbers (`double` in Java, `f64` in Rust), but slight differences
/// in their implementations (e.g., rounding modes, intermediate representations) can lead to small deviations.
///
/// 2. **Math Libraries**:
/// Java uses Apache Commons Math for statistical computations, which may implement certain calculations
/// (e.g., Z-scores and normal distribution CDFs) differently compared to the `statrs` crate used in Rust.
///
/// 3. **Bit Accuracy**:
/// The Java implementation notes the significance of bit accuracy in floating-point computations,
/// as defined in the IEEE 754 standard. Differences in handling edge cases (e.g., subnormal values,
/// precision limits) could lead to slight variations.
///
/// These differences are generally negligible for practical purposes and do not affect the overall functionality or
/// statistical significance of the test.
pub fn pi(samples: &[u8]) -> (f64, f64, f64) { pub fn pi(samples: &[u8]) -> (f64, f64, f64) {
if samples.len() < 4 {
return (-1.0, 0.0, 1.0); // Not enough data
}
if samples.is_empty() { let normalized_samples: Vec<f32> = get_floats(samples);
if normalized_samples.is_empty() {
return (-1.0, 0.0, 1.0); return (-1.0, 0.0, 1.0);
} }
// Normalize samples to [0.0, 1.0)
let normalized_samples: Vec<f64> = samples.iter().map(|&x| x as f64 / 255.0).collect();
// Initialize variables for summary statistics
let mut sum_y = 0.0; let mut sum_y = 0.0;
let mut count = 0.0; let count = normalized_samples.len() as f64;
// Compute y-values (sqrt(1 - x^2)) and update summary statistics
for &x in &normalized_samples { for &x in &normalized_samples {
let y = (1.0 - x * x).sqrt(); let y = (1.0 - x.powi(2)).sqrt();
sum_y += y; sum_y += y as f64;
count += 1.0;
} }
// Calculate mean of y-values
let mean_y = sum_y / count; let mean_y = sum_y / count;
// Calculate the test statistic
let test_statistic = 4.0 * mean_y; let test_statistic = 4.0 * mean_y;
let variance = compute_variance(count);
// Calculate variance and standard deviation
let variance = (16.0 / count) * ((2.0 / 3.0) - (std::f64::consts::PI / 4.0).powi(2));
let std_dev = variance.sqrt(); let std_dev = variance.sqrt();
// Calculate Z-score
let z_score = (test_statistic - std::f64::consts::PI) / std_dev; let z_score = (test_statistic - std::f64::consts::PI) / std_dev;
// Use normal distribution to calculate p-value
let normal_dist = Normal::new(0.0, 1.0).expect("Failed to create Normal distribution"); let normal_dist = Normal::new(0.0, 1.0).expect("Failed to create Normal distribution");
let p_value = 2.0 * (1.0 - normal_dist.cdf(z_score.abs())); let p_value = 2.0 * (1.0 - normal_dist.cdf(z_score.abs()));
@ -65,3 +39,19 @@ impl Onod {
} }
} }
fn get_floats(samples: &[u8]) -> Vec<f32> {
let mut floats = Vec::new();
for chunk in samples.chunks_exact(4) {
let int_val = i32::from_be_bytes([chunk[0], chunk[1], chunk[2], chunk[3]]);
let unsigned_val = (int_val as u32) >> 1; // Discard sign bit
let normalized = unsigned_val as f32 / i32::MAX as f32;
floats.push(normalized);
}
floats
}
fn compute_variance(n: f64) -> f64 {
let term = (2.0 / 3.0) - (std::f64::consts::PI / 4.0).powi(2);
(16.0 / n) * term
}

71
src/uniformity/prediction.rs
View file

@ -1,51 +1,66 @@
use statrs::distribution::{Normal, ContinuousCDF}; use statrs::distribution::{ChiSquared, ContinuousCDF};
/*
* Blatantly copied from David Sexton's battery.
*
* An algorithm is used to predict the value of each byte of the sequence from
* the beginning of the sequence to the end. In a random sequence the
* probability of success of any such algorithm is 1/256. The number of successes
* is counted. A chi-squared statistic is calculated. The degrees-of-freedom is 1.
* The algorithm is as follows: the next byte is predicted to be equal to all the
* previous bytes bitwise XORed together.
*/
use crate::Onod; use crate::Onod;
impl Onod { impl Onod {
/// Prediction randomness test /// Prediction randomness test
/// Evaluates the predictability of the next bit based on current data and returns a p-value. /// Evaluates the predictability of the next byte based on XORing the previous bytes
/// and returns the total predictions, z-score, and p-value.
pub fn prediction(samples: &[u8]) -> (f64, f64, f64) { pub fn prediction(samples: &[u8]) -> (f64, f64, f64) {
if samples.len() < 3 {
if samples.is_empty() { return (-1.0, 0.0, 1.0); // Not enough data for meaningful calculation
return (-1.0, 0.0, 1.0);
} }
let mut correct_predictions = 0; let mut correct_predictions = 0;
let mut total_predictions = 0; let mut total_predictions = 0;
for window in samples.windows(2) { let mut prediction = samples[0]; // Start with the first byte
if let [current, next] = window { for i in 2..samples.len() {
let predicted = if current & 0x01 == 0 { 0 } else { 1 }; // Predict next bit based on LSB prediction ^= samples[i - 1]; // XOR all preceding bytes
let actual = next & 0x01; // Check LSB of the next byte
if predicted == actual { if prediction == samples[i] {
correct_predictions += 1; correct_predictions += 1;
} }
total_predictions += 1; total_predictions += 1;
} }
}
if total_predictions == 0 { // Calculate expected and observed frequencies
return (-1.0, 0.0, 1.0); // No predictions possible let expected = vec![
} (1.0 / 256.0) * samples.len() as f64, // Probability of correct prediction
(255.0 / 256.0) * samples.len() as f64, // Probability of incorrect prediction
];
let observed = vec![
correct_predictions as f64, // Actual correct predictions
(samples.len() - correct_predictions) as f64, // Actual incorrect predictions
];
// Calculate observed proportion of correct predictions // Calculate chi-squared statistic
let observed_proportion = correct_predictions as f64 / total_predictions as f64; let chi_squared_stat: f64 = observed
.iter()
.zip(expected.iter())
.map(|(o, e)| (o - e).powi(2) / e)
.sum();
// Expected proportion for random data // Use Chi-Squared distribution to calculate p-value
let expected_proportion = 0.5; let chi_squared_dist = ChiSquared::new(1.0).unwrap(); // Degrees of freedom = 1
let std_dev = (0.5 * 0.5 / total_predictions as f64).sqrt(); // Standard deviation for a binomial distribution let p_value = 1.0 - chi_squared_dist.cdf(chi_squared_stat);
// Calculate the z-score // Calculate z-score (optional, for diagnostics)
let z_score = (observed_proportion - expected_proportion) / std_dev; let mean = 1.0; // Mean of chi-squared distribution
let std_dev = (2.0 as f64).sqrt(); // Standard deviation of chi-squared distribution
let z_score = (chi_squared_stat - mean) / std_dev;
// Use normal distribution to calculate p-value (total_predictions as f64, z_score, p_value)
let normal_dist = Normal::new(0.0, 1.0).expect("Failed to create Normal distribution");
let p_value = 2.0 * (1.0 - normal_dist.cdf(z_score.abs()));
(observed_proportion, z_score, p_value)
} }
} }

103
src/uniformity/shells.rs
View file

@ -1,87 +1,84 @@
use statrs::distribution::{ChiSquared, ContinuousCDF};
use crate::Onod; use crate::Onod;
use crate::chisquaretest::chi_square_test;
impl Onod { impl Onod {
/// Shells randomness test /// Shells randomness test
/// Evaluates the uniformity of distances between identical byte values and returns a p-value. /// Evaluates the uniformity of distances in a 3D sphere and returns the chi-squared statistic, z-score, and p-value.
pub fn shells(input: &[u8]) -> (f64, f64, f64) { pub fn shells(input: &[u8]) -> (f64, f64, f64) {
// Define shell radii (precomputed to ensure equal volumes) // Define shell radii (precomputed to ensure equal volumes)
const SHELL_RADII: [f64; 35] = [ const SHELL_RADII: [f64; 35] = [
1., 0.990384019787941, 0.980577593308067, 0.970571001281035, 0.960353705642329, 1., 0.990384019787941, 0.980577593308067,
0.949914251592996, 0.939240154232372, 0.928317766722556, 0.91713212619864, 0.970571001281035, 0.960353705642329, 0.949914251592996,
0.905666772691187, 0.893903535096568, 0.881822276616739, 0.869400589952457, 0.939240154232372, 0.928317766722556, 0.91713212619864,
0.856613429672063, 0.843432665301749, 0.829826533366243, 0.815758959214771, 0.905666772691187, 0.893903535096568, 0.881822276616739,
0.801188709029197, 0.786068317431936, 0.770342714221672, 0.753947441129154, 0.869400589952457, 0.856613429672063, 0.843432665301749,
0.736806299728077, 0.718828193851318, 0.699902804775202, 0.67989452969576, 0.829826533366243, 0.815758959214771, 0.801188709029197,
0.65863375600835, 0.635903899768996, 0.61142141746576, 0.584803547642573, 0.786068317431936, 0.770342714221672, 0.753947441129154,
0.555513224287824, 0.52275795857471, 0.485285500640517, 0.440911138308369, 0.736806299728077, 0.718828193851318, 0.699902804775202,
0.385171357110836, 0.30571070873288 0.67989452969576, 0.65863375600835, 0.635903899768996,
0.61142141746576, 0.584803547642573, 0.555513224287824,
0.52275795857471, 0.485285500640517, 0.440911138308369,
0.385171357110836, 0.30571070873288,
]; ];
let samples = convert_to_3d_points(input); let samples = convert_to_3d_points(&input);
// if samples.len() < 25000 {
// // eprintln!("---------------------------------------------------------------");
// // eprintln!("ERROR: Shells test requires at least 25,000 points for statistical validity. Skipping.");
// // eprintln!("---------------------------------------------------------------");
// return (-1.0, 0.0, 1.0); // Skip the test for small datasets
// }
let sphere_radius = SHELL_RADII[0]; let sphere_radius = SHELL_RADII[0];
let no_shells = SHELL_RADII.len(); let num_shells = SHELL_RADII.len();
// Calculate sphere and cube volume proportions // Calculate sphere and cube proportions
let cube_side = 2.0 * sphere_radius; let cube_side = 2.0 * sphere_radius;
let cube_volume = cube_side.powi(3); let cube_volume = cube_side.powi(3);
let sphere_volume = (4.0 / 3.0) * std::f64::consts::PI * sphere_radius.powi(3); let sphere_volume = (4.0 / 3.0) * std::f64::consts::PI * sphere_radius.powi(3);
let sphere_proportion = sphere_volume / cube_volume; // Theoretical value: π/6 let sphere_proportion = sphere_volume / cube_volume;
let no_points = samples.len() as f64; let num_points = samples.len() as f64;
let no_points_per_shell = sphere_proportion * no_points / no_shells as f64; let num_points_per_shell = sphere_proportion * num_points / num_shells as f64;
let mut observed = vec![0u64; no_shells]; // Initialize observed and expected frequencies
let expected: Vec<f64> = vec![no_points_per_shell; no_shells]; let mut observed = vec![0u64; num_shells];
let expected: Vec<f64> = vec![num_points_per_shell; num_shells];
// Assign points to shells
for (x, y, z) in samples { for (x, y, z) in samples {
// Compute radius from origin
let radius = (x.powi(2) + y.powi(2) + z.powi(2)).sqrt(); let radius = (x.powi(2) + y.powi(2) + z.powi(2)).sqrt();
// Ignore points outside the sphere // Skip points outside the sphere
if radius > sphere_radius { if radius > sphere_radius {
continue; continue;
} }
// Assign to the correct shell for j in 1..num_shells {
for j in 1..SHELL_RADII.len() {
if radius > SHELL_RADII[j] { if radius > SHELL_RADII[j] {
observed[j - 1] += 1; observed[j - 1] += 1;
break; break;
} }
} }
if radius < SHELL_RADII[no_shells - 1] { // Assign to the last shell if radius <= SHELL_RADII[num_shells - 1]
observed[no_shells - 1] += 1; if radius < SHELL_RADII[num_shells - 1] {
} observed[num_shells - 1] += 1;
} }
// Perform Chi-Square Test }
let chi_squared_stat: f64 = observed.iter()
// Calculate chi-squared statistic
let chi_squared_stat: f64 = observed
.iter()
.zip(expected.iter()) .zip(expected.iter())
.map(|(&o, &e)| (o as f64 - e).powi(2) / e) .map(|(&o, &e)| (o as f64 - e).powi(2) / e)
.sum(); .sum();
let degrees_of_freedom = no_shells as f64 - 1.0; // Perform chi-squared test
let chi_squared_dist = ChiSquared::new(degrees_of_freedom).expect("Failed to create ChiSquared distribution"); let degrees_of_freedom = num_shells as f64 - 1.0;
let p_value = 1.0 - chi_squared_dist.cdf(chi_squared_stat); let p_value = chi_square_test(&observed, &expected);
// Z-score calculation (standardization of the chi-squared statistic) // Calculate z-score
let mean = degrees_of_freedom; // Mean of the chi-squared distribution let mean = degrees_of_freedom;
let std_dev = (2.0 * degrees_of_freedom).sqrt(); // Standard deviation of the chi-squared distribution let std_dev = (2.0 * degrees_of_freedom).sqrt();
let z_score = (chi_squared_stat - mean) / std_dev; let z_score = (chi_squared_stat - mean) / std_dev;
// Return the results
(chi_squared_stat, z_score, p_value) (chi_squared_stat, z_score, p_value)
} }
} }
@ -89,13 +86,17 @@ impl Onod {
fn convert_to_3d_points(data: &[u8]) -> Vec<(f64, f64, f64)> { fn convert_to_3d_points(data: &[u8]) -> Vec<(f64, f64, f64)> {
let mut points = Vec::new(); let mut points = Vec::new();
for chunk in data.chunks(3) { for chunk in data.chunks(12) {
if chunk.len() == 3 { if chunk.len() == 12 {
// Normalize the bytes to [0.0, 1.0) range let x = u32::from_be_bytes([chunk[0], chunk[1], chunk[2], chunk[3]]) >> 1;
let x = chunk[0] as f64 / 255.0; let y = u32::from_be_bytes([chunk[4], chunk[5], chunk[6], chunk[7]]) >> 1;
let y = chunk[1] as f64 / 255.0; let z = u32::from_be_bytes([chunk[8], chunk[9], chunk[10], chunk[11]]) >> 1;
let z = chunk[2] as f64 / 255.0;
points.push((x, y, z)); points.push((
x as f64 / (i32::MAX as f64),
y as f64 / (i32::MAX as f64),
z as f64 / (i32::MAX as f64),
));
} }
} }

62
src/well19937c.rs Normal file
View file

@ -0,0 +1,62 @@
pub struct Well19937c {
state: [u32; 624],
index: usize,
}
impl Well19937c {
/// Creates a new instance of Well19937c with a given seed.
pub fn new(seed: u32) -> Self {
let mut state = [0u32; 624];
state[0] = seed;
for i in 1..624 {
state[i] = 1812433253u32
.wrapping_mul(state[i - 1] ^ (state[i - 1] >> 30))
.wrapping_add(i as u32);
}
Well19937c { state, index: 0 }
}
/// Updates the internal state.
fn twist(&mut self) {
const M: usize = 397;
const MATRIX_A: u32 = 0x9908b0df; // Constant matrix A
const UPPER_MASK: u32 = 0x80000000; // Most significant w-r bits
const LOWER_MASK: u32 = 0x7fffffff; // Least significant r bits
for i in 0..624 {
let x = (self.state[i] & UPPER_MASK) + (self.state[(i + 1) % 624] & LOWER_MASK);
let mut x_a = x >> 1;
if x % 2 != 0 {
x_a ^= MATRIX_A;
}
self.state[i] = self.state[(i + M) % 624] ^ x_a;
}
self.index = 0;
}
/// Generates the next random number in the sequence.
pub fn next_u32(&mut self) -> u32 {
if self.index == 0 {
self.twist();
}
let mut y = self.state[self.index];
self.index = (self.index + 1) % 624;
// Matsumoto-Kurita tempering
y ^= (y << 7) & 0xe46e1700;
y ^= (y << 15) & 0x9b868000;
y
}
/// Generates the next random `f64` in [0, 1).
pub fn next_f64(&mut self) -> f64 {
self.next_u32() as f64 / u32::MAX as f64
}
}