Wrapping all allocation within opus_alloc() and opus_free()
[opus.git] / libcelt / mathops.h
index 2218429..896fc77 100644 (file)
@@ -1,4 +1,7 @@
-/* Copyright (C) 2002-2008 Jean-Marc Valin */
+/* Copyright (c) 2002-2008 Jean-Marc Valin
+   Copyright (c) 2007-2008 CSIRO
+   Copyright (c) 2007-2009 Xiph.Org Foundation
+   Written by Jean-Marc Valin */
 /**
    @file mathops.h
    @brief Various math functions
    Redistribution and use in source and binary forms, with or without
    modification, are permitted provided that the following conditions
    are met:
-   
+
    - Redistributions of source code must retain the above copyright
    notice, this list of conditions and the following disclaimer.
-   
+
    - Redistributions in binary form must reproduce the above copyright
    notice, this list of conditions and the following disclaimer in the
    documentation and/or other materials provided with the distribution.
-   
-   - Neither the name of the Xiph.org Foundation nor the names of its
-   contributors may be used to endorse or promote products derived from
-   this software without specific prior written permission.
-   
+
    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 #define MATHOPS_H
 
 #include "arch.h"
+#include "entcode.h"
+#include "os_support.h"
 
-#ifndef OVERRIDE_FIND_MAX16
-static inline int find_max16(celt_word16_t *x, int len)
-{
-   celt_word16_t 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_t *x, int len)
-{
-   celt_word32_t 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+((opus_int32)(opus_int16)(a)*(opus_int16)(b)))>>15)
 
+unsigned isqrt32(opus_uint32 _val);
 
 #ifndef FIXED_POINT
 
+#define PI 3.141592653f
 #define celt_sqrt(x) ((float)sqrt(x))
 #define celt_rsqrt(x) (1.f/celt_sqrt(x))
-#define celt_acos acos
+#define celt_rsqrt_norm(x) (celt_rsqrt(x))
 #define celt_exp exp
-#define celt_cos_norm(x) (cos((.5f*M_PI)*(x)))
+#define celt_cos_norm(x) ((float)cos((.5f*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))
 
-#endif
+#ifdef FLOAT_APPROX
+
+/* Note: This assumes radix-2 floating point with the exponent at bits 23..30 and an offset of 127
+         denorm, +/- inf and NaN are *not* handled */
+
+/** Base-2 log approximation (log2(x)). */
+static inline float celt_log2(float x)
+{
+   int integer;
+   float frac;
+   union {
+      float f;
+      opus_uint32 i;
+   } in;
+   in.f = x;
+   integer = (in.i>>23)-127;
+   in.i -= integer<<23;
+   frac = in.f - 1.5f;
+   frac = -0.41445418f + frac*(0.95909232f
+          + frac*(-0.33951290f + frac*0.16541097f));
+   return 1+integer+frac;
+}
 
+/** Base-2 exponential approximation (2^x). */
+static inline float celt_exp2(float x)
+{
+   int integer;
+   float frac;
+   union {
+      float f;
+      opus_uint32 i;
+   } res;
+   integer = floor(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));
+   res.i = (res.i + (integer<<23)) & 0x7fffffff;
+   return res.f;
+}
 
+#else
+#define celt_log2(x) ((float)(1.442695040888963387*log(x)))
+#define celt_exp2(x) ((float)exp(0.6931471805599453094*(x)))
+#endif
+
+#endif
 
 #ifdef FIXED_POINT
 
-#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_t celt_ilog2(celt_word32_t x)
+static inline opus_int16 celt_ilog2(opus_int32 x)
 {
    celt_assert2(x>0, "celt_ilog2() only defined for strictly positive numbers");
    return EC_ILOG(x)-1;
@@ -102,10 +120,10 @@ static inline celt_int16_t celt_ilog2(celt_word32_t x)
 #endif
 
 #ifndef OVERRIDE_CELT_MAXABS16
-static inline celt_word16_t celt_maxabs16(celt_word16_t *x, int len)
+static inline opus_val16 celt_maxabs16(opus_val16 *x, int len)
 {
    int i;
-   celt_word16_t maxval = 0;
+   opus_val16 maxval = 0;
    for (i=0;i<len;i++)
       maxval = MAX16(maxval, ABS16(x[i]));
    return maxval;
@@ -113,102 +131,30 @@ static inline celt_word16_t celt_maxabs16(celt_word16_t *x, int len)
 #endif
 
 /** Integer log in base2. Defined for zero, but not for negative numbers */
-static inline celt_int16_t celt_zlog2(celt_word32_t x)
+static inline opus_int16 celt_zlog2(opus_val32 x)
 {
-   return EC_ILOG(x)-1;
-}
-
-/** Reciprocal sqrt approximation (Q30 input, Q0 output or equivalent) */
-static inline celt_word32_t celt_rsqrt(celt_word32_t x)
-{
-   int k;
-   celt_word16_t n;
-   celt_word32_t rt;
-   const celt_word16_t C[5] = {23126, -11496, 9812, -9097, 4100};
-   k = celt_ilog2(x)>>1;
-   x = VSHR32(x, (k-7)<<1);
-   /* Range of n is [-16384,32767] */
-   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,k);
-   return rt;
-}
-
-/** Sqrt approximation (QX input, QX/2 output) */
-static inline celt_word32_t celt_sqrt(celt_word32_t x)
-{
-   int k;
-   celt_word16_t n;
-   celt_word32_t rt;
-   const celt_word16_t C[5] = {23174, 11584, -3011, 1570, -557};
-   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;
+   return x <= 0 ? 0 : celt_ilog2(x);
 }
 
+opus_val16 celt_rsqrt_norm(opus_val32 x);
 
-#define L1 32767
-#define L2 -7651
-#define L3 8277
-#define L4 -626
-
-static inline celt_word16_t _celt_cos_pi_2(celt_word16_t x)
-{
-   celt_word16_t 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
-                                                                                ))))))));
-}
-
-#undef L1
-#undef L2
-#undef L3
-#undef L4
+opus_val32 celt_sqrt(opus_val32 x);
 
-static inline celt_word16_t celt_cos_norm(celt_word32_t 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;
-   }
-}
+opus_val16 celt_cos_norm(opus_val32 x);
 
-static inline celt_word16_t celt_log2(celt_word32_t x)
+static inline opus_val16 celt_log2(opus_val32 x)
 {
    int i;
-   celt_word16_t n, frac;
-   /*-0.41446   0.96093  -0.33981   0.15600 */
-   const celt_word16_t C[4] = {-6791, 7872, -1392, 319};
+   opus_val16 n, frac;
+   /* -0.41509302963303146, 0.9609890551383969, -0.31836011537636605,
+       0.15530808010959576, -0.08556153059057618 */
+   static const opus_val16 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_Q14(n, ADD16(C[1], MULT16_16_Q14(n, ADD16(C[2], MULT16_16_Q14(n, (C[3])))))));
-   /*printf ("%d %d %d %d\n", x, n, ret, SHL16(i-13,8)+SHR16(ret,14-8));*/
-   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);
 }
 
 /*
@@ -217,42 +163,66 @@ static inline celt_word16_t celt_log2(celt_word32_t x)
  K2 = 3-4*log(2)
  K3 = 3*log(2) - 2
 */
-#define D0 16384
-#define D1 11356
-#define D2 3726
-#define D3 1301
-/** Base-2 exponential approximation (2^x). (Q11 input, Q16 output) */
-static inline celt_word32_t celt_exp2(celt_word16_t x)
+#define D0 16383
+#define D1 22804
+#define D2 14819
+#define D3 10204
+/** Base-2 exponential approximation (2^x). (Q10 input, Q16 output) */
+static inline opus_val32 celt_exp2(opus_val16 x)
 {
    int integer;
-   celt_word16_t frac;
-   integer = SHR16(x,11);
+   opus_val16 frac;
+   integer = SHR16(x,10);
    if (integer>14)
-      return 0x7fffffff;
+      return 0x7f000000;
    else if (integer < -15)
       return 0;
-   frac = SHL16(x-SHL16(integer,11),3);
-   frac = ADD16(D0, MULT16_16_Q14(frac, ADD16(D1, MULT16_16_Q14(frac, ADD16(D2 , MULT16_16_Q14(D3,frac))))));
+   frac = SHL16(x-SHL16(integer,10),4);
+   frac = ADD16(D0, MULT16_16_Q15(frac, ADD16(D1, MULT16_16_Q15(frac, ADD16(D2 , MULT16_16_Q15(D3,frac))))));
    return VSHR32(EXTEND32(frac), -integer-2);
 }
 
-/** Reciprocal approximation (Q15 input, Q16 output) */
-static inline celt_word32_t celt_rcp(celt_word32_t x)
+opus_val32 celt_rcp(opus_val32 x);
+
+#define celt_div(a,b) MULT32_32_Q31((opus_val32)(a),celt_rcp(b))
+
+opus_val32 frac_div32(opus_val32 a, opus_val32 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 opus_val16 celt_atan01(opus_val16 x)
 {
-   int i;
-   celt_word16_t n, frac;
-   const celt_word16_t C[5] = {21848, -7251, 2403, -934, 327};
-   celt_assert2(x>0, "celt_rcp() only defined for positive values");
-   i = celt_ilog2(x);
-   n = VSHR32(x,i-16)-SHL32(EXTEND32(3),15);
-   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 VSHR32(EXTEND32(frac),i-16);
+   return MULT16_16_P15(x, ADD32(M1, MULT16_16_P15(x, ADD32(M2, MULT16_16_P15(x, ADD32(M3, MULT16_16_P15(M4, x)))))));
 }
 
-#define celt_div(a,b) MULT32_32_Q31((celt_word32_t)(a),celt_rcp(b))
-
-#endif /* FIXED_POINT */
+#undef M1
+#undef M2
+#undef M3
+#undef M4
 
+/* atan2() approximation valid for positive input values */
+static inline opus_val16 celt_atan2p(opus_val16 y, opus_val16 x)
+{
+   if (y < x)
+   {
+      opus_val32 arg;
+      arg = celt_div(SHL32(EXTEND32(y),15),x);
+      if (arg >= 32767)
+         arg = 32767;
+      return SHR16(celt_atan01(EXTRACT16(arg)),1);
+   } else {
+      opus_val32 arg;
+      arg = celt_div(SHL32(EXTEND32(x),15),y);
+      if (arg >= 32767)
+         arg = 32767;
+      return 25736-SHR16(celt_atan01(EXTRACT16(arg)),1);
+   }
+}
 
+#endif /* FIXED_POINT */
 #endif /* MATHOPS_H */