fixed-point: defined celt_maxabs16() as basic operator
[opus.git] / libcelt / mathops.h
1 /* Copyright (C) 2002-2008 Jean-Marc Valin */
2 /**
3    @file mathops.h
4    @brief Various math functions
5 */
6 /*
7    Redistribution and use in source and binary forms, with or without
8    modification, are permitted provided that the following conditions
9    are met:
10    
11    - Redistributions of source code must retain the above copyright
12    notice, this list of conditions and the following disclaimer.
13    
14    - Redistributions in binary form must reproduce the above copyright
15    notice, this list of conditions and the following disclaimer in the
16    documentation and/or other materials provided with the distribution.
17    
18    - Neither the name of the Xiph.org Foundation nor the names of its
19    contributors may be used to endorse or promote products derived from
20    this software without specific prior written permission.
21    
22    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
23    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
24    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
25    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
26    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
27    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
28    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
29    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
30    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
31    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
32    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 */
34
35 #ifndef MATHOPS_H
36 #define MATHOPS_H
37
38 #include "arch.h"
39
40 #ifndef FIXED_POINT
41
42 #define celt_sqrt sqrt
43 #define celt_acos acos
44 #define celt_exp exp
45 #define celt_cos_norm(x) (cos((.5f*M_PI)*(x)))
46 #define celt_atan atan
47 #define celt_rcp(x) (1.f/(x))
48
49 #endif
50
51
52
53 #ifdef FIXED_POINT
54
55 #include "entcode.h"
56 #include "os_support.h"
57
58 #ifndef OVERRIDE_CELT_ILOG2
59 /** Integer log in base2. Undefined for zero and negative numbers */
60 static inline celt_int16_t celt_ilog2(celt_word32_t x)
61 {
62    celt_assert2(x>0, "celt_ilog2() only defined for strictly positive numbers");
63    return EC_ILOG(x)-1;
64 }
65 #endif
66
67 #ifndef OVERRIDE_CELT_MAXABS16
68 static inline celt_word16_t celt_maxabs16(celt_word16_t *x, int len)
69 {
70    int i;
71    celt_word16_t maxval = 0;
72    for (i=0;i<len;i++)
73       maxval = MAX16(maxval, ABS16(x[i]));
74    return maxval;
75 }
76 #endif
77
78 /** Integer log in base2. Defined for zero, but not for negative numbers */
79 static inline celt_int16_t celt_zlog2(celt_word32_t x)
80 {
81    return EC_ILOG(x)-1;
82 }
83
84 /** Sqrt approximation (QX input, QX/2 output) */
85 static inline celt_word32_t celt_sqrt(celt_word32_t x)
86 {
87    int k;
88    celt_word16_t n;
89    celt_word32_t rt;
90    const celt_word16_t C[5] = {23174, 11584, -3011, 1570, -557};
91    if (x==0)
92       return 0;
93    k = (celt_ilog2(x)>>1)-7;
94    x = VSHR32(x, (k<<1));
95    n = x-32768;
96    rt = ADD16(C[0], MULT16_16_Q15(n, ADD16(C[1], MULT16_16_Q15(n, ADD16(C[2], 
97               MULT16_16_Q15(n, ADD16(C[3], MULT16_16_Q15(n, (C[4])))))))));
98    rt = VSHR32(rt,7-k);
99    return rt;
100 }
101
102
103 #define L1 32767
104 #define L2 -7651
105 #define L3 8277
106 #define L4 -626
107
108 static inline celt_word16_t _celt_cos_pi_2(celt_word16_t x)
109 {
110    celt_word16_t x2;
111    
112    x2 = MULT16_16_P15(x,x);
113    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
114                                                                                 ))))))));
115 }
116
117 #undef L1
118 #undef L2
119 #undef L3
120 #undef L4
121
122 static inline celt_word16_t celt_cos_norm(celt_word32_t x)
123 {
124    x = x&0x0001ffff;
125    if (x>SHL32(EXTEND32(1), 16))
126       x = SUB32(SHL32(EXTEND32(1), 17),x);
127    if (x&0x00007fff)
128    {
129       if (x<SHL32(EXTEND32(1), 15))
130       {
131          return _celt_cos_pi_2(EXTRACT16(x));
132       } else {
133          return NEG32(_celt_cos_pi_2(EXTRACT16(65536-x)));
134       }
135    } else {
136       if (x&0x0000ffff)
137          return 0;
138       else if (x&0x0001ffff)
139          return -32767;
140       else
141          return 32767;
142    }
143 }
144
145 static inline celt_word16_t celt_log2(celt_word32_t x)
146 {
147    int i;
148    celt_word16_t n, frac;
149    /*-0.41446   0.96093  -0.33981   0.15600 */
150    const celt_word16_t C[4] = {-6791, 7872, -1392, 319};
151    if (x==0)
152       return -32767;
153    i = celt_ilog2(x);
154    n = VSHR32(x,i-15)-32768-16384;
155    frac = ADD16(C[0], MULT16_16_Q14(n, ADD16(C[1], MULT16_16_Q14(n, ADD16(C[2], MULT16_16_Q14(n, (C[3])))))));
156    /*printf ("%d %d %d %d\n", x, n, ret, SHL16(i-13,8)+SHR16(ret,14-8));*/
157    return SHL16(i-13,8)+SHR16(frac,14-8);
158 }
159
160 /*
161  K0 = 1
162  K1 = log(2)
163  K2 = 3-4*log(2)
164  K3 = 3*log(2) - 2
165 */
166 #define D0 16384
167 #define D1 11356
168 #define D2 3726
169 #define D3 1301
170 /** Base-2 exponential approximation (2^x). (Q11 input, Q16 output) */
171 static inline celt_word32_t celt_exp2(celt_word16_t x)
172 {
173    int integer;
174    celt_word16_t frac;
175    integer = SHR16(x,11);
176    if (integer>14)
177       return 0x7fffffff;
178    else if (integer < -15)
179       return 0;
180    frac = SHL16(x-SHL16(integer,11),3);
181    frac = ADD16(D0, MULT16_16_Q14(frac, ADD16(D1, MULT16_16_Q14(frac, ADD16(D2 , MULT16_16_Q14(D3,frac))))));
182    return VSHR32(EXTEND32(frac), -integer-2);
183 }
184
185 /** Reciprocal approximation (Q15 input, Q16 output) */
186 static inline celt_word32_t celt_rcp(celt_word32_t x)
187 {
188    int i, neg=0;
189    celt_word16_t n, frac;
190    const celt_word16_t C[5] = {21848, -7251, 2403, -934, 327};
191    if (x<0)
192    {
193       neg = 1;
194       x = NEG16(x);
195    }
196    i = celt_ilog2(x);
197    n = VSHR32(x,i-16)-SHL32(EXTEND32(3),15);
198    frac = ADD16(C[0], MULT16_16_Q15(n, ADD16(C[1], MULT16_16_Q15(n, ADD16(C[2], 
199                 MULT16_16_Q15(n, ADD16(C[3], MULT16_16_Q15(n, (C[4])))))))));
200    if (neg)
201       frac = -frac;
202    return VSHR32(EXTEND32(frac),i-16);
203 }
204
205 #endif /* FIXED_POINT */
206
207
208 #endif /* MATHOPS_H */