Line data Source code
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 : #define PI 3.141592653f
42 :
43 : /* Multiplies two 16-bit fractional values. Bit-exactness of this macro is important */
44 : #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 : #define cA 0.43157974f
51 : #define cB 0.67848403f
52 : #define cC 0.08595542f
53 : #define cE ((float)PI/2)
54 0 : static OPUS_INLINE float fast_atan2f(float y, float x) {
55 : 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 : {
62 0 : return 0;
63 : }
64 0 : if(x2<y2){
65 0 : float den = (y2 + cB*x2) * (y2 + cC*x2);
66 0 : return -x*y*(y2 + cA*x2) / den + (y<0 ? -cE : cE);
67 : }else{
68 0 : float den = (x2 + cB*y2) * (x2 + cC*y2);
69 0 : return x*y*(x2 + cA*y2) / den + (y<0 ? -cE : cE) - (x*y<0 ? -cE : cE);
70 : }
71 : }
72 : #undef cA
73 : #undef cB
74 : #undef cC
75 : #undef cD
76 : #endif
77 :
78 :
79 : #ifndef OVERRIDE_CELT_MAXABS16
80 0 : static OPUS_INLINE opus_val32 celt_maxabs16(const opus_val16 *x, int len)
81 : {
82 : int i;
83 0 : opus_val16 maxval = 0;
84 0 : opus_val16 minval = 0;
85 0 : for (i=0;i<len;i++)
86 : {
87 0 : maxval = MAX16(maxval, x[i]);
88 0 : minval = MIN16(minval, x[i]);
89 : }
90 0 : return MAX32(EXTEND32(maxval),-EXTEND32(minval));
91 : }
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 : #define celt_sqrt(x) ((float)sqrt(x))
117 : #define celt_rsqrt(x) (1.f/celt_sqrt(x))
118 : #define celt_rsqrt_norm(x) (celt_rsqrt(x))
119 : #define celt_cos_norm(x) ((float)cos((.5f*PI)*(x)))
120 : #define celt_rcp(x) (1.f/(x))
121 : #define celt_div(a,b) ((a)/(b))
122 : #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 : {
132 : int integer;
133 : float frac;
134 : union {
135 : float f;
136 : opus_uint32 i;
137 : } in;
138 : in.f = x;
139 : integer = (in.i>>23)-127;
140 : in.i -= integer<<23;
141 : frac = in.f - 1.5f;
142 : frac = -0.41445418f + frac*(0.95909232f
143 : + frac*(-0.33951290f + frac*0.16541097f));
144 : return 1+integer+frac;
145 : }
146 :
147 : /** Base-2 exponential approximation (2^x). */
148 : static OPUS_INLINE float celt_exp2(float x)
149 : {
150 : int integer;
151 : float frac;
152 : union {
153 : float f;
154 : opus_uint32 i;
155 : } res;
156 : integer = floor(x);
157 : if (integer < -50)
158 : return 0;
159 : frac = x-integer;
160 : /* K0 = 1, K1 = log(2), K2 = 3-4*log(2), K3 = 3*log(2) - 2 */
161 : res.f = 0.99992522f + frac * (0.69583354f
162 : + frac * (0.22606716f + 0.078024523f*frac));
163 : res.i = (res.i + (integer<<23)) & 0x7fffffff;
164 : return res.f;
165 : }
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_assert2(x>0, "celt_ilog2() only defined for strictly positive numbers");
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 : /** Base-2 exponential approximation (2^x). (Q10 input, Q16 output) */
234 : static OPUS_INLINE opus_val32 celt_exp2(opus_val16 x)
235 : {
236 : int integer;
237 : opus_val16 frac;
238 : integer = SHR16(x,10);
239 : if (integer>14)
240 : return 0x7f000000;
241 : else if (integer < -15)
242 : return 0;
243 : frac = celt_exp2_frac(x-SHL16(integer,10));
244 : return VSHR32(EXTEND32(frac), -integer-2);
245 : }
246 :
247 : opus_val32 celt_rcp(opus_val32 x);
248 :
249 : #define celt_div(a,b) MULT32_32_Q31((opus_val32)(a),celt_rcp(b))
250 :
251 : opus_val32 frac_div32(opus_val32 a, opus_val32 b);
252 :
253 : #define M1 32767
254 : #define M2 -21
255 : #define M3 -11943
256 : #define M4 4936
257 :
258 : /* Atan approximation using a 4th order polynomial. Input is in Q15 format
259 : and normalized by pi/4. Output is in Q15 format */
260 : static OPUS_INLINE opus_val16 celt_atan01(opus_val16 x)
261 : {
262 : return MULT16_16_P15(x, ADD32(M1, MULT16_16_P15(x, ADD32(M2, MULT16_16_P15(x, ADD32(M3, MULT16_16_P15(M4, x)))))));
263 : }
264 :
265 : #undef M1
266 : #undef M2
267 : #undef M3
268 : #undef M4
269 :
270 : /* atan2() approximation valid for positive input values */
271 : static OPUS_INLINE opus_val16 celt_atan2p(opus_val16 y, opus_val16 x)
272 : {
273 : if (y < x)
274 : {
275 : opus_val32 arg;
276 : arg = celt_div(SHL32(EXTEND32(y),15),x);
277 : if (arg >= 32767)
278 : arg = 32767;
279 : return SHR16(celt_atan01(EXTRACT16(arg)),1);
280 : } else {
281 : opus_val32 arg;
282 : arg = celt_div(SHL32(EXTEND32(x),15),y);
283 : if (arg >= 32767)
284 : arg = 32767;
285 : return 25736-SHR16(celt_atan01(EXTRACT16(arg)),1);
286 : }
287 : }
288 :
289 : #endif /* FIXED_POINT */
290 : #endif /* MATHOPS_H */
|
