idsp/
num.rs

1use num_traits::{AsPrimitive, Float, Num};
2
3/// Helper trait unifying fixed point and floating point coefficients/samples
4pub trait Coefficient: 'static + Copy + Num + AsPrimitive<Self::ACCU> {
5    /// Multiplicative identity
6    const ONE: Self;
7    /// Negative multiplicative identity, equal to `-Self::ONE`.
8    const NEG_ONE: Self;
9    /// Additive identity
10    const ZERO: Self;
11    /// Lowest value
12    const MIN: Self;
13    /// Highest value
14    const MAX: Self;
15    /// Accumulator type
16    type ACCU: AsPrimitive<Self> + Num;
17
18    /// Proper scaling and potentially using a wide accumulator.
19    /// Clamp `self` such that `min <= self <= max`.
20    /// Undefined result if `max < min`.
21    fn macc(self, s: Self::ACCU, min: Self, max: Self, e1: Self) -> (Self, Self);
22
23    /// Clamp to between min and max
24    ///
25    /// Undefined if `min > max`.
26    fn clip(self, min: Self, max: Self) -> Self;
27
28    /// Multiplication (scaled)
29    fn mul_scaled(self, other: Self) -> Self;
30
31    /// Division (scaled)
32    fn div_scaled(self, other: Self) -> Self;
33
34    /// Scale and quantize a floating point value.
35    fn quantize<C>(value: C) -> Self
36    where
37        Self: AsPrimitive<C>,
38        C: Float + AsPrimitive<Self>;
39    // TODO: range check and Result
40}
41
42macro_rules! impl_float {
43    ($T:ty) => {
44        impl Coefficient for $T {
45            const ONE: Self = 1.0;
46            const NEG_ONE: Self = -1.0;
47            const ZERO: Self = 0.0;
48            const MIN: Self = <$T>::NEG_INFINITY;
49            const MAX: Self = <$T>::INFINITY;
50            type ACCU = Self;
51
52            #[inline]
53            fn macc(self, s: Self::ACCU, min: Self, max: Self, _e1: Self) -> (Self, Self) {
54                ((self + s).clip(min, max), 0.0)
55            }
56
57            #[inline]
58            fn clip(self, min: Self, max: Self) -> Self {
59                // <$T>::clamp() is slow and checks
60                self.max(min).min(max)
61            }
62
63            #[inline]
64            fn div_scaled(self, other: Self) -> Self {
65                self / other
66            }
67
68            #[inline]
69            fn mul_scaled(self, other: Self) -> Self {
70                self * other
71            }
72
73            #[inline]
74            fn quantize<C: Float + AsPrimitive<Self>>(value: C) -> Self {
75                value.as_()
76            }
77        }
78    };
79}
80impl_float!(f32);
81impl_float!(f64);
82
83macro_rules! impl_int {
84    ($T:ty, $U:ty, $A:ty, $Q:literal) => {
85        impl Coefficient for $T {
86            const ONE: Self = 1 << $Q;
87            const NEG_ONE: Self = -1 << $Q;
88            const ZERO: Self = 0;
89            const MIN: Self = <$T>::MIN;
90            const MAX: Self = <$T>::MAX;
91            type ACCU = $A;
92
93            #[inline]
94            fn macc(self, mut s: Self::ACCU, min: Self, max: Self, e1: Self) -> (Self, Self) {
95                const S: usize = core::mem::size_of::<$T>() * 8;
96                // Guard bits
97                const G: usize = S - $Q;
98                // Combine offset (u << $Q) with previous quantization error e1
99                s += (((self >> G) as $A) << S) | (((self << $Q) | e1) as $U as $A);
100                // Ord::clamp() is slow and checks
101                // This clamping truncates the lowest G bits of the value and the limits.
102                debug_assert_eq!(min & ((1 << G) - 1), 0);
103                debug_assert_eq!(max & ((1 << G) - 1), (1 << G) - 1);
104                let y0 = if (s >> S) as $T < (min >> G) {
105                    min
106                } else if (s >> S) as $T > (max >> G) {
107                    max
108                } else {
109                    (s >> $Q) as $T
110                };
111                // Quantization error
112                let e0 = s as $T & ((1 << $Q) - 1);
113                (y0, e0)
114            }
115
116            #[inline]
117            fn clip(self, min: Self, max: Self) -> Self {
118                // Ord::clamp() is slow and checks
119                if self < min {
120                    min
121                } else if self > max {
122                    max
123                } else {
124                    self
125                }
126            }
127
128            #[inline]
129            fn div_scaled(self, other: Self) -> Self {
130                (((self as $A) << $Q) / other as $A) as $T
131            }
132
133            #[inline]
134            fn mul_scaled(self, other: Self) -> Self {
135                (((1 << ($Q - 1)) + self as $A * other as $A) >> $Q) as $T
136            }
137
138            #[inline]
139            fn quantize<C>(value: C) -> Self
140            where
141                Self: AsPrimitive<C>,
142                C: Float + AsPrimitive<Self>,
143            {
144                (value * (1 << $Q).as_()).round().as_()
145            }
146        }
147    };
148}
149// Q2.X chosen to be able to exactly and inclusively represent -2 as `-1 << X + 1`
150// This is necessary to meet a1 = -2
151// It also create 2 guard bits for clamping in the accumulator which is often enough.
152impl_int!(i8, u8, i16, 6);
153impl_int!(i16, u16, i32, 14);
154impl_int!(i32, u32, i64, 30);
155impl_int!(i64, u64, i128, 62);