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 */
7 @brief Various math functions
10 Redistribution and use in source and binary forms, with or without
11 modification, are permitted provided that the following conditions
14 - Redistributions of source code must retain the above copyright
15 notice, this list of conditions and the following disclaimer.
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.
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.
39 #include "os_support.h"
41 /* Multiplies two 16-bit fractional values. Bit-exactness of this macro is important */
42 #define FRAC_MUL16(a,b) ((16384+((opus_int32)(opus_int16)(a)*(opus_int16)(b)))>>15)
44 unsigned isqrt32(opus_uint32 _val
);
46 #ifndef OVERRIDE_CELT_MAXABS16
47 static OPUS_INLINE opus_val32
celt_maxabs16(const opus_val16
*x
, int len
)
50 opus_val16 maxval
= 0;
51 opus_val16 minval
= 0;
54 maxval
= MAX16(maxval
, x
[i
]);
55 minval
= MIN16(minval
, x
[i
]);
57 return MAX32(EXTEND32(maxval
),-EXTEND32(minval
));
61 #ifndef OVERRIDE_CELT_MAXABS32
63 static OPUS_INLINE opus_val32
celt_maxabs32(const opus_val32
*x
, int len
)
66 opus_val32 maxval
= 0;
67 opus_val32 minval
= 0;
70 maxval
= MAX32(maxval
, x
[i
]);
71 minval
= MIN32(minval
, x
[i
]);
73 return MAX32(maxval
, -minval
);
76 #define celt_maxabs32(x,len) celt_maxabs16(x,len)
83 #define PI 3.141592653f
84 #define celt_sqrt(x) ((float)sqrt(x))
85 #define celt_rsqrt(x) (1.f/celt_sqrt(x))
86 #define celt_rsqrt_norm(x) (celt_rsqrt(x))
87 #define celt_cos_norm(x) ((float)cos((.5f*PI)*(x)))
88 #define celt_rcp(x) (1.f/(x))
89 #define celt_div(a,b) ((a)/(b))
90 #define frac_div32(a,b) ((float)(a)/(b))
94 /* Note: This assumes radix-2 floating point with the exponent at bits 23..30 and an offset of 127
95 denorm, +/- inf and NaN are *not* handled */
97 /** Base-2 log approximation (log2(x)). */
98 static OPUS_INLINE
float celt_log2(float x
)
107 integer
= (in
.i
>>23)-127;
110 frac
= -0.41445418f
+ frac
*(0.95909232f
111 + frac
*(-0.33951290f
+ frac
*0.16541097f
));
112 return 1+integer
+frac
;
115 /** Base-2 exponential approximation (2^x). */
116 static OPUS_INLINE
float celt_exp2(float x
)
128 /* K0 = 1, K1 = log(2), K2 = 3-4*log(2), K3 = 3*log(2) - 2 */
129 res
.f
= 0.99992522f
+ frac
* (0.69583354f
130 + frac
* (0.22606716f
+ 0.078024523f
*frac
));
131 res
.i
= (res
.i
+ (integer
<<23)) & 0x7fffffff;
136 #define celt_log2(x) ((float)(1.442695040888963387*log(x)))
137 #define celt_exp2(x) ((float)exp(0.6931471805599453094*(x)))
144 #include "os_support.h"
146 #ifndef OVERRIDE_CELT_ILOG2
147 /** Integer log in base2. Undefined for zero and negative numbers */
148 static OPUS_INLINE opus_int16
celt_ilog2(opus_int32 x
)
150 celt_assert2(x
>0, "celt_ilog2() only defined for strictly positive numbers");
156 /** Integer log in base2. Defined for zero, but not for negative numbers */
157 static OPUS_INLINE opus_int16
celt_zlog2(opus_val32 x
)
159 return x
<= 0 ? 0 : celt_ilog2(x
);
162 opus_val16
celt_rsqrt_norm(opus_val32 x
);
164 opus_val32
celt_sqrt(opus_val32 x
);
166 opus_val16
celt_cos_norm(opus_val32 x
);
168 /** Base-2 logarithm approximation (log2(x)). (Q14 input, Q10 output) */
169 static OPUS_INLINE opus_val16
celt_log2(opus_val32 x
)
173 /* -0.41509302963303146, 0.9609890551383969, -0.31836011537636605,
174 0.15530808010959576, -0.08556153059057618 */
175 static const opus_val16 C
[5] = {-6801+(1<<(13-DB_SHIFT
)), 15746, -5217, 2545, -1401};
179 n
= VSHR32(x
,i
-15)-32768-16384;
180 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]))))))));
181 return SHL16(i
-13,DB_SHIFT
)+SHR16(frac
,14-DB_SHIFT
);
195 static OPUS_INLINE opus_val32
celt_exp2_frac(opus_val16 x
)
199 return ADD16(D0
, MULT16_16_Q15(frac
, ADD16(D1
, MULT16_16_Q15(frac
, ADD16(D2
, MULT16_16_Q15(D3
,frac
))))));
201 /** Base-2 exponential approximation (2^x). (Q10 input, Q16 output) */
202 static OPUS_INLINE opus_val32
celt_exp2(opus_val16 x
)
206 integer
= SHR16(x
,10);
209 else if (integer
< -15)
211 frac
= celt_exp2_frac(x
-SHL16(integer
,10));
212 return VSHR32(EXTEND32(frac
), -integer
-2);
215 opus_val32
celt_rcp(opus_val32 x
);
217 #define celt_div(a,b) MULT32_32_Q31((opus_val32)(a),celt_rcp(b))
219 opus_val32
frac_div32(opus_val32 a
, opus_val32 b
);
226 /* Atan approximation using a 4th order polynomial. Input is in Q15 format
227 and normalized by pi/4. Output is in Q15 format */
228 static OPUS_INLINE opus_val16
celt_atan01(opus_val16 x
)
230 return MULT16_16_P15(x
, ADD32(M1
, MULT16_16_P15(x
, ADD32(M2
, MULT16_16_P15(x
, ADD32(M3
, MULT16_16_P15(M4
, x
)))))));
238 /* atan2() approximation valid for positive input values */
239 static OPUS_INLINE opus_val16
celt_atan2p(opus_val16 y
, opus_val16 x
)
244 arg
= celt_div(SHL32(EXTEND32(y
),15),x
);
247 return SHR16(celt_atan01(EXTRACT16(arg
)),1);
250 arg
= celt_div(SHL32(EXTEND32(x
),15),y
);
253 return 25736-SHR16(celt_atan01(EXTRACT16(arg
)),1);
257 #endif /* FIXED_POINT */
258 #endif /* MATHOPS_H */