sva_samples/
grid.rs

1// Concern: the sample grid: a sample's instant, the sample an instant or edge lands on, spans of samples | Non-concern: what a sample holds | IO: (n or t) -> t or n
2
3fn decimal(x: f64) -> Option<(i128, i128)> {
4    const LIMIT: i128 = 1 << 120;
5    let text = format!("{x:e}");
6    let (mantissa, exp) = text.split_once('e')?;
7    let (whole, frac) = mantissa.split_once('.').unwrap_or((mantissa, ""));
8    let digits: i128 = format!("{whole}{frac}").parse().ok()?;
9    let shift = exp.parse::<i32>().ok()? - frac.len() as i32;
10    let ten = 10i128
11        .checked_pow(shift.unsigned_abs())
12        .filter(|p| *p < LIMIT)?;
13    match shift >= 0 {
14        true => Some((digits.checked_mul(ten)?, 1)),
15        false => Some((digits, ten)),
16    }
17}
18
19/// No pull, stored chunk or machine block crosses a multiple.
20pub const BLOCK: i64 = 1 << 12;
21
22pub fn block_end(n: i64) -> i64 {
23    (n.div_euclid(BLOCK) + 1).saturating_mul(BLOCK)
24}
25
26/// Samples `[start, end)` of the grid, whose sample 0 is t = 0.
27#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
28pub struct Extent {
29    pub start: i64,
30    pub end: i64,
31}
32
33impl Extent {
34    /// `i64::MIN` and `i64::MAX` stand for no edge.
35    pub const EVERYWHERE: Extent = Extent {
36        start: i64::MIN,
37        end: i64::MAX,
38    };
39
40    pub const NOWHERE: Extent = Extent { start: 0, end: 0 };
41
42    pub fn new(start: i64, end: i64) -> Extent {
43        assert!(start <= end, "an extent [{start}, {end}) runs backwards");
44        Extent { start, end }
45    }
46
47    pub fn from(start: i64) -> Extent {
48        Extent::new(start, i64::MAX)
49    }
50
51    pub fn is_bounded(&self) -> bool {
52        self.start != i64::MIN && self.end != i64::MAX
53    }
54
55    pub fn contains(&self, n: i64) -> bool {
56        self.start <= n && n < self.end
57    }
58
59    pub fn intersect(self, other: Extent) -> Extent {
60        let start = self.start.max(other.start);
61        let end = self.end.min(other.end);
62        match start < end {
63            true => Extent { start, end },
64            false => Extent::NOWHERE,
65        }
66    }
67
68    pub fn hull(self, other: Extent) -> Extent {
69        match (self.is_empty(), other.is_empty()) {
70            (true, _) => other,
71            (_, true) => self,
72            _ => Extent {
73                start: self.start.min(other.start),
74                end: self.end.max(other.end),
75            },
76        }
77    }
78
79    /// Every sample moved `by` later; an unbounded edge stays unbounded.
80    pub fn shifted(self, by: i64) -> Extent {
81        if self.is_empty() {
82            return self;
83        }
84        let edge = |n: i64| match n {
85            i64::MIN | i64::MAX => n,
86            n => n.saturating_add(by).clamp(i64::MIN + 1, i64::MAX - 1),
87        };
88        Extent {
89            start: edge(self.start),
90            end: edge(self.end),
91        }
92    }
93
94    pub fn secs(rate: u32, start_secs: f64, end_secs: f64) -> Extent {
95        let at = |secs: f64| (secs * f64::from(rate)).round() as i64;
96        Extent::new(at(start_secs), at(end_secs))
97    }
98
99    pub fn len(&self) -> usize {
100        debug_assert!(self.is_bounded(), "an unbounded extent has no length");
101        (self.end - self.start) as usize
102    }
103
104    pub fn is_empty(&self) -> bool {
105        self.end == self.start
106    }
107
108    pub fn start_secs(&self, rate: u32) -> f64 {
109        self.start as f64 / f64::from(rate)
110    }
111
112    pub fn span_secs(&self, rate: u32) -> f64 {
113        self.len() as f64 / f64::from(rate)
114    }
115}
116
117/// Sample `n` stands at `a*n/d` samples of `rate`, in lowest terms, `a, d > 0`: every grid
118/// starts at t = 0, so a sample index is the one clock every read and key counts in.
119#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord)]
120pub struct Grid {
121    pub rate: u32,
122    pub a: i128,
123    pub d: i128,
124}
125
126impl Grid {
127    pub const fn of(rate: u32) -> Grid {
128        Grid { rate, a: 1, d: 1 }
129    }
130
131    pub fn finer(rate: u32, scale: usize) -> Grid {
132        Grid {
133            rate,
134            a: 1,
135            d: scale as i128,
136        }
137    }
138
139    pub fn is_rate(&self) -> bool {
140        (self.a, self.d) == (1, 1)
141    }
142
143    pub(crate) fn exact(&self, n: i64) -> Option<(i128, i128)> {
144        let num = self.a.checked_mul(i128::from(n))?;
145        Some((num, self.d.checked_mul(i128::from(self.rate))?))
146    }
147
148    /// One quotient: correctly rounded while both integers are under 2^53; past that each
149    /// integer rounds once converting and the quotient once more.
150    pub fn instant(&self, n: i64) -> f64 {
151        if self.is_rate() {
152            return n as f64 / f64::from(self.rate);
153        }
154        let num = self.a.saturating_mul(i128::from(n));
155        num as f64 / self.d.saturating_mul(i128::from(self.rate)) as f64
156    }
157
158    /// The first sample whose exact instant is at or past `edge` read as the decimal it prints
159    /// as: the one rule every crop and indicator edge meets the grid by.
160    pub fn first_at(&self, edge: f64) -> i64 {
161        let beyond = if edge < 0.0 { i64::MIN } else { i64::MAX };
162        if !edge.is_finite() {
163            return beyond;
164        }
165        match self.decimal_steps(edge) {
166            Some(n) => i64::try_from(n).unwrap_or(beyond),
167            None => self.step_at(edge, Round::Ceil).unwrap_or(beyond),
168        }
169    }
170
171    /// An edge `first_at` meets at sample `n` exactly, beside `n`'s instant.
172    pub fn edge_at(&self, n: i64) -> Option<f64> {
173        let t = self.instant(n);
174        [t, t.next_down(), t.next_up()]
175            .into_iter()
176            .find(|edge| self.first_at(*edge) == n)
177    }
178
179    /// `ceil(edge * d * rate / a)`, `edge` as the shortest decimal that prints it.
180    fn decimal_steps(&self, edge: f64) -> Option<i128> {
181        self.steps(decimal(edge)?)
182    }
183
184    fn steps(&self, (num, den): (i128, i128)) -> Option<i128> {
185        let top = num.checked_mul(self.d.checked_mul(i128::from(self.rate))?)?;
186        let bottom = den.checked_mul(self.a)?;
187        Some(top.div_euclid(bottom) + i128::from(top.rem_euclid(bottom) != 0))
188    }
189
190    /// `first_at` of `edge - by`, the two decimals summed exactly rather than their doubles.
191    pub fn first_past(&self, edge: f64, by: f64) -> i64 {
192        let beyond = if edge - by < 0.0 { i64::MIN } else { i64::MAX };
193        let exact = || {
194            let ((a, p), (b, q)) = (decimal(edge)?, decimal(by)?);
195            let den = p.max(q);
196            let num = a
197                .checked_mul(den / p)?
198                .checked_sub(b.checked_mul(den / q)?)?;
199            self.steps((num, den))
200        };
201        match (edge.is_finite() && by.is_finite()).then(exact).flatten() {
202            Some(n) => i64::try_from(n).unwrap_or(beyond),
203            None => self.first_at(edge - by),
204        }
205    }
206
207    /// Whether sample `n` lies in `[l, r)` by `first_at`, which only a tie of its instant needs.
208    pub(crate) fn inside(&self, n: i64, l: f64, r: f64) -> bool {
209        let t = self.instant(n);
210        let reached = |edge: f64| match t == edge {
211            true => n >= self.first_at(edge),
212            false => t > edge,
213        };
214        reached(l) && !reached(r)
215    }
216
217    pub fn position(&self, n: i64) -> f64 {
218        self.a.saturating_mul(i128::from(n)) as f64 / self.d as f64
219    }
220
221    pub fn sr(&self) -> f64 {
222        f64::from(self.rate) * self.d as f64 / self.a as f64
223    }
224
225    pub fn count(&self, t: f64) -> f64 {
226        t * f64::from(self.rate) * self.d as f64 / self.a as f64
227    }
228
229    /// The step `t` falls nearest, rounded exactly from `t`'s own binary value; `None` where
230    /// that is past what the integers hold.
231    pub fn step_at(&self, t: f64, round: Round) -> Option<i64> {
232        match self.is_rate() {
233            true => rated(t, self.rate, round).or_else(|| self.wide(t, round)),
234            false => self.wide(t, round),
235        }
236    }
237
238    fn wide(&self, t: f64, round: Round) -> Option<i64> {
239        #[cfg(test)]
240        WIDE.with(|w| w.set(w.get() + 1));
241        if !t.is_finite() {
242            return None;
243        }
244        let bits = t.abs().to_bits();
245        let (exp, frac) = ((bits >> 52) as i32, (bits & ((1 << 52) - 1)) as i128);
246        let (mantissa, shift) = match exp {
247            0 => (frac, 1074),
248            e => (frac | (1 << 52), 1075 - e),
249        };
250        let zeros = mantissa.trailing_zeros().min(127) as i32;
251        let (mantissa, shift) = match mantissa {
252            0 => (0, 0),
253            m => (m >> zeros, shift - zeros),
254        };
255        let mantissa = if t < 0.0 { -mantissa } else { mantissa };
256        let (num, den) = match shift {
257            s if s <= 0 => (mantissa.checked_mul(1i128.checked_shl((-s) as u32)?)?, 1),
258            s if s < 127 => (mantissa, 1i128 << s),
259            _ => return None,
260        };
261        let top = num.checked_mul(self.d.checked_mul(i128::from(self.rate))?)?;
262        let bottom = self.a.checked_mul(den)?;
263        let (floor, rem) = (top.div_euclid(bottom), top.rem_euclid(bottom));
264        let k = match round {
265            Round::Floor => floor,
266            Round::Ceil => floor + i128::from(rem != 0),
267            Round::Even => match (2 * rem).cmp(&bottom) {
268                std::cmp::Ordering::Less => floor,
269                std::cmp::Ordering::Greater => floor + 1,
270                std::cmp::Ordering::Equal => floor + floor.rem_euclid(2),
271            },
272        };
273        i64::try_from(k).ok()
274    }
275}
276
277/// `t*rate` rounded from the double nearest it and its exact error: halves of `t` times `rate`
278/// exactly, summed exactly; the error decides only a sum on a whole or half step. `None` where
279/// `wide` may refuse or a product could round.
280fn rated(t: f64, rate: u32, round: Round) -> Option<i64> {
281    let held = t == 0.0 || (2f64.powi(-74)..2f64.powi(51)).contains(&t.abs());
282    if !held || rate >= 1 << 26 {
283        return None;
284    }
285    let r = f64::from(rate);
286    let split = 134_217_729.0 * t;
287    let hi = split - (split - t);
288    let (x, y) = (hi * r, (t - hi) * r);
289    let p = x + y;
290    let back = p - x;
291    let e = (x - (p - back)) + (y - back);
292    if p.abs() >= 2f64.powi(51) {
293        return None;
294    }
295    let f = p.floor();
296    let k = match round {
297        Round::Floor if p == f && e < 0.0 => f - 1.0,
298        Round::Floor => f,
299        Round::Ceil if p == p.ceil() && e > 0.0 => p + 1.0,
300        Round::Ceil => p.ceil(),
301        Round::Even => match (p - f).total_cmp(&0.5) {
302            std::cmp::Ordering::Less => f,
303            std::cmp::Ordering::Greater => f + 1.0,
304            std::cmp::Ordering::Equal if e > 0.0 => f + 1.0,
305            std::cmp::Ordering::Equal if e < 0.0 => f,
306            std::cmp::Ordering::Equal => f + f.rem_euclid(2.0),
307        },
308    };
309    Some(k as i64)
310}
311
312#[cfg(test)]
313thread_local! {
314    pub(crate) static WIDE: std::cell::Cell<u64> = const { std::cell::Cell::new(0) };
315}
316
317#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord)]
318pub enum Round {
319    Even,
320    Floor,
321    Ceil,
322}
323
324#[cfg(test)]
325mod tests {
326    use super::{Grid, Round, rated};
327
328    /// 0.8333333333333334 s lies past 5/6 s, sample 40000's at 48 kHz, though the doubles tie.
329    #[test]
330    fn an_edge_in_seconds_tying_an_instant_is_decided_by_its_decimal() {
331        let cd = Grid::of(44_100);
332        assert_eq!(cd.first_at(0.1), 4410);
333        assert!(cd.inside(4410, 0.1, 1.0));
334        let grid = Grid::of(48_000);
335        let edge = 60.0 / 72.0;
336        assert_eq!(grid.instant(40_000), edge);
337        assert_eq!(grid.first_at(edge), 40_001);
338        assert!(!grid.inside(40_000, edge, 2.0));
339        assert!(grid.inside(40_000, 0.0, edge));
340    }
341
342    fn instants(rate: u32) -> Vec<f64> {
343        let r = f64::from(rate);
344        let mut out = vec![0.0, -0.0, 1e-300, -1e-300, 1e300, f64::MIN_POSITIVE];
345        let mut seed = 0x9e37_79b9_7f4a_7c15u64;
346        for k in -2_000i64..2_000 {
347            for half in [0.0, 0.5] {
348                let at = (k as f64 + half) / r;
349                let mut near = at;
350                for _ in 0..3 {
351                    near = near.next_up();
352                    out.push(near);
353                }
354                near = at;
355                for _ in 0..3 {
356                    near = near.next_down();
357                    out.push(near);
358                }
359                out.push(at);
360                out.push(at * 1e9);
361            }
362            seed = seed.wrapping_mul(6_364_136_223_846_793_005).wrapping_add(1);
363            out.push(f64::from_bits(seed >> 2) * if seed & 1 == 0 { 1.0 } else { -1.0 });
364        }
365        out
366    }
367
368    #[test]
369    fn a_step_rounds_alike_in_doubles_and_in_wide_integers() {
370        for rate in [1, 8_000, 44_100, 48_000, 88_200, 96_000, 192_000] {
371            let grid = Grid::of(rate);
372            for t in instants(rate) {
373                for round in [Round::Floor, Round::Ceil, Round::Even] {
374                    if let Some(k) = rated(t, rate, round) {
375                        assert_eq!(Some(k), grid.wide(t, round), "{t:e} at {rate} {round:?}");
376                    }
377                }
378            }
379        }
380    }
381}
sva_samples/grid.rs

sva_samples/
grid.rs
