/src/opus/celt/mathops.h

Source (jump to first uncovered line)
/* 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.

   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
   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#ifndef MATHOPS_H
#define MATHOPS_H

#include "arch.h"
#include "entcode.h"
#include "os_support.h"

#define PI 3.141592653f

/* 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);

/* CELT doesn't need it for fixed-point, by analysis.c does. */
#if !defined(FIXED_POINT) || defined(ANALYSIS_C)
#define cA 0.43157974f
#define cB 0.67848403f
#define cC 0.08595542f
#define cE ((float)PI/2)
static OPUS_INLINE float fast_atan2f(float y, float x) {
   float x2, y2;
   x2 = x*x;
   y2 = y*y;
   /* For very small values, we don't care about the answer, so
      we can just return 0. */
   if (x2 + y2 < 1e-18f)
   {
      return 0;
   }
   if(x2<y2){
      float den = (y2 + cB*x2) * (y2 + cC*x2);
      return -x*y*(y2 + cA*x2) / den + (y<0 ? -cE : cE);
   }else{
      float den = (x2 + cB*y2) * (x2 + cC*y2);
      return  x*y*(x2 + cA*y2) / den + (y<0 ? -cE : cE) - (x*y<0 ? -cE : cE);
   }
}
#undef cA
#undef cB
#undef cC
#undef cE
#endif


#ifndef OVERRIDE_CELT_MAXABS16
static OPUS_INLINE opus_val32 celt_maxabs16(const opus_val16 *x, int len)
{
   int i;
   opus_val16 maxval = 0;
   opus_val16 minval = 0;
   for (i=0;i<len;i++)
   {
      maxval = MAX16(maxval, x[i]);
      minval = MIN16(minval, x[i]);
   }
   return MAX32(EXTEND32(maxval),-EXTEND32(minval));
}
#endif

#ifndef OVERRIDE_CELT_MAXABS32
#ifdef FIXED_POINT
static OPUS_INLINE opus_val32 celt_maxabs32(const opus_val32 *x, int len)
{
   int i;
   opus_val32 maxval = 0;
   opus_val32 minval = 0;
   for (i=0;i<len;i++)
   {
      maxval = MAX32(maxval, x[i]);
      minval = MIN32(minval, x[i]);
   }
   return MAX32(maxval, -minval);
}
#else
#define celt_maxabs32(x,len) celt_maxabs16(x,len)
#endif
#endif


#ifndef FIXED_POINT

#define celt_sqrt(x) ((float)sqrt(x))
#define celt_rsqrt(x) (1.f/celt_sqrt(x))
#define celt_rsqrt_norm(x) (celt_rsqrt(x))
#define celt_cos_norm(x) ((float)cos((.5f*PI)*(x)))
#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

/* 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 OPUS_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 -= (opus_uint32)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 OPUS_INLINE float celt_exp2(float x)
{
   int integer;
   float frac;
   union {
      float f;
      opus_uint32 i;
   } res;
   integer = (int)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 + ((opus_uint32)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 "os_support.h"

#ifndef OVERRIDE_CELT_ILOG2
/** Integer log in base2. Undefined for zero and negative numbers */
static OPUS_INLINE opus_int16 celt_ilog2(opus_int32 x)
{
   celt_sig_assert(x>0);
   return EC_ILOG(x)-1;
}
#endif


/** Integer log in base2. Defined for zero, but not for negative numbers */
static OPUS_INLINE opus_int16 celt_zlog2(opus_val32 x)
{
   return x <= 0 ? 0 : celt_ilog2(x);
}

opus_val16 celt_rsqrt_norm(opus_val32 x);

opus_val32 celt_sqrt(opus_val32 x);

opus_val16 celt_cos_norm(opus_val32 x);

/** Base-2 logarithm approximation (log2(x)). (Q14 input, Q10 output) */
static OPUS_INLINE opus_val16 celt_log2(opus_val32 x)
{
   int i;
   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_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);
}

/*
 K0 = 1
 K1 = log(2)
 K2 = 3-4*log(2)
 K3 = 3*log(2) - 2
*/
#define D0 16383
#define D1 22804
#define D2 14819
#define D3 10204

