Move skip coding into interp_bits2pulses().
[opus.git] / libcelt / mathops.h
index 2ec3a12..160c220 100644 (file)
 #include "entcode.h"
 #include "os_support.h"
 
-#ifndef OVERRIDE_CELT_ILOG2
-/** Integer log in base2. Undefined for zero and negative numbers */
-static inline celt_int16 celt_ilog2(celt_word32 x)
-{
-   celt_assert2(x>0, "celt_ilog2() only defined for strictly positive numbers");
-   return EC_ILOG(x)-1;
-}
-#endif
-
-#ifndef OVERRIDE_FIND_MAX16
-static inline int find_max16(celt_word16 *x, int len)
-{
-   celt_word16 max_corr=-VERY_LARGE16;
-   int i, id = 0;
-   for (i=0;i<len;i++)
-   {
-      if (x[i] > max_corr)
-      {
-         id = i;
-         max_corr = x[i];
-      }
-   }
-   return id;
-}
-#endif
-
-#ifndef OVERRIDE_FIND_MAX32
-static inline int find_max32(celt_word32 *x, int len)
-{
-   celt_word32 max_corr=-VERY_LARGE32;
-   int i, id = 0;
-   for (i=0;i<len;i++)
-   {
-      if (x[i] > max_corr)
-      {
-         id = i;
-         max_corr = x[i];
-      }
-   }
-   return id;
-}
-#endif
-
+/* Multiplies two 16-bit fractional values. Bit-exactness of this macro is important */
 #define FRAC_MUL16(a,b) ((16384+((celt_int32)(celt_int16)(a)*(celt_int16)(b)))>>15)
-static inline celt_int16 bitexact_cos(celt_int16 x)
-{
-   celt_int32 tmp;
-   celt_int16 x2;
-   tmp = (4096+((celt_int32)(x)*(x)))>>13;
-   if (tmp > 32767)
-      tmp = 32767;
-   x2 = tmp;
-   x2 = (32767-x2) + FRAC_MUL16(x2, (-7651 + FRAC_MUL16(x2, (8277 + FRAC_MUL16(-626, x2)))));
-   if (x2 > 32766)
-      x2 = 32766;
-   return 1+x2;
-}
 
+unsigned isqrt32(celt_uint32 _val);
 
 #ifndef FIXED_POINT
 
 #define celt_sqrt(x) ((float)sqrt(x))
-#define celt_psqrt(x) ((float)sqrt(x))
 #define celt_rsqrt(x) (1.f/celt_sqrt(x))
+#define celt_rsqrt_norm(x) (celt_rsqrt(x))
 #define celt_acos acos
 #define celt_exp exp
-#define celt_cos_norm(x) (cos((.5f*M_PI)*(x)))
+#define celt_cos_norm(x) ((float)cos((.5f*M_PI)*(x)))
 #define celt_atan atan
 #define celt_rcp(x) (1.f/(x))
 #define celt_div(a,b) ((a)/(b))
+#define frac_div32(a,b) ((float)(a)/(b))
 
 #ifdef FLOAT_APPROX
 
@@ -157,8 +104,8 @@ static inline float celt_exp2(float x)
 }
 
 #else
-#define celt_log2(x) (1.442695040888963387*log(x))
-#define celt_exp2(x) (exp(0.6931471805599453094*(x)))
+#define celt_log2(x) ((float)(1.442695040888963387*log(x)))
+#define celt_exp2(x) ((float)exp(0.6931471805599453094*(x)))
 #endif
 
 #endif
@@ -169,6 +116,16 @@ static inline float celt_exp2(float x)
 
 #include "os_support.h"
 
+#ifndef OVERRIDE_CELT_ILOG2
+/** Integer log in base2. Undefined for zero and negative numbers */
+static inline celt_int16 celt_ilog2(celt_int32 x)
+{
+   celt_assert2(x>0, "celt_ilog2() only defined for strictly positive numbers");
+   return EC_ILOG(x)-1;
+}
+#endif
+
+
 #ifndef OVERRIDE_CELT_MAXABS16
 static inline celt_word16 celt_maxabs16(celt_word16 *x, int len)
 {
@@ -186,132 +143,26 @@ static inline celt_int16 celt_zlog2(celt_word32 x)
    return x <= 0 ? 0 : celt_ilog2(x);
 }
 
