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 */ 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 */ 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 */ 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 */ 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