1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/media/libopus/celt/mathops.h Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,258 @@
1.4 +/* Copyright (c) 2002-2008 Jean-Marc Valin
1.5 + Copyright (c) 2007-2008 CSIRO
1.6 + Copyright (c) 2007-2009 Xiph.Org Foundation
1.7 + Written by Jean-Marc Valin */
1.8 +/**
1.9 + @file mathops.h
1.10 + @brief Various math functions
1.11 +*/
1.12 +/*
1.13 + Redistribution and use in source and binary forms, with or without
1.14 + modification, are permitted provided that the following conditions
1.15 + are met:
1.16 +
1.17 + - Redistributions of source code must retain the above copyright
1.18 + notice, this list of conditions and the following disclaimer.
1.19 +
1.20 + - Redistributions in binary form must reproduce the above copyright
1.21 + notice, this list of conditions and the following disclaimer in the
1.22 + documentation and/or other materials provided with the distribution.
1.23 +
1.24 + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
1.25 + ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
1.26 + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
1.27 + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
1.28 + OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
1.29 + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
1.30 + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
1.31 + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
1.32 + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
1.33 + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
1.34 + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
1.35 +*/
1.36 +
1.37 +#ifndef MATHOPS_H
1.38 +#define MATHOPS_H
1.39 +
1.40 +#include "arch.h"
1.41 +#include "entcode.h"
1.42 +#include "os_support.h"
1.43 +
1.44 +/* Multiplies two 16-bit fractional values. Bit-exactness of this macro is important */
1.45 +#define FRAC_MUL16(a,b) ((16384+((opus_int32)(opus_int16)(a)*(opus_int16)(b)))>>15)
1.46 +
1.47 +unsigned isqrt32(opus_uint32 _val);
1.48 +
1.49 +#ifndef OVERRIDE_CELT_MAXABS16
1.50 +static OPUS_INLINE opus_val32 celt_maxabs16(const opus_val16 *x, int len)
1.51 +{
1.52 + int i;
1.53 + opus_val16 maxval = 0;
1.54 + opus_val16 minval = 0;
1.55 + for (i=0;i<len;i++)
1.56 + {
1.57 + maxval = MAX16(maxval, x[i]);
1.58 + minval = MIN16(minval, x[i]);
1.59 + }
1.60 + return MAX32(EXTEND32(maxval),-EXTEND32(minval));
1.61 +}
1.62 +#endif
1.63 +
1.64 +#ifndef OVERRIDE_CELT_MAXABS32
1.65 +#ifdef FIXED_POINT
1.66 +static OPUS_INLINE opus_val32 celt_maxabs32(const opus_val32 *x, int len)
1.67 +{
1.68 + int i;
1.69 + opus_val32 maxval = 0;
1.70 + opus_val32 minval = 0;
1.71 + for (i=0;i<len;i++)
1.72 + {
1.73 + maxval = MAX32(maxval, x[i]);
1.74 + minval = MIN32(minval, x[i]);
1.75 + }
1.76 + return MAX32(maxval, -minval);
1.77 +}
1.78 +#else
1.79 +#define celt_maxabs32(x,len) celt_maxabs16(x,len)
1.80 +#endif
1.81 +#endif
1.82 +
1.83 +
1.84 +#ifndef FIXED_POINT
1.85 +
1.86 +#define PI 3.141592653f
1.87 +#define celt_sqrt(x) ((float)sqrt(x))
1.88 +#define celt_rsqrt(x) (1.f/celt_sqrt(x))
1.89 +#define celt_rsqrt_norm(x) (celt_rsqrt(x))
1.90 +#define celt_cos_norm(x) ((float)cos((.5f*PI)*(x)))
1.91 +#define celt_rcp(x) (1.f/(x))
1.92 +#define celt_div(a,b) ((a)/(b))
1.93 +#define frac_div32(a,b) ((float)(a)/(b))
1.94 +
1.95 +#ifdef FLOAT_APPROX
1.96 +
1.97 +/* Note: This assumes radix-2 floating point with the exponent at bits 23..30 and an offset of 127
1.98 + denorm, +/- inf and NaN are *not* handled */
1.99 +
1.100 +/** Base-2 log approximation (log2(x)). */
1.101 +static OPUS_INLINE float celt_log2(float x)
1.102 +{
1.103 + int integer;
1.104 + float frac;
1.105 + union {
1.106 + float f;
1.107 + opus_uint32 i;
1.108 + } in;
1.109 + in.f = x;
1.110 + integer = (in.i>>23)-127;
1.111 + in.i -= integer<<23;
1.112 + frac = in.f - 1.5f;
1.113 + frac = -0.41445418f + frac*(0.95909232f
1.114 + + frac*(-0.33951290f + frac*0.16541097f));
1.115 + return 1+integer+frac;
1.116 +}
1.117 +
1.118 +/** Base-2 exponential approximation (2^x). */
1.119 +static OPUS_INLINE float celt_exp2(float x)
1.120 +{
1.121 + int integer;
1.122 + float frac;
1.123 + union {
1.124 + float f;
1.125 + opus_uint32 i;
1.126 + } res;
1.127 + integer = floor(x);
1.128 + if (integer < -50)
1.129 + return 0;
1.130 + frac = x-integer;
1.131 + /* K0 = 1, K1 = log(2), K2 = 3-4*log(2), K3 = 3*log(2) - 2 */
1.132 + res.f = 0.99992522f + frac * (0.69583354f
1.133 + + frac * (0.22606716f + 0.078024523f*frac));
1.134 + res.i = (res.i + (integer<<23)) & 0x7fffffff;
1.135 + return res.f;
1.136 +}
1.137 +
1.138 +#else
1.139 +#define celt_log2(x) ((float)(1.442695040888963387*log(x)))
1.140 +#define celt_exp2(x) ((float)exp(0.6931471805599453094*(x)))
1.141 +#endif
1.142 +
1.143 +#endif
1.144 +
1.145 +#ifdef FIXED_POINT
1.146 +
1.147 +#include "os_support.h"
1.148 +
1.149 +#ifndef OVERRIDE_CELT_ILOG2
1.150 +/** Integer log in base2. Undefined for zero and negative numbers */
1.151 +static OPUS_INLINE opus_int16 celt_ilog2(opus_int32 x)
1.152 +{
1.153 + celt_assert2(x>0, "celt_ilog2() only defined for strictly positive numbers");
1.154 + return EC_ILOG(x)-1;
1.155 +}
1.156 +#endif
1.157 +
1.158 +
1.159 +/** Integer log in base2. Defined for zero, but not for negative numbers */
1.160 +static OPUS_INLINE opus_int16 celt_zlog2(opus_val32 x)
1.161 +{
1.162 + return x <= 0 ? 0 : celt_ilog2(x);
1.163 +}
1.164 +
1.165 +opus_val16 celt_rsqrt_norm(opus_val32 x);
1.166 +
1.167 +opus_val32 celt_sqrt(opus_val32 x);
1.168 +
1.169 +opus_val16 celt_cos_norm(opus_val32 x);
1.170 +
1.171 +/** Base-2 logarithm approximation (log2(x)). (Q14 input, Q10 output) */
1.172 +static OPUS_INLINE opus_val16 celt_log2(opus_val32 x)
1.173 +{
1.174 + int i;
1.175 + opus_val16 n, frac;
1.176 + /* -0.41509302963303146, 0.9609890551383969, -0.31836011537636605,
1.177 + 0.15530808010959576, -0.08556153059057618 */
1.178 + static const opus_val16 C[5] = {-6801+(1<<(13-DB_SHIFT)), 15746, -5217, 2545, -1401};
1.179 + if (x==0)
1.180 + return -32767;
1.181 + i = celt_ilog2(x);
1.182 + n = VSHR32(x,i-15)-32768-16384;
1.183 + 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]))))))));
1.184 + return SHL16(i-13,DB_SHIFT)+SHR16(frac,14-DB_SHIFT);
1.185 +}
1.186 +
1.187 +/*
1.188 + K0 = 1
1.189 + K1 = log(2)
1.190 + K2 = 3-4*log(2)
1.191 + K3 = 3*log(2) - 2
1.192 +*/
1.193 +#define D0 16383
1.194 +#define D1 22804
1.195 +#define D2 14819
1.196 +#define D3 10204
1.197 +
1.198 +static OPUS_INLINE opus_val32 celt_exp2_frac(opus_val16 x)
1.199 +{
1.200 + opus_val16 frac;
1.201 + frac = SHL16(x, 4);
1.202 + return ADD16(D0, MULT16_16_Q15(frac, ADD16(D1, MULT16_16_Q15(frac, ADD16(D2 , MULT16_16_Q15(D3,frac))))));
1.203 +}
1.204 +/** Base-2 exponential approximation (2^x). (Q10 input, Q16 output) */
1.205 +static OPUS_INLINE opus_val32 celt_exp2(opus_val16 x)
1.206 +{
1.207 + int integer;
1.208 + opus_val16 frac;
1.209 + integer = SHR16(x,10);
1.210 + if (integer>14)
1.211 + return 0x7f000000;
1.212 + else if (integer < -15)
1.213 + return 0;
1.214 + frac = celt_exp2_frac(x-SHL16(integer,10));
1.215 + return VSHR32(EXTEND32(frac), -integer-2);
1.216 +}
1.217 +
1.218 +opus_val32 celt_rcp(opus_val32 x);
1.219 +
1.220 +#define celt_div(a,b) MULT32_32_Q31((opus_val32)(a),celt_rcp(b))
1.221 +
1.222 +opus_val32 frac_div32(opus_val32 a, opus_val32 b);
1.223 +
1.224 +#define M1 32767
1.225 +#define M2 -21
1.226 +#define M3 -11943
1.227 +#define M4 4936
1.228 +
1.229 +/* Atan approximation using a 4th order polynomial. Input is in Q15 format
1.230 + and normalized by pi/4. Output is in Q15 format */
1.231 +static OPUS_INLINE opus_val16 celt_atan01(opus_val16 x)
1.232 +{
1.233 + return MULT16_16_P15(x, ADD32(M1, MULT16_16_P15(x, ADD32(M2, MULT16_16_P15(x, ADD32(M3, MULT16_16_P15(M4, x)))))));
1.234 +}
1.235 +
1.236 +#undef M1
1.237 +#undef M2
1.238 +#undef M3
1.239 +#undef M4
1.240 +
1.241 +/* atan2() approximation valid for positive input values */
1.242 +static OPUS_INLINE opus_val16 celt_atan2p(opus_val16 y, opus_val16 x)
1.243 +{
1.244 + if (y < x)
1.245 + {
1.246 + opus_val32 arg;
1.247 + arg = celt_div(SHL32(EXTEND32(y),15),x);
1.248 + if (arg >= 32767)
1.249 + arg = 32767;
1.250 + return SHR16(celt_atan01(EXTRACT16(arg)),1);
1.251 + } else {
1.252 + opus_val32 arg;
1.253 + arg = celt_div(SHL32(EXTEND32(x),15),y);
1.254 + if (arg >= 32767)
1.255 + arg = 32767;
1.256 + return 25736-SHR16(celt_atan01(EXTRACT16(arg)),1);
1.257 + }
1.258 +}
1.259 +
1.260 +#endif /* FIXED_POINT */
1.261 +#endif /* MATHOPS_H */