-/** Reciprocal sqrt approximation (Q30 input, Q0 output or equivalent) */
-static inline celt_word32 celt_rsqrt(celt_word32 x)
-{
-   int k;
-   celt_word16 n;
-   celt_word16 r;
-   celt_word16 r2;
-   celt_word16 y;
-   celt_word32 rt;
-   k = celt_ilog2(x)>>1;
-   x = VSHR32(x, (k-7)<<1);
-   /* Range of n is [-16384,32767] ([-0.5,1) in Q15). */
-   n = x-32768;
-   /* Get a rough initial guess for the root.
-      The optimal minimax quadratic approximation (using relative error) is
-       r = 1.437799046117536+n*(-0.823394375837328+n*0.4096419668459485).
-      Coefficients here, and the final result r, are Q14.*/
-   r = ADD16(23557, MULT16_16_Q15(n, ADD16(-13490, MULT16_16_Q15(n, 6713))));
-   /* We want y = x*r*r-1 in Q15, but x is 32-bit Q16 and r is Q14.
-      We can compute the result from n and r using Q15 multiplies with some
-       adjustment, carefully done to avoid overflow.
-      Range of y is [-1564,1594]. */
-   r2 = MULT16_16_Q15(r, r);
-   y = SHL16(SUB16(ADD16(MULT16_16_Q15(r2, n), r2), 16384), 1);
-   /* Apply a 2nd-order Householder iteration: r += r*y*(y*0.375-0.5). */
-   rt = ADD16(r, MULT16_16_Q15(r, MULT16_16_Q15(y,
-              SUB16(MULT16_16_Q15(y, 12288), 16384))));
-   /* rt is now the Q14 reciprocal square root of the Q16 x, with a maximum
-       relative error of 1.04956E-4, a (relative) RMSE of 2.80979E-5, and a
-       peak absolute error of 2.26591/16384. */
-   /* Most of the error in this function comes from the following shift: */
-   rt = PSHR32(rt,k);
-   return rt;
-}
-
-/** Sqrt approximation (QX input, QX/2 output) */
-static inline celt_word32 celt_sqrt(celt_word32 x)
-{
-   int k;
-   celt_word16 n;
-   celt_word32 rt;
-   const celt_word16 C[5] = {23175, 11561, -3011, 1699, -664};
-   if (x==0)
-      return 0;
-   k = (celt_ilog2(x)>>1)-7;
-   x = VSHR32(x, (k<<1));
-   n = x-32768;
-   rt = ADD16(C[0], MULT16_16_Q15(n, ADD16(C[1], MULT16_16_Q15(n, ADD16(C[2], 
-              MULT16_16_Q15(n, ADD16(C[3], MULT16_16_Q15(n, (C[4])))))))));
-   rt = VSHR32(rt,7-k);
-   return rt;
-}
-
-/** Sqrt approximation (QX input, QX/2 output) that assumes that the input is
-    strictly positive */
-static inline celt_word32 celt_psqrt(celt_word32 x)
-{
-   int k;
-   celt_word16 n;
-   celt_word32 rt;
-   const celt_word16 C[5] = {23175, 11561, -3011, 1699, -664};
-   k = (celt_ilog2(x)>>1)-7;
-   x = VSHR32(x, (k<<1));
-   n = x-32768;
-   rt = ADD16(C[0], MULT16_16_Q15(n, ADD16(C[1], MULT16_16_Q15(n, ADD16(C[2], 
-              MULT16_16_Q15(n, ADD16(C[3], MULT16_16_Q15(n, (C[4])))))))));
-   rt = VSHR32(rt,7-k);
-   return rt;
-}
-
-#define L1 32767
-#define L2 -7651
-#define L3 8277
-#define L4 -626
+celt_word16 celt_rsqrt_norm(celt_word32 x);
 
-static inline celt_word16 _celt_cos_pi_2(celt_word16 x)
-{
-   celt_word16 x2;
-   
-   x2 = MULT16_16_P15(x,x);
-   return ADD16(1,MIN16(32766,ADD32(SUB16(L1,x2), MULT16_16_P15(x2, ADD32(L2, MULT16_16_P15(x2, ADD32(L3, MULT16_16_P15(L4, x2
-                                                                                ))))))));
-}
+celt_word32 celt_sqrt(celt_word32 x);
 
-#undef L1
-#undef L2
-#undef L3
-#undef L4
+celt_word16 celt_cos_norm(celt_word32 x);
 
-static inline celt_word16 celt_cos_norm(celt_word32 x)
-{
-   x = x&0x0001ffff;
-   if (x>SHL32(EXTEND32(1), 16))
-      x = SUB32(SHL32(EXTEND32(1), 17),x);
-   if (x&0x00007fff)
-   {
-      if (x<SHL32(EXTEND32(1), 15))
-      {
-         return _celt_cos_pi_2(EXTRACT16(x));
-      } else {
-         return NEG32(_celt_cos_pi_2(EXTRACT16(65536-x)));
-      }
-   } else {
-      if (x&0x0000ffff)
-         return 0;
-      else if (x&0x0001ffff)
-         return -32767;
-      else
-         return 32767;
-   }
-}
 
 static inline celt_word16 celt_log2(celt_word32 x)
 {
    int i;
    celt_word16 n, frac;
-   /*-0.41446   0.96093  -0.33981   0.15600 */
-   /* -0.4144541824871411+32/16384, 0.9590923197873218, -0.3395129038105771,
-       0.16541096501128538 */
-   const celt_word16 C[4] = {-6758, 15715, -5563, 2708};
+   /* -0.41509302963303146, 0.9609890551383969, -0.31836011537636605,
+       0.15530808010959576, -0.08556153059057618 */
+   static const celt_word16 C[5] = {-6801+(1<<13-DB_SHIFT), 15746, -5217, 2545, -1401};
    if (x==0)
       return -32767;
    i = celt_ilog2(x);
    n = VSHR32(x,i-15)-32768-16384;
-   frac = ADD16(C[0], MULT16_16_Q15(n, ADD16(C[1], MULT16_16_Q15(n, ADD16(C[2], MULT16_16_Q15(n, (C[3])))))));
-   return SHL16(i-13,8)+SHR16(frac,14-8);
+   frac = ADD16(C[0], MULT16_16_Q15(n, ADD16(C[1], MULT16_16_Q15(n, ADD16(C[2], MULT16_16_Q15(n, ADD16(C[3], MULT16_16_Q15(n, C[4]))))))));
+   return SHL16(i-13,DB_SHIFT)+SHR16(frac,14-DB_SHIFT);
 }
 
 /*
@@ -339,43 +190,19 @@ static inline celt_word32 celt_exp2(celt_word16 x)
    return VSHR32(EXTEND32(frac), -integer-2);
 }
 
-/** Reciprocal approximation (Q15 input, Q16 output) */
-static inline celt_word32 celt_rcp(celt_word32 x)
-{
-   int i;
-   celt_word16 n;
-   celt_word16 r;
-   celt_assert2(x>0, "celt_rcp() only defined for positive values");
-   i = celt_ilog2(x);
-   /* n is Q15 with range [0,1). */
-   n = VSHR32(x,i-15)-32768;
-   /* Start with a linear approximation:
-      r = 1.8823529411764706-0.9411764705882353*n.
-      The coefficients and the result are Q14 in the range [15420,30840].*/
-   r = ADD16(30840, MULT16_16_Q15(-15420, n));
-   /* Perform two Newton iterations:
-      r -= r*((r*n)-1.Q15)
-         = r*((r*n)+(r-1.Q15)). */
-   r = SUB16(r, MULT16_16_Q15(r,
-             ADD16(MULT16_16_Q15(r, n), ADD16(r, -32768))));
-   /* We subtract an extra 1 in the second iteration to avoid overflow; it also
-       neatly compensates for truncation error in the rest of the process. */
-   r = SUB16(r, ADD16(1, MULT16_16_Q15(r,
-             ADD16(MULT16_16_Q15(r, n), ADD16(r, -32768)))));
-   /* r is now the Q15 solution to 2/(n+1), with a maximum relative error
-       of 7.05346E-5, a (relative) RMSE of 2.14418E-5, and a peak absolute
-       error of 1.24665/32768. */
-   return VSHR32(EXTEND32(r),i-16);
-}
+celt_word32 celt_rcp(celt_word32 x);
 
 #define celt_div(a,b) MULT32_32_Q31((celt_word32)(a),celt_rcp(b))
 
+celt_word32 frac_div32(celt_word32 a, celt_word32 b);
 
 #define M1 32767
 #define M2 -21
 #define M3 -11943
 #define M4 4936
 
+/* Atan approximation using a 4th order polynomial. Input is in Q15 format
+   and normalized by pi/4. Output is in Q15 format */
 static inline celt_word16 celt_atan01(celt_word16 x)
 {
    return MULT16_16_P15(x, ADD32(M1, MULT16_16_P15(x, ADD32(M2, MULT16_16_P15(x, ADD32(M3, MULT16_16_P15(M4, x)))))));
@@ -386,6 +213,7 @@ static inline celt_word16 celt_atan01(celt_word16 x)
 #undef M3
 #undef M4
 
+/* atan2() approximation valid for positive input values */
 static inline celt_word16 celt_atan2p(celt_word16 y, celt_word16 x)
 {
    if (y < x)