Alexandre Lision | 744f742 | 2013-09-25 11:39:37 -0400 | [diff] [blame] | 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 | /* 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) |
| 43 | |
| 44 | unsigned isqrt32(opus_uint32 _val); |
| 45 | |
| 46 | #ifndef FIXED_POINT |
| 47 | |
| 48 | #define PI 3.141592653f |
| 49 | #define celt_sqrt(x) ((float)sqrt(x)) |
| 50 | #define celt_rsqrt(x) (1.f/celt_sqrt(x)) |
| 51 | #define celt_rsqrt_norm(x) (celt_rsqrt(x)) |
| 52 | #define celt_cos_norm(x) ((float)cos((.5f*PI)*(x))) |
| 53 | #define celt_rcp(x) (1.f/(x)) |
| 54 | #define celt_div(a,b) ((a)/(b)) |
| 55 | #define frac_div32(a,b) ((float)(a)/(b)) |
| 56 | |
| 57 | #ifdef FLOAT_APPROX |
| 58 | |
| 59 | /* Note: This assumes radix-2 floating point with the exponent at bits 23..30 and an offset of 127 |
| 60 | denorm, +/- inf and NaN are *not* handled */ |
| 61 | |
| 62 | /** Base-2 log approximation (log2(x)). */ |
| 63 | static inline float celt_log2(float x) |
| 64 | { |
| 65 | int integer; |
| 66 | float frac; |
| 67 | union { |
| 68 | float f; |
| 69 | opus_uint32 i; |
| 70 | } in; |
| 71 | in.f = x; |
| 72 | integer = (in.i>>23)-127; |
| 73 | in.i -= integer<<23; |
| 74 | frac = in.f - 1.5f; |
| 75 | frac = -0.41445418f + frac*(0.95909232f |
| 76 | + frac*(-0.33951290f + frac*0.16541097f)); |
| 77 | return 1+integer+frac; |
| 78 | } |
| 79 | |
| 80 | /** Base-2 exponential approximation (2^x). */ |
| 81 | static inline float celt_exp2(float x) |
| 82 | { |
| 83 | int integer; |
| 84 | float frac; |
| 85 | union { |
| 86 | float f; |
| 87 | opus_uint32 i; |
| 88 | } res; |
| 89 | integer = floor(x); |
| 90 | if (integer < -50) |
| 91 | return 0; |
| 92 | frac = x-integer; |
| 93 | /* K0 = 1, K1 = log(2), K2 = 3-4*log(2), K3 = 3*log(2) - 2 */ |
| 94 | res.f = 0.99992522f + frac * (0.69583354f |
| 95 | + frac * (0.22606716f + 0.078024523f*frac)); |
| 96 | res.i = (res.i + (integer<<23)) & 0x7fffffff; |
| 97 | return res.f; |
| 98 | } |
| 99 | |
| 100 | #else |
| 101 | #define celt_log2(x) ((float)(1.442695040888963387*log(x))) |
| 102 | #define celt_exp2(x) ((float)exp(0.6931471805599453094*(x))) |
| 103 | #endif |
| 104 | |
| 105 | #endif |
| 106 | |
| 107 | #ifdef FIXED_POINT |
| 108 | |
| 109 | #include "os_support.h" |
| 110 | |
| 111 | #ifndef OVERRIDE_CELT_ILOG2 |
| 112 | /** Integer log in base2. Undefined for zero and negative numbers */ |
| 113 | static inline opus_int16 celt_ilog2(opus_int32 x) |
| 114 | { |
| 115 | celt_assert2(x>0, "celt_ilog2() only defined for strictly positive numbers"); |
| 116 | return EC_ILOG(x)-1; |
| 117 | } |
| 118 | #endif |
| 119 | |
| 120 | #ifndef OVERRIDE_CELT_MAXABS16 |
| 121 | static inline opus_val32 celt_maxabs16(const opus_val16 *x, int len) |
| 122 | { |
| 123 | int i; |
| 124 | opus_val16 maxval = 0; |
| 125 | opus_val16 minval = 0; |
| 126 | for (i=0;i<len;i++) |
| 127 | { |
| 128 | maxval = MAX16(maxval, x[i]); |
| 129 | minval = MIN16(minval, x[i]); |
| 130 | } |
| 131 | return MAX32(EXTEND32(maxval),-EXTEND32(minval)); |
| 132 | } |
| 133 | #endif |
| 134 | |
| 135 | #ifndef OVERRIDE_CELT_MAXABS32 |
| 136 | static inline opus_val32 celt_maxabs32(opus_val32 *x, int len) |
| 137 | { |
| 138 | int i; |
| 139 | opus_val32 maxval = 0; |
| 140 | for (i=0;i<len;i++) |
| 141 | maxval = MAX32(maxval, ABS32(x[i])); |
| 142 | return maxval; |
| 143 | } |
| 144 | #endif |
| 145 | |
| 146 | /** Integer log in base2. Defined for zero, but not for negative numbers */ |
| 147 | static inline opus_int16 celt_zlog2(opus_val32 x) |
| 148 | { |
| 149 | return x <= 0 ? 0 : celt_ilog2(x); |
| 150 | } |
| 151 | |
| 152 | opus_val16 celt_rsqrt_norm(opus_val32 x); |
| 153 | |
| 154 | opus_val32 celt_sqrt(opus_val32 x); |
| 155 | |
| 156 | opus_val16 celt_cos_norm(opus_val32 x); |
| 157 | |
| 158 | static inline opus_val16 celt_log2(opus_val32 x) |
| 159 | { |
| 160 | int i; |
| 161 | opus_val16 n, frac; |
| 162 | /* -0.41509302963303146, 0.9609890551383969, -0.31836011537636605, |
| 163 | 0.15530808010959576, -0.08556153059057618 */ |
| 164 | static const opus_val16 C[5] = {-6801+(1<<(13-DB_SHIFT)), 15746, -5217, 2545, -1401}; |
| 165 | if (x==0) |
| 166 | return -32767; |
| 167 | i = celt_ilog2(x); |
| 168 | n = VSHR32(x,i-15)-32768-16384; |
| 169 | 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])))))))); |
| 170 | return SHL16(i-13,DB_SHIFT)+SHR16(frac,14-DB_SHIFT); |
| 171 | } |
| 172 | |
| 173 | /* |
| 174 | K0 = 1 |
| 175 | K1 = log(2) |
| 176 | K2 = 3-4*log(2) |
| 177 | K3 = 3*log(2) - 2 |
| 178 | */ |
| 179 | #define D0 16383 |
| 180 | #define D1 22804 |
| 181 | #define D2 14819 |
| 182 | #define D3 10204 |
| 183 | /** Base-2 exponential approximation (2^x). (Q10 input, Q16 output) */ |
| 184 | static inline opus_val32 celt_exp2(opus_val16 x) |
| 185 | { |
| 186 | int integer; |
| 187 | opus_val16 frac; |
| 188 | integer = SHR16(x,10); |
| 189 | if (integer>14) |
| 190 | return 0x7f000000; |
| 191 | else if (integer < -15) |
| 192 | return 0; |
| 193 | frac = SHL16(x-SHL16(integer,10),4); |
| 194 | frac = ADD16(D0, MULT16_16_Q15(frac, ADD16(D1, MULT16_16_Q15(frac, ADD16(D2 , MULT16_16_Q15(D3,frac)))))); |
| 195 | return VSHR32(EXTEND32(frac), -integer-2); |
| 196 | } |
| 197 | |
| 198 | opus_val32 celt_rcp(opus_val32 x); |
| 199 | |
| 200 | #define celt_div(a,b) MULT32_32_Q31((opus_val32)(a),celt_rcp(b)) |
| 201 | |
| 202 | opus_val32 frac_div32(opus_val32 a, opus_val32 b); |
| 203 | |
| 204 | #define M1 32767 |
| 205 | #define M2 -21 |
| 206 | #define M3 -11943 |
| 207 | #define M4 4936 |
| 208 | |
| 209 | /* Atan approximation using a 4th order polynomial. Input is in Q15 format |
| 210 | and normalized by pi/4. Output is in Q15 format */ |
| 211 | static inline opus_val16 celt_atan01(opus_val16 x) |
| 212 | { |
| 213 | return MULT16_16_P15(x, ADD32(M1, MULT16_16_P15(x, ADD32(M2, MULT16_16_P15(x, ADD32(M3, MULT16_16_P15(M4, x))))))); |
| 214 | } |
| 215 | |
| 216 | #undef M1 |
| 217 | #undef M2 |
| 218 | #undef M3 |
| 219 | #undef M4 |
| 220 | |
| 221 | /* atan2() approximation valid for positive input values */ |
| 222 | static inline opus_val16 celt_atan2p(opus_val16 y, opus_val16 x) |
| 223 | { |
| 224 | if (y < x) |
| 225 | { |
| 226 | opus_val32 arg; |
| 227 | arg = celt_div(SHL32(EXTEND32(y),15),x); |
| 228 | if (arg >= 32767) |
| 229 | arg = 32767; |
| 230 | return SHR16(celt_atan01(EXTRACT16(arg)),1); |
| 231 | } else { |
| 232 | opus_val32 arg; |
| 233 | arg = celt_div(SHL32(EXTEND32(x),15),y); |
| 234 | if (arg >= 32767) |
| 235 | arg = 32767; |
| 236 | return 25736-SHR16(celt_atan01(EXTRACT16(arg)),1); |
| 237 | } |
| 238 | } |
| 239 | |
| 240 | #endif /* FIXED_POINT */ |
| 241 | #endif /* MATHOPS_H */ |