diff --git a/celt/mathops.h b/celt/mathops.h index 0e6112e8..060a8c79 100644 --- 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 @@ static OPUS_INLINE float celt_log2(float x) 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 @@ static OPUS_INLINE float celt_exp2(float x) 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; } diff --git a/celt/tests/test_unit_mathops.c b/celt/tests/test_unit_mathops.c index 7e8f69c4..8849d8ed 100644 --- 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 testlog2(void) 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) {