Enhancements the fixed-point approximations of non-linear functions.

Accuracy for rsqrt, rcp, cos, and log2 is now at the level of truncation error
 for the current output resolution of these functions.
sqrt and exp2 still have non-trivial algebraic error, but this cannot be
 reduced much further using the current method without additional computation.
Also updates the fast float approximations for log2 and exp2 with coefficients
 that give slightly lower maximum relative error.

Patch modified by Jean-Marc Valin to leave the cos approximation as is and
leave the check for x<-15 in exp2 as is.

1 files changed, 48 insertions, 6 deletions

diff --git a/tests/mathops-test.c b/tests/mathops-test.c
index 423d1e9..5f086ae 100644
--- a/tests/mathops-test.c
+++ b/tests/mathops-test.c
@@ -30,7 +30,7 @@ void testdiv(void)
 #else
       prod = val*i;
 #endif
-      if (fabs(prod-1) > .001)
+      if (fabs(prod-1) > .00025)
       {
          fprintf (stderr, "div failed: 1/%d="WORD" (product = %f)\n", i, val, prod);
          ret = 1;
@@ -47,7 +47,7 @@ void testsqrt(void)
       celt_word16 val;
       val = celt_sqrt(i);
       ratio = val/sqrt(i);
-      if (fabs(ratio - 1) > .001 && fabs(val-sqrt(i)) > 2)
+      if (fabs(ratio - 1) > .0005 && fabs(val-sqrt(i)) > 2)
       {
          fprintf (stderr, "sqrt failed: sqrt(%d)="WORD" (ratio = %f)\n", i, val, ratio);
          ret = 1;
@@ -81,7 +81,7 @@ void testlog2(void)
    for (x=0.001;x<1677700.0;x+=(x/8.0))
    {
       float error = fabs((1.442695040888963387*log(x))-celt_log2(x));
-      if (error>0.001)
+      if (error>0.0009)
       {
          fprintf (stderr, "celt_log2 failed: fabs((1.442695040888963387*log(x))-celt_log2(x))>0.001 (x = %f, error = %f)\n", x,error);
          ret = 1;    
@@ -95,7 +95,7 @@ void testexp2(void)
    for (x=-11.0;x<24.0;x+=0.0007)
    {
       float error = fabs(x-(1.442695040888963387*log(celt_exp2(x))));
-      if (error>0.0005)
+      if (error>0.0002)
       {
          fprintf (stderr, "celt_exp2 failed: fabs(x-(1.442695040888963387*log(celt_exp2(x))))>0.0005 (x = %f, error = %f)\n", x,error);
          ret = 1;    
@@ -117,6 +117,49 @@ void testexp2log2(void)
    }
 }
 #else
+void testlog2(void)
+{
+   celt_word32 x;
+   for (x=8;x<1073741824;x+=(x>>3))
+   {
+      float error = fabs((1.442695040888963387*log(x/16384.0))-celt_log2(x)/256.0);
+      if (error>0.003)
+      {
+         fprintf (stderr, "celt_log2 failed: x = %ld, error = %f\n", (long)x,error);
+         ret = 1;
+      }
+   }
+}
+
+void testexp2(void)
+{
+   celt_word16 x;
+   for (x=-32768;x<-30720;x++)
+   {
+      float error1 = fabs(x/2048.0-(1.442695040888963387*log(celt_exp2(x)/65536.0)));
+      float error2 = fabs(exp(0.6931471805599453094*x/2048.0)-celt_exp2(x)/65536.0);
+      if (error1>0.0002&&error2>0.00004)
+      {
+         fprintf (stderr, "celt_exp2 failed: x = "WORD", error1 = %f, error2 = %f\n", x,error1,error2);
+         ret = 1;
+      }
+   }
+}
+
+void testexp2log2(void)
+{
+   celt_word32 x;
+   for (x=8;x<65536;x+=(x>>3))
+   {
+      float error = fabs(x-0.25*celt_exp2(celt_log2(x)<<3))/16384;
+      if (error>0.004)
+      {
+         fprintf (stderr, "celt_log2/celt_exp2 failed: fabs(x-(celt_log2(celt_exp2(x))))>0.001 (x = %ld, error = %f)\n", (long)x,error);
+         ret = 1;
+      }
+   }
+}
+
 void testilog2(void)
 {
    celt_word32 x;
@@ -137,11 +180,10 @@ int main(void)
    testdiv();
    testsqrt();
    testrsqrt();
-#ifndef FIXED_POINT
    testlog2();
    testexp2();
    testexp2log2();
-#else
+#ifdef FIXED_POINT
    testilog2();
 #endif
    return ret;