Skip to main content

sc_neurocore_engine/
cordiv.rs

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