ref: 255f0130353213d2b58af7b521fb88747ff6c74b
parent: 7fa23b4ed0101fa4cf1afc8657fab581d8fd6fe0
author: Yunho Huh <yunho@google.com>
date: Sat Nov 2 09:02:16 EDT 2024
Improve exp2's precision to 20-24 bits. The exp2 function was approximated using lolremez, achieving an accuracy of less than 2*10^-7 within the range of 0 to 1. Signed-off-by: Jean-Marc Valin <jeanmarcv@google.com>
--- a/celt/mathops.h
+++ b/celt/mathops.h
@@ -1,7 +1,7 @@
/* Copyright (c) 2002-2008 Jean-Marc Valin
Copyright (c) 2007-2008 CSIRO
Copyright (c) 2007-2009 Xiph.Org Foundation
- Written by Jean-Marc Valin */
+ Written by Jean-Marc Valin, and Yunho Huh */
/**
@file mathops.h
@brief Various math functions
@@ -214,7 +214,11 @@
return integer + in.f + log2_y_norm_coeff[range_idx];
}
-/** Base-2 exponential approximation (2^x). */
+/* Calculates an approximation of 2^x. The approximation was achieved by
+ * employing a base-2 exponential function and utilizing a Remez approximation
+ * of order 5, ensuring a controlled relative error.
+ * exp2(x) = exp2(integer + fraction)
+ * ~ exp2(integer) + exp2(fraction) */
static OPUS_INLINE float celt_exp2(float x)
{
int integer;
@@ -227,9 +231,22 @@
if (integer < -50)
return 0;
frac = x-integer;
- /* K0 = 1, K1 = log(2), K2 = 3-4*log(2), K3 = 3*log(2) - 2 */
- res.f = 0.99992522f + frac * (0.69583354f
- + frac * (0.22606716f + 0.078024523f*frac));
+
+ /* Polynomial coefficients approximated in the [0, 1] range.
+ * Lolremez command: lolremez --degree 5 --range 0:1
+ * "exp(x*0.693147180559945)" "exp(x*0.693147180559945)"
+ * NOTE: log(2) ~ 0.693147180559945 */
+ #define EXP2_COEFF_A0 9.999999403953552246093750000000e-01f
+ #define EXP2_COEFF_A1 6.931530833244323730468750000000e-01f
+ #define EXP2_COEFF_A2 2.401536107063293457031250000000e-01f
+ #define EXP2_COEFF_A3 5.582631751894950866699218750000e-02f
+ #define EXP2_COEFF_A4 8.989339694380760192871093750000e-03f
+ #define EXP2_COEFF_A5 1.877576694823801517486572265625e-03f
+ res.f = EXP2_COEFF_A0 + frac * (EXP2_COEFF_A1
+ + frac * (EXP2_COEFF_A2
+ + frac * (EXP2_COEFF_A3
+ + frac * (EXP2_COEFF_A4
+ + frac * (EXP2_COEFF_A5)))));
res.i = (res.i + ((opus_uint32)integer<<23)) & 0x7fffffff;
return res.f;
}
--- a/celt/tests/test_unit_mathops.c
+++ b/celt/tests/test_unit_mathops.c
@@ -1,6 +1,7 @@
/* Copyright (c) 2008-2011 Xiph.Org Foundation, Mozilla Corporation,
Gregory Maxwell
- Written by Jean-Marc Valin, Gregory Maxwell, and Timothy B. Terriberry */
+ Written by Jean-Marc Valin, Gregory Maxwell, Timothy B. Terriberry,
+ and Yunho Huh */
/*
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
@@ -170,21 +171,32 @@
void testexp2(void)
{
float x;
+ float error_threshold = 2.3e-07;
+ float max_error = 0;
for (x=-11.0;x<24.0;x+=0.0007f)
{
float error = fabs(x-(1.442695040888963387*log(celt_exp2(x))));
- if (error>0.0002)
+ if (max_error < error)
{
- fprintf (stderr, "celt_exp2 failed: fabs(x-(1.442695040888963387*log(celt_exp2(x))))>0.0005 (x = %f, error = %f)\n", x,error);
+ max_error = error;
+ }
+
+ if (error > error_threshold)
+ {
+ fprintf (stderr,
+ "celt_exp2 failed: "
+ "fabs(x-(1.442695040888963387*log(celt_exp2(x))))>%15.25e "
+ "(x = %f, error = %15.25e)\n", error_threshold, x, error);
ret = 1;
}
}
+ fprintf (stdout, "celt_exp2 max_error: %15.25e\n", max_error);
}
void testexp2log2(void)
{
float x;
- float error_threshold = 5.0e-04;
+ float error_threshold = 2.0e-06;
float max_error = 0;
for (x=-11.0;x<24.0;x+=0.0007f)
{
--
⑨