ref: e75503be9419dcc6e4dc3c998948dc2374ae6d07
parent: 4ba06d9e324a76d1ae28c9b80cae80f5d2aa66bb
author: Yunho Huh <yunho@google.com>
date: Sat Nov 2 08:27:09 EDT 2024
Improve log2's precision to 20-24 bits. The log2 function was approximated using lolremez, achieving an accuracy of less than 1.4*10^-8 within the range of 1 to 2. Signed-off-by: Jean-Marc Valin <jeanmarcv@google.com>
--- a/celt/mathops.h
+++ b/celt/mathops.h
@@ -143,15 +143,45 @@
#define frac_div32_q29(a,b) frac_div32(a,b)
#ifdef FLOAT_APPROX
+/* Calculates the base-2 logarithm (log2(x)) of a number. It is designed for
+ * systems using radix-2 floating-point representation, with the exponent
+ * located at bits 23 to 30 and an offset of 127. Note that special cases like
+ * denormalized numbers, positive/negative infinity, and NaN are not handled.
+ * log2(x) = log2(x^exponent * mantissa)
+ * = exponent + log2(mantissa) */
-/* Note: This assumes radix-2 floating point with the exponent at bits 23..30 and an offset of 127
- denorm, +/- inf and NaN are *not* handled */
+/* Log2 x normalization single precision coefficients calculated by
+ * 1 / (1 + 0.125 * index).
+ * Coefficients in Double Precision
+ * double log2_x_norm_coeff[8] = {
+ * 1.0000000000000000000, 8.888888888888888e-01,
+ * 8.000000000000000e-01, 7.272727272727273e-01,
+ * 6.666666666666666e-01, 6.153846153846154e-01,
+ * 5.714285714285714e-01, 5.333333333333333e-01} */
+static const float log2_x_norm_coeff[8] = {
+ 1.000000000000000000000000000f, 8.88888895511627197265625e-01f,
+ 8.00000000000000000000000e-01f, 7.27272748947143554687500e-01f,
+ 6.66666686534881591796875e-01f, 6.15384638309478759765625e-01f,
+ 5.71428596973419189453125e-01f, 5.33333361148834228515625e-01f};
-/** Base-2 log approximation (log2(x)). */
+/* Log2 y normalization single precision coefficients calculated by
+ * log2(1 + 0.125 * index).
+ * Coefficients in Double Precision
+ * double log2_y_norm_coeff[8] = {
+ * 0.0000000000000000000, 1.699250014423124e-01,
+ * 3.219280948873623e-01, 4.594316186372973e-01,
+ * 5.849625007211562e-01, 7.004397181410922e-01,
+ * 8.073549220576041e-01, 9.068905956085185e-01}; */
+static const float log2_y_norm_coeff[8] = {
+ 0.0000000000000000000000000000f, 1.699250042438507080078125e-01f,
+ 3.219280838966369628906250e-01f, 4.594316184520721435546875e-01f,
+ 5.849624872207641601562500e-01f, 7.004396915435791015625000e-01f,
+ 8.073549270629882812500000e-01f, 9.068905711174011230468750e-01f};
+
static OPUS_INLINE float celt_log2(float x)
{
int integer;
- float frac;
+ opus_int32 range_idx;
union {
float f;
opus_uint32 i;
@@ -159,10 +189,29 @@
in.f = x;
integer = (in.i>>23)-127;
in.i -= (opus_uint32)integer<<23;
- frac = in.f - 1.5f;
- frac = -0.41445418f + frac*(0.95909232f
- + frac*(-0.33951290f + frac*0.16541097f));
- return 1+integer+frac;
+
+ /* Normalize the mantissa range from [1, 2] to [1,1.125], and then shift x
+ * by 1.0625 to [-0.0625, 0.0625]. */
+ range_idx = (in.i >> 20) & 0x7;
+ in.f = in.f * log2_x_norm_coeff[range_idx] - 1.0625f;
+
+ /* Polynomial coefficients approximated in the [1, 1.125] range.
+ * Lolremez command: lolremez --degree 4 --range -0.0625:0.0625
+ * "log(x+1.0625)/log(2)"
+ * Coefficients in Double Precision
+ * A0: 8.7462840624502679e-2 A1: 1.3578296070972002
+ * A2: -6.3897703690210047e-1 A3: 4.0197125617419959e-1
+ * A4: -2.8415445877832832e-1 */
+ #define LOG2_COEFF_A0 8.74628424644470214843750000e-02f
+ #define LOG2_COEFF_A1 1.357829570770263671875000000000f
+ #define LOG2_COEFF_A2 -6.3897705078125000000000000e-01f
+ #define LOG2_COEFF_A3 4.01971250772476196289062500e-01f
+ #define LOG2_COEFF_A4 -2.8415444493293762207031250e-01f
+ in.f = LOG2_COEFF_A0 + in.f * (LOG2_COEFF_A1
+ + in.f * (LOG2_COEFF_A2
+ + in.f * (LOG2_COEFF_A3
+ + in.f * (LOG2_COEFF_A4))));
+ return integer + in.f + log2_y_norm_coeff[range_idx];
}
/** Base-2 exponential approximation (2^x). */
--- a/celt/tests/test_unit_mathops.c
+++ b/celt/tests/test_unit_mathops.c
@@ -143,15 +143,26 @@
void testlog2(void)
{
float x;
+ float error_threshold = 2.e-06;
+ float max_error = 0;
for (x=0.001f;x<1677700.0;x+=(x/8.0))
{
float error = fabs((1.442695040888963387*log(x))-celt_log2(x));
- if (error>0.0009)
+ if (max_error < error)
{
- fprintf (stderr, "celt_log2 failed: fabs((1.442695040888963387*log(x))-celt_log2(x))>0.001 (x = %f, error = %f)\n", x,error);
+ max_error = error;
+ }
+
+ if (error > error_threshold)
+ {
+ fprintf (stderr,
+ "celt_log2 failed: "
+ "fabs((1.442695040888963387*log(x))-celt_log2(x))>%15.25e "
+ "(x = %f, error = %15.25e)\n", error_threshold, x, error);
ret = 1;
}
}
+ fprintf (stdout, "celt_log2 max_error: %15.25e\n", max_error);
}
void testexp2(void)
@@ -171,15 +182,26 @@
void testexp2log2(void)
{
float x;
+ float error_threshold = 5.0e-04;
+ float max_error = 0;
for (x=-11.0;x<24.0;x+=0.0007f)
{
float error = fabs(x-(celt_log2(celt_exp2(x))));
- if (error>0.001)
+ if (max_error < error)
{
- fprintf (stderr, "celt_log2/celt_exp2 failed: fabs(x-(celt_log2(celt_exp2(x))))>0.001 (x = %f, error = %f)\n", x,error);
+ max_error = error;
+ }
+
+ if (error > error_threshold)
+ {
+ fprintf (stderr,
+ "celt_log2/celt_exp2 failed: "
+ "fabs(x-(celt_log2(celt_exp2(x))))>%15.25e "
+ "(x = %f, error = %15.25e)\n", error_threshold, x, error);
ret = 1;
}
}
+ fprintf (stdout, "celt_exp2, celt_log2 max_error: %15.25e\n", max_error);
}
#else
void testlog2(void)
--
⑨