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>

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
+   /* Polynomial coefficients approximated in the [0, 1] range.
-           + frac * (0.22606716f + 0.078024523f*frac));
+    * 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;
 }

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)
    {