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