idsp/
num.rs

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
use num_traits::{AsPrimitive, Float, Num};

/// Helper trait unifying fixed point and floating point coefficients/samples
pub trait Coefficient: 'static + Copy + Num + AsPrimitive<Self::ACCU> {
    /// Multiplicative identity
    const ONE: Self;
    /// Negative multiplicative identity, equal to `-Self::ONE`.
    const NEG_ONE: Self;
    /// Additive identity
    const ZERO: Self;
    /// Lowest value
    const MIN: Self;
    /// Highest value
    const MAX: Self;
    /// Accumulator type
    type ACCU: AsPrimitive<Self> + Num;

    /// Proper scaling and potentially using a wide accumulator.
    /// Clamp `self` such that `min <= self <= max`.
    /// Undefined result if `max < min`.
    fn macc(self, s: Self::ACCU, min: Self, max: Self, e1: Self) -> (Self, Self);

    /// Clamp to between min and max
    ///
    /// Undefined if `min > max`.
    fn clip(self, min: Self, max: Self) -> Self;

    /// Multiplication (scaled)
    fn mul_scaled(self, other: Self) -> Self;

    /// Division (scaled)
    fn div_scaled(self, other: Self) -> Self;

    /// Scale and quantize a floating point value.
    fn quantize<C>(value: C) -> Self
    where
        Self: AsPrimitive<C>,
        C: Float + AsPrimitive<Self>;
    // TODO: range check and Result
}

macro_rules! impl_float {
    ($T:ty) => {
        impl Coefficient for $T {
            const ONE: Self = 1.0;
            const NEG_ONE: Self = -1.0;
            const ZERO: Self = 0.0;
            const MIN: Self = <$T>::NEG_INFINITY;
            const MAX: Self = <$T>::INFINITY;
            type ACCU = Self;

            #[inline]
            fn macc(self, s: Self::ACCU, min: Self, max: Self, _e1: Self) -> (Self, Self) {
                ((self + s).clip(min, max), 0.0)
            }

            #[inline]
            fn clip(self, min: Self, max: Self) -> Self {
                // <$T>::clamp() is slow and checks
                self.max(min).min(max)
            }

            #[inline]
            fn div_scaled(self, other: Self) -> Self {
                self / other
            }

            #[inline]
            fn mul_scaled(self, other: Self) -> Self {
                self * other
            }

            #[inline]
            fn quantize<C: Float + AsPrimitive<Self>>(value: C) -> Self {
                value.as_()
            }
        }
    };
}
impl_float!(f32);
impl_float!(f64);

macro_rules! impl_int {
    ($T:ty, $U:ty, $A:ty, $Q:literal) => {
        impl Coefficient for $T {
            const ONE: Self = 1 << $Q;
            const NEG_ONE: Self = -1 << $Q;
            const ZERO: Self = 0;
            const MIN: Self = <$T>::MIN;
            const MAX: Self = <$T>::MAX;
            type ACCU = $A;

            #[inline]
            fn macc(self, mut s: Self::ACCU, min: Self, max: Self, e1: Self) -> (Self, Self) {
                const S: usize = core::mem::size_of::<$T>() * 8;
                // Guard bits
                const G: usize = S - $Q;
                // Combine offset (u << $Q) with previous quantization error e1
                s += (((self >> G) as $A) << S) | (((self << $Q) | e1) as $U as $A);
                // Ord::clamp() is slow and checks
                // This clamping truncates the lowest G bits of the value and the limits.
                debug_assert_eq!(min & ((1 << G) - 1), 0);
                debug_assert_eq!(max & ((1 << G) - 1), (1 << G) - 1);
                let y0 = if (s >> S) as $T < (min >> G) {
                    min
                } else if (s >> S) as $T > (max >> G) {
                    max
                } else {
                    (s >> $Q) as $T
                };
                // Quantization error
                let e0 = s as $T & ((1 << $Q) - 1);
                (y0, e0)
            }

            #[inline]
            fn clip(self, min: Self, max: Self) -> Self {
                // Ord::clamp() is slow and checks
                if self < min {
                    min
                } else if self > max {
                    max
                } else {
                    self
                }
            }

            #[inline]
            fn div_scaled(self, other: Self) -> Self {
                (((self as $A) << $Q) / other as $A) as $T
            }

            #[inline]
            fn mul_scaled(self, other: Self) -> Self {
                (((1 << ($Q - 1)) + self as $A * other as $A) >> $Q) as $T
            }

            #[inline]
            fn quantize<C>(value: C) -> Self
            where
                Self: AsPrimitive<C>,
                C: Float + AsPrimitive<Self>,
            {
                (value * (1 << $Q).as_()).round().as_()
            }
        }
    };
}
// Q2.X chosen to be able to exactly and inclusively represent -2 as `-1 << X + 1`
// This is necessary to meet a1 = -2
// It also create 2 guard bits for clamping in the accumulator which is often enough.
impl_int!(i8, u8, i16, 6);
impl_int!(i16, u16, i32, 14);
impl_int!(i32, u32, i64, 30);
impl_int!(i64, u64, i128, 62);