static OPUS_INLINE opus_val32 celt_exp2_frac(opus_val16 x)
{
   opus_val16 frac;
   frac = SHL16(x, 4);
   return ADD16(D0, MULT16_16_Q15(frac, ADD16(D1, MULT16_16_Q15(frac, ADD16(D2 , MULT16_16_Q15(D3,frac))))));
}

#undef D0
#undef D1
#undef D2
#undef D3

/** Base-2 exponential approximation (2^x). (Q10 input, Q16 output) */
static OPUS_INLINE opus_val32 celt_exp2(opus_val16 x)
{
   int integer;
   opus_val16 frac;
   integer = SHR16(x,10);
   if (integer>14)
      return 0x7f000000;
   else if (integer < -15)
      return 0;
   frac = celt_exp2_frac(x-SHL16(integer,10));
   return VSHR32(EXTEND32(frac), -integer-2);
}

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 OPUS_INLINE opus_val16 celt_atan01(opus_val16 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)))))));
}

#undef M1
#undef M2
#undef M3
#undef M4

/* atan2() approximation valid for positive input values */
static OPUS_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 */

Coverage Report

Created: 2024-09-06 07:53

Line	Count	Source (jump to first uncovered line)
1		/* Copyright (c) 2002-2008 Jean-Marc Valin
2		Copyright (c) 2007-2008 CSIRO
3		Copyright (c) 2007-2009 Xiph.Org Foundation
4		Written by Jean-Marc Valin */
5		/**
6		@file mathops.h
7		@brief Various math functions
8		*/
9		/*
10		Redistribution and use in source and binary forms, with or without
11		modification, are permitted provided that the following conditions
12		are met:
13
14		- Redistributions of source code must retain the above copyright
15		notice, this list of conditions and the following disclaimer.
16
17		- Redistributions in binary form must reproduce the above copyright
18		notice, this list of conditions and the following disclaimer in the
19		documentation and/or other materials provided with the distribution.
20
21		THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22		``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23		LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24		A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
25		OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
26		EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
27		PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
28		PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
29		LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
30		NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31		SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32		*/
33
34		#ifndef MATHOPS_H
35		#define MATHOPS_H
36
37		#include "arch.h"
38		#include "entcode.h"
39		#include "os_support.h"
40
41	0	#define PI 3.141592653f
42
43		/* Multiplies two 16-bit fractional values. Bit-exactness of this macro is important */
44	0	#define FRAC_MUL16(a,b) ((16384+((opus_int32)(opus_int16)(a)*(opus_int16)(b)))>>15)
45
46		unsigned isqrt32(opus_uint32 _val);
47
48		/* CELT doesn't need it for fixed-point, by analysis.c does. */
49		#if !defined(FIXED_POINT) \|\| defined(ANALYSIS_C)
50	0	#define cA 0.43157974f
51	0	#define cB 0.67848403f
52	0	#define cC 0.08595542f
53	0	#define cE ((float)PI/2)
54	0	static OPUS_INLINE float fast_atan2f(float y, float x) {
55	0	float x2, y2;
56	0	x2 = x*x;
57	0	y2 = y*y;
58		/* For very small values, we don't care about the answer, so
59		we can just return 0. */
60	0	if (x2 + y2 < 1e-18f)
61	0	{
62	0	return 0;
63	0	}
64	0	if(x2<y2){
65	0	float den = (y2 + cBx2) (y2 + cC*x2);
66	0	return -xy(y2 + cA*x2) / den + (y<0 ? -cE : cE);
67	0	}else{
68	0	float den = (x2 + cBy2) (x2 + cC*y2);
69	0	return xy(x2 + cAy2) / den + (y<0 ? -cE : cE) - (xy<0 ? -cE : cE);
70	0	}
71	0	} Unexecuted instantiation: celt.c:fast_atan2f Unexecuted instantiation: pitch_sse.c:fast_atan2f Unexecuted instantiation: opus_multistream_encoder.c:fast_atan2f Unexecuted instantiation: bands.c:fast_atan2f Unexecuted instantiation: celt_encoder.c:fast_atan2f Unexecuted instantiation: mathops.c:fast_atan2f Unexecuted instantiation: mdct.c:fast_atan2f Unexecuted instantiation: modes.c:fast_atan2f Unexecuted instantiation: pitch.c:fast_atan2f Unexecuted instantiation: celt_lpc.c:fast_atan2f Unexecuted instantiation: quant_bands.c:fast_atan2f Unexecuted instantiation: vq.c:fast_atan2f Unexecuted instantiation: vq_sse2.c:fast_atan2f Unexecuted instantiation: opus_decoder.c:fast_atan2f Unexecuted instantiation: opus_encoder.c:fast_atan2f Unexecuted instantiation: analysis.c:fast_atan2f Unexecuted instantiation: celt_decoder.c:fast_atan2f Unexecuted instantiation: cwrs.c:fast_atan2f Unexecuted instantiation: kiss_fft.c:fast_atan2f Unexecuted instantiation: laplace.c:fast_atan2f
72		#undef cA
73		#undef cB
74		#undef cC
75		#undef cE
76		#endif
77
78
79		#ifndef OVERRIDE_CELT_MAXABS16
80		static OPUS_INLINE opus_val32 celt_maxabs16(const opus_val16 *x, int len)
81	0	{
82	0	int i;
83	0	opus_val16 maxval = 0;
84	0	opus_val16 minval = 0;
85	0	for (i=0;i<len;i++)
86	0	{
87	0	maxval = MAX16(maxval, x[i]);
88	0	minval = MIN16(minval, x[i]);
89	0	}
90	0	return MAX32(EXTEND32(maxval),-EXTEND32(minval));
91	0	} Unexecuted instantiation: celt.c:celt_maxabs16 Unexecuted instantiation: pitch_sse.c:celt_maxabs16 Unexecuted instantiation: opus_multistream_encoder.c:celt_maxabs16 Unexecuted instantiation: bands.c:celt_maxabs16 Unexecuted instantiation: celt_encoder.c:celt_maxabs16 Unexecuted instantiation: mathops.c:celt_maxabs16 Unexecuted instantiation: mdct.c:celt_maxabs16 Unexecuted instantiation: modes.c:celt_maxabs16 Unexecuted instantiation: pitch.c:celt_maxabs16 Unexecuted instantiation: celt_lpc.c:celt_maxabs16 Unexecuted instantiation: quant_bands.c:celt_maxabs16 Unexecuted instantiation: vq.c:celt_maxabs16 Unexecuted instantiation: vq_sse2.c:celt_maxabs16 Unexecuted instantiation: opus_decoder.c:celt_maxabs16 Unexecuted instantiation: opus_encoder.c:celt_maxabs16 Unexecuted instantiation: analysis.c:celt_maxabs16 Unexecuted instantiation: celt_decoder.c:celt_maxabs16 Unexecuted instantiation: cwrs.c:celt_maxabs16 Unexecuted instantiation: kiss_fft.c:celt_maxabs16 Unexecuted instantiation: laplace.c:celt_maxabs16
92		#endif
93
94		#ifndef OVERRIDE_CELT_MAXABS32
95		#ifdef FIXED_POINT
96		static OPUS_INLINE opus_val32 celt_maxabs32(const opus_val32 *x, int len)
97		{
98		int i;
99		opus_val32 maxval = 0;
100		opus_val32 minval = 0;
101		for (i=0;i<len;i++)
102		{
103		maxval = MAX32(maxval, x[i]);
104		minval = MIN32(minval, x[i]);
105		}
106		return MAX32(maxval, -minval);
107		}
108		#else
109		#define celt_maxabs32(x,len) celt_maxabs16(x,len)
110		#endif
111		#endif
112
113
114		#ifndef FIXED_POINT
115
116	0	#define celt_sqrt(x) ((float)sqrt(x))
117	0	#define celt_rsqrt(x) (1.f/celt_sqrt(x))
118	0	#define celt_rsqrt_norm(x) (celt_rsqrt(x))
119	0	#define celt_cos_norm(x) ((float)cos((.5fPI)(x)))
120	0	#define celt_rcp(x) (1.f/(x))
121	0	#define celt_div(a,b) ((a)/(b))
122	0	#define frac_div32(a,b) ((float)(a)/(b))
123
124		#ifdef FLOAT_APPROX
125
126		/* Note: This assumes radix-2 floating point with the exponent at bits 23..30 and an offset of 127
127		denorm, +/- inf and NaN are not handled */
128
129		/** Base-2 log approximation (log2(x)). */
130		static OPUS_INLINE float celt_log2(float x)
131	0	{
132	0	int integer;
133	0	float frac;
134	0	union {
135	0	float f;
136	0	opus_uint32 i;
137	0	} in;
138	0	in.f = x;
139	0	integer = (in.i>>23)-127;
140	0	in.i -= (opus_uint32)integer<<23;
141	0	frac = in.f - 1.5f;
142	0	frac = -0.41445418f + frac*(0.95909232f
143	0	+ frac(-0.33951290f + frac0.16541097f));
144	0	return 1+integer+frac;
145	0	} Unexecuted instantiation: celt.c:celt_log2 Unexecuted instantiation: pitch_sse.c:celt_log2 Unexecuted instantiation: opus_multistream_encoder.c:celt_log2 Unexecuted instantiation: bands.c:celt_log2 Unexecuted instantiation: celt_encoder.c:celt_log2 Unexecuted instantiation: mathops.c:celt_log2 Unexecuted instantiation: mdct.c:celt_log2 Unexecuted instantiation: modes.c:celt_log2 Unexecuted instantiation: pitch.c:celt_log2 Unexecuted instantiation: celt_lpc.c:celt_log2 Unexecuted instantiation: quant_bands.c:celt_log2 Unexecuted instantiation: vq.c:celt_log2 Unexecuted instantiation: vq_sse2.c:celt_log2 Unexecuted instantiation: opus_decoder.c:celt_log2 Unexecuted instantiation: opus_encoder.c:celt_log2 Unexecuted instantiation: analysis.c:celt_log2 Unexecuted instantiation: celt_decoder.c:celt_log2 Unexecuted instantiation: cwrs.c:celt_log2 Unexecuted instantiation: kiss_fft.c:celt_log2 Unexecuted instantiation: laplace.c:celt_log2
146
147		/** Base-2 exponential approximation (2^x). */
148		static OPUS_INLINE float celt_exp2(float x)
149	0	{
150	0	int integer;
151	0	float frac;
152	0	union {
153	0	float f;
154	0	opus_uint32 i;
155	0	} res;
156	0	integer = (int)floor(x);
157	0	if (integer < -50)
158	0	return 0;
159	0	frac = x-integer;
160		/* K0 = 1, K1 = log(2), K2 = 3-4log(2), K3 = 3log(2) - 2 */
161	0	res.f = 0.99992522f + frac * (0.69583354f
162	0	+ frac * (0.22606716f + 0.078024523f*frac));
163	0	res.i = (res.i + ((opus_uint32)integer<<23)) & 0x7fffffff;
164	0	return res.f;
165	0	} Unexecuted instantiation: celt.c:celt_exp2 Unexecuted instantiation: pitch_sse.c:celt_exp2 Unexecuted instantiation: opus_multistream_encoder.c:celt_exp2 Unexecuted instantiation: bands.c:celt_exp2 Unexecuted instantiation: celt_encoder.c:celt_exp2 Unexecuted instantiation: mathops.c:celt_exp2 Unexecuted instantiation: mdct.c:celt_exp2 Unexecuted instantiation: modes.c:celt_exp2 Unexecuted instantiation: pitch.c:celt_exp2 Unexecuted instantiation: celt_lpc.c:celt_exp2 Unexecuted instantiation: quant_bands.c:celt_exp2 Unexecuted instantiation: vq.c:celt_exp2 Unexecuted instantiation: vq_sse2.c:celt_exp2 Unexecuted instantiation: opus_decoder.c:celt_exp2 Unexecuted instantiation: opus_encoder.c:celt_exp2 Unexecuted instantiation: analysis.c:celt_exp2 Unexecuted instantiation: celt_decoder.c:celt_exp2 Unexecuted instantiation: cwrs.c:celt_exp2 Unexecuted instantiation: kiss_fft.c:celt_exp2 Unexecuted instantiation: laplace.c:celt_exp2
166
167		#else
168		#define celt_log2(x) ((float)(1.442695040888963387*log(x)))
169		#define celt_exp2(x) ((float)exp(0.6931471805599453094*(x)))
170		#endif
171
172		#endif
173
174		#ifdef FIXED_POINT
175
176		#include "os_support.h"
177
178		#ifndef OVERRIDE_CELT_ILOG2
179		/** Integer log in base2. Undefined for zero and negative numbers */
180		static OPUS_INLINE opus_int16 celt_ilog2(opus_int32 x)
181		{
182		celt_sig_assert(x>0);
183		return EC_ILOG(x)-1;
184		}
185		#endif
186
187
188		/** Integer log in base2. Defined for zero, but not for negative numbers */
189		static OPUS_INLINE opus_int16 celt_zlog2(opus_val32 x)
190		{
191		return x <= 0 ? 0 : celt_ilog2(x);
192		}
193
194		opus_val16 celt_rsqrt_norm(opus_val32 x);
195
196		opus_val32 celt_sqrt(opus_val32 x);
197
198		opus_val16 celt_cos_norm(opus_val32 x);
199
200		/** Base-2 logarithm approximation (log2(x)). (Q14 input, Q10 output) */
201		static OPUS_INLINE opus_val16 celt_log2(opus_val32 x)
202		{
203		int i;
204		opus_val16 n, frac;
205		/* -0.41509302963303146, 0.9609890551383969, -0.31836011537636605,
206		0.15530808010959576, -0.08556153059057618 */
207		static const opus_val16 C[5] = {-6801+(1<<(13-DB_SHIFT)), 15746, -5217, 2545, -1401};
208		if (x==0)
209		return -32767;
210		i = celt_ilog2(x);
211		n = VSHR32(x,i-15)-32768-16384;
212		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]))))))));
213		return SHL16(i-13,DB_SHIFT)+SHR16(frac,14-DB_SHIFT);
214		}
215
216		/*
217		K0 = 1
218		K1 = log(2)
219		K2 = 3-4*log(2)
220		K3 = 3*log(2) - 2
221		*/
222		#define D0 16383
223		#define D1 22804
224		#define D2 14819
225		#define D3 10204
226
227		static OPUS_INLINE opus_val32 celt_exp2_frac(opus_val16 x)
228		{
229		opus_val16 frac;
230		frac = SHL16(x, 4);
231		return ADD16(D0, MULT16_16_Q15(frac, ADD16(D1, MULT16_16_Q15(frac, ADD16(D2 , MULT16_16_Q15(D3,frac))))));
232		}
233
234		#undef D0
235		#undef D1
236		#undef D2
237		#undef D3
238
239		/** Base-2 exponential approximation (2^x). (Q10 input, Q16 output) */
240		static OPUS_INLINE opus_val32 celt_exp2(opus_val16 x)
241		{
242		int integer;
243		opus_val16 frac;
244		integer = SHR16(x,10);
245		if (integer>14)
246		return 0x7f000000;
247		else if (integer < -15)
248		return 0;
249		frac = celt_exp2_frac(x-SHL16(integer,10));
250		return VSHR32(EXTEND32(frac), -integer-2);
251		}
252
253		opus_val32 celt_rcp(opus_val32 x);
254
255		#define celt_div(a,b) MULT32_32_Q31((opus_val32)(a),celt_rcp(b))
256
257		opus_val32 frac_div32(opus_val32 a, opus_val32 b);
258
259		#define M1 32767
260		#define M2 -21
261		#define M3 -11943
262		#define M4 4936
263
264		/* Atan approximation using a 4th order polynomial. Input is in Q15 format
265		and normalized by pi/4. Output is in Q15 format */
266		static OPUS_INLINE opus_val16 celt_atan01(opus_val16 x)
267		{
268		return MULT16_16_P15(x, ADD32(M1, MULT16_16_P15(x, ADD32(M2, MULT16_16_P15(x, ADD32(M3, MULT16_16_P15(M4, x)))))));
269		}
270
271		#undef M1
272		#undef M2
273		#undef M3
274		#undef M4
275
276		/* atan2() approximation valid for positive input values */
277		static OPUS_INLINE opus_val16 celt_atan2p(opus_val16 y, opus_val16 x)
278		{
279		if (y < x)
280		{
281		opus_val32 arg;
282		arg = celt_div(SHL32(EXTEND32(y),15),x);
283		if (arg >= 32767)
284		arg = 32767;
285		return SHR16(celt_atan01(EXTRACT16(arg)),1);
286		} else {
287		opus_val32 arg;
288		arg = celt_div(SHL32(EXTEND32(x),15),y);
289		if (arg >= 32767)
290		arg = 32767;
291		return 25736-SHR16(celt_atan01(EXTRACT16(arg)),1);
292		}
293		}
294
295		#endif /* FIXED_POINT */
296		#endif /* MATHOPS_H */