pi_compression/
lib.rs

1#![feature(let_chains)]
2// This is a port of the BBP formula, converted from C code provided by
3// David H. Bailey.
4//
5// Info: https://www.experimentalmath.info/bbp-codes/
6// Direct C link: http://www.experimentalmath.info/bbp-codes/piqpr8.c
7
8use lazy_static::lazy_static;
9
10const EPSILON: f64 = 1e-17;
11const SIXTEEN: f64 = 16.0;
12
13#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
14pub struct Location {
15    pub offset: usize,
16    pub length: usize,
17}
18
19pub struct PiIterator {
20    offset: usize,
21    digits: Vec<u8>,
22}
23
24impl PiIterator {
25    pub fn new() -> PiIterator {
26        PiIterator {
27            offset: 0,
28            digits: vec![],
29        }
30    }
31
32    pub fn from(offset: usize) -> PiIterator {
33        PiIterator {
34            offset,
35            digits: vec![],
36        }
37    }
38
39    pub fn refill(&mut self) {
40        self.digits = pi_digits(self.offset);
41        self.digits.reverse();
42        self.offset += self.digits.len();
43    }
44}
45
46impl Iterator for PiIterator {
47    type Item = u8;
48
49    fn next(&mut self) -> Option<Self::Item> {
50        if self.digits.is_empty() {
51            self.refill();
52        }
53        self.digits.pop()
54    }
55}
56
57pub fn nibbles(msg: String) -> Vec<u8> {
58    msg.bytes().flat_map(|b| [b / 16, b & 0x0F]).collect()
59}
60
61pub fn compress(msg: &[u8], limit: usize) -> Vec<Location> {
62    let mut locs = vec![];
63
64    let mut index = 0;
65    while index < msg.len() {
66        let m = find_pi_match(&msg[index..], limit).expect("should find some match");
67        index += m.length;
68        locs.push(m);
69    }
70
71    locs
72}
73
74pub fn find_pi_match(msg: &[u8], limit: usize) -> Option<Location> {
75    let mut best_match: Option<Location> = None;
76    let mut offset = 0;
77
78    'outer: while offset < limit {
79        let mut pi_digits = PiIterator::from(offset);
80
81        while let Some(b) = pi_digits.next()
82            && b != msg[0]
83        {
84            offset += 1;
85            if offset >= limit {
86                break 'outer;
87            }
88        }
89
90        let length = pi_digits
91            .zip(msg.iter().skip(1))
92            .take_while(|(a, &b)| *a == b)
93            .count() + 1;
94
95        if length == msg.len() {
96            return Some(Location { offset, length });
97        } else if let Some(m) = &best_match {
98            if length > m.length {
99                best_match = Some(Location { offset, length });
100            }
101        } else {
102            best_match = Some(Location { offset, length });
103        }
104
105        offset += 1;
106    }
107
108    best_match
109}
110
111pub fn pi_digits(id: usize) -> Vec<u8> {
112    let s1 = series(1, id as i64);
113    let s2 = series(4, id as i64);
114    let s3 = series(5, id as i64);
115    let s4 = series(6, id as i64);
116    let pid = 4. * s1 - 2. * s2 - s3 - s4;
117    let pid = pid - pid.floor() + 1.;
118
119    to_hex(pid)
120}
121
122pub fn to_hex(x: f64) -> Vec<u8> {
123    let mut chx = vec![0; 8];
124
125    let mut y = x.abs();
126
127    for i in 0..8 {
128        y = 16.0 * (y - y.floor());
129        chx[i] = y as u8;
130    }
131
132    chx
133}
134
135pub fn series(m: i64, id: i64) -> f64 {
136    let mut s: f64 = 0.0;
137
138    for k in 0..id {
139        let ak = (8 * k + m) as f64;
140        let p = (id - k) as f64;
141        let t = expm(p, ak);
142        s = s + t / ak;
143        s = s - f64::floor(s);
144    }
145
146    for k in id..=(id + 100) {
147        let ak = (8 * k + m) as f64;
148        let t = SIXTEEN.powf((id - k) as f64) / ak;
149        if t < EPSILON {
150            break;
151        }
152
153        s = s + t;
154        s = s - s.floor();
155    }
156
157    s
158}
159
160const NUM_POWERS: usize = 25;
161lazy_static! {
162    static ref TWO_TO_THE: [f64; NUM_POWERS] = {
163        let mut t = [0.0; NUM_POWERS];
164
165        t[0] = 1.0;
166        for i in 1..NUM_POWERS {
167            t[i] = t[i - 1] * 2.;
168        }
169
170        t
171    };
172}
173
174pub fn expm(p: f64, ak: f64) -> f64 {
175    if ak == 1.0 {
176        return 0.0;
177    }
178
179    let mut i = 0;
180    for idx in 1..NUM_POWERS {
181        if TWO_TO_THE[idx] > p {
182            i = idx;
183            break;
184        }
185    }
186    let mut pt = TWO_TO_THE[i - 1];
187
188    let mut p1 = p;
189    let mut r: f64 = 1.0;
190
191    for _ in 1..=i {
192        if p1 >= pt {
193            r = 16.0 * r;
194            r = r - (r / ak).floor() * ak;
195            p1 = p1 - pt;
196        }
197        pt = 0.5 * pt;
198        if pt >= 1.0 {
199            r = r * r;
200            r = r - (r / ak).floor() * ak;
201        }
202    }
203
204    r
205}
pi_compression/lib.rs

pi_compression/
lib.rs
