Skip to main content

sc_neurocore_engine/
cordiv.rs

1// SPDX-License-Identifier: AGPL-3.0-or-later | Commercial license available
2// © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
3// © Code 2020–2026 Miroslav Šotek. All rights reserved.
4// ORCID: 0009-0009-3560-0851
5// Contact: www.anulum.li | protoscience@anulum.li
6// SC-NeuroCore — CORDIV stochastic division (Li et al. 2014)
7
8//! CORDIV: stochastic computing divider.
9//!
10//! Sequential circuit: P(z=1) ≈ P(x=1) / P(y=1) when P(x) <= P(y).
11//!
12//! Truth table per bit:
13//!   x=1         → z = 1
14//!   x=0, y=1    → z = 0
15//!   x=0, y=0    → z = z_prev (hold)
16
17/// Compute SC division on unpacked bitstreams.
18/// Returns quotient bitstream of same length.
19pub fn cordiv(numerator: &[u8], denominator: &[u8]) -> Vec<u8> {
20    let len = numerator.len().min(denominator.len());
21    let mut out = vec![0u8; len];
22    let mut prev = 0u8;
23    for t in 0..len {
24        if numerator[t] == 1 {
25            out[t] = 1;
26        } else if denominator[t] == 1 {
27            out[t] = 0;
28        } else {
29            out[t] = prev;
30        }
31        prev = out[t];
32    }
33    out
34}
35
36/// Compute SC division on packed u64 bitstreams (bit-serial).
37/// Each u64 word contains 64 sequential bits.
38pub fn cordiv_packed(num_packed: &[u64], den_packed: &[u64], length: usize) -> Vec<u64> {
39    let n_words = num_packed.len().min(den_packed.len());
40    let mut out = vec![0u64; n_words];
41    let mut prev_bit = 0u8;
42
43    let total_bits = length.min(n_words * 64);
44    for bit_idx in 0..total_bits {
45        let word = bit_idx / 64;
46        let bit = bit_idx % 64;
47        let x = ((num_packed[word] >> bit) & 1) as u8;
48        let y = ((den_packed[word] >> bit) & 1) as u8;
49
50        let z = if x == 1 {
51            1u8
52        } else if y == 1 {
53            0u8
54        } else {
55            prev_bit
56        };
57
58        if z == 1 {
59            out[word] |= 1u64 << bit;
60        }
61        prev_bit = z;
62    }
63    out
64}
65
66/// Compute adaptive bitstream length via Hoeffding bound.
67/// Returns minimum length (rounded up to power of 2) for target precision.
68pub fn adaptive_length_hoeffding(epsilon: f64, confidence: f64) -> usize {
69    let delta = 1.0 - confidence;
70    if delta <= 0.0 || epsilon <= 0.0 {
71        return 65536;
72    }
73    let l = -(delta / 2.0).ln() / (2.0 * epsilon * epsilon);
74    let l_int = l.ceil() as usize;
75    let l_int = l_int.max(64);
76    // Round up to power of 2
77    let mut p = 1usize;
78    while p < l_int {
79        p *= 2;
80    }
81    p.min(65536)
82}
83
84#[cfg(test)]
85mod tests {
86    use super::*;
87
88    #[test]
89    fn test_cordiv_half_div_one() {
90        // x = alternating 0,1 (p=0.5), y = all 1s (p=1.0)
91        // result should be ~0.5
92        let x: Vec<u8> = (0..1000).map(|i| if i % 2 == 0 { 1 } else { 0 }).collect();
93        let y = vec![1u8; 1000];
94        let z = cordiv(&x, &y);
95        let p: f64 = z.iter().map(|&b| b as f64).sum::<f64>() / z.len() as f64;
96        assert!((p - 0.5).abs() < 0.05, "Expected ~0.5, got {p}");
97    }
98
99    #[test]
100    fn test_cordiv_zero_numerator() {
101        let x = vec![0u8; 100];
102        let y = vec![1u8; 100];
103        let z = cordiv(&x, &y);
104        assert!(z.iter().all(|&b| b == 0));
105    }
106
107    #[test]
108    fn test_cordiv_identity() {
109        // x/x should be ~1.0
110        let x: Vec<u8> = (0..1000).map(|i| if i % 3 != 0 { 1 } else { 0 }).collect();
111        let z = cordiv(&x, &x);
112        let p: f64 = z.iter().map(|&b| b as f64).sum::<f64>() / z.len() as f64;
113        assert!(p > 0.9, "Expected ~1.0, got {p}");
114    }
115
116    #[test]
117    fn test_adaptive_length_hoeffding() {
118        let l = adaptive_length_hoeffding(0.01, 0.95);
119        assert!(l >= 64);
120        assert!(l <= 65536);
121        assert!(l & (l - 1) == 0); // power of 2
122    }
123
124    #[test]
125    fn test_adaptive_tighter_needs_more() {
126        let l_coarse = adaptive_length_hoeffding(0.1, 0.95);
127        let l_fine = adaptive_length_hoeffding(0.01, 0.95);
128        assert!(l_fine > l_coarse);
129    }
130}