1 /******************************************************************************
2 *
3 * Copyright 2022 Google LLC
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at:
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 ******************************************************************************/
18
19 /**
20 * LC3 - Mathematics function approximation
21 */
22
23 #ifndef __LC3_FASTMATH_H
24 #define __LC3_FASTMATH_H
25
26 #include <stdint.h>
27 #include <math.h>
28
29
30 /**
31 * Fast 2^n approximation
32 * x Operand, range -8 to 8
33 * return 2^x approximation (max relative error ~ 7e-6)
34 */
fast_exp2f(float x)35 static inline float fast_exp2f(float x)
36 {
37 float y;
38
39 /* --- Polynomial approx in range -0.5 to 0.5 --- */
40
41 static const float c[] = { 1.27191277e-09, 1.47415221e-07,
42 1.35510312e-05, 9.38375815e-04, 4.33216946e-02 };
43
44 y = ( c[0]) * x;
45 y = (y + c[1]) * x;
46 y = (y + c[2]) * x;
47 y = (y + c[3]) * x;
48 y = (y + c[4]) * x;
49 y = (y + 1.f);
50
51 /* --- Raise to the power of 16 --- */
52
53 y = y*y;
54 y = y*y;
55 y = y*y;
56 y = y*y;
57
58 return y;
59 }
60
61 /**
62 * Fast log2(x) approximation
63 * x Operand, greater than 0
64 * return log2(x) approximation (max absolute error ~ 1e-4)
65 */
fast_log2f(float x)66 static inline float fast_log2f(float x)
67 {
68 float y;
69 int e;
70
71 /* --- Polynomial approx in range 0.5 to 1 --- */
72
73 static const float c[] = {
74 -1.29479677, 5.11769018, -8.42295281, 8.10557963, -3.50567360 };
75
76 x = frexpf(x, &e);
77
78 y = ( c[0]) * x;
79 y = (y + c[1]) * x;
80 y = (y + c[2]) * x;
81 y = (y + c[3]) * x;
82 y = (y + c[4]);
83
84 /* --- Add log2f(2^e) and return --- */
85
86 return e + y;
87 }
88
89 /**
90 * Fast log10(x) approximation
91 * x Operand, greater than 0
92 * return log10(x) approximation (max absolute error ~ 1e-4)
93 */
fast_log10f(float x)94 static inline float fast_log10f(float x)
95 {
96 return log10f(2) * fast_log2f(x);
97 }
98
99 /**
100 * Fast `10 * log10(x)` (or dB) approximation in fixed Q16
101 * x Operand, in range 2^-63 to 2^63 (1e-19 to 1e19)
102 * return 10 * log10(x) in fixed Q16 (-190 to 192 dB)
103 *
104 * - The 0 value is accepted and return the minimum value ~ -191dB
105 * - This function assumed that float 32 bits is coded IEEE 754
106 */
fast_db_q16(float x)107 static inline int32_t fast_db_q16(float x)
108 {
109 /* --- Table in Q15 --- */
110
111 static const uint16_t t[][2] = {
112
113 /* [n][0] = 10 * log10(2) * log2(1 + n/32), with n = [0..15] */
114 /* [n][1] = [n+1][0] - [n][0] (while defining [16][0]) */
115
116 { 0, 4379 }, { 4379, 4248 }, { 8627, 4125 }, { 12753, 4009 },
117 { 16762, 3899 }, { 20661, 3795 }, { 24456, 3697 }, { 28153, 3603 },
118 { 31755, 3514 }, { 35269, 3429 }, { 38699, 3349 }, { 42047, 3272 },
119 { 45319, 3198 }, { 48517, 3128 }, { 51645, 3061 }, { 54705, 2996 },
120
121 /* [n][0] = 10 * log10(2) * log2(1 + n/32) - 10 * log10(2) / 2, */
122 /* with n = [16..31] */
123 /* [n][1] = [n+1][0] - [n][0] (while defining [32][0]) */
124
125 { 8381, 2934 }, { 11315, 2875 }, { 14190, 2818 }, { 17008, 2763 },
126 { 19772, 2711 }, { 22482, 2660 }, { 25142, 2611 }, { 27754, 2564 },
127 { 30318, 2519 }, { 32837, 2475 }, { 35312, 2433 }, { 37744, 2392 },
128 { 40136, 2352 }, { 42489, 2314 }, { 44803, 2277 }, { 47080, 2241 },
129
130 };
131
132 /* --- Approximation ---
133 *
134 * 10 * log10(x^2) = 10 * log10(2) * log2(x^2)
135 *
136 * And log2(x^2) = 2 * log2( (1 + m) * 2^e )
137 * = 2 * (e + log2(1 + m)) , with m in range [0..1]
138 *
139 * Split the float values in :
140 * e2 Double value of the exponent (2 * e + k)
141 * hi High 5 bits of mantissa, for precalculated result `t[hi][0]`
142 * lo Low 16 bits of mantissa, for linear interpolation `t[hi][1]`
143 *
144 * Two cases, from the range of the mantissa :
145 * 0 to 0.5 `k = 0`, use 1st part of the table
146 * 0.5 to 1 `k = 1`, use 2nd part of the table */
147
148 union { float f; uint32_t u; } x2 = { .f = x*x };
149
150 int e2 = (int)(x2.u >> 22) - 2*127;
151 int hi = (x2.u >> 18) & 0x1f;
152 int lo = (x2.u >> 2) & 0xffff;
153
154 return e2 * 49321 + t[hi][0] + ((t[hi][1] * lo) >> 16);
155 }
156
157
158 #endif /* __LC3_FASTMATH_H */
159
