cannam@154: /* Copyright (c) 2002-2008 Jean-Marc Valin
cannam@154:    Copyright (c) 2007-2008 CSIRO
cannam@154:    Copyright (c) 2007-2009 Xiph.Org Foundation
cannam@154:    Written by Jean-Marc Valin */
cannam@154: /**
cannam@154:    @file mathops.h
cannam@154:    @brief Various math functions
cannam@154: */
cannam@154: /*
cannam@154:    Redistribution and use in source and binary forms, with or without
cannam@154:    modification, are permitted provided that the following conditions
cannam@154:    are met:
cannam@154: 
cannam@154:    - Redistributions of source code must retain the above copyright
cannam@154:    notice, this list of conditions and the following disclaimer.
cannam@154: 
cannam@154:    - Redistributions in binary form must reproduce the above copyright
cannam@154:    notice, this list of conditions and the following disclaimer in the
cannam@154:    documentation and/or other materials provided with the distribution.
cannam@154: 
cannam@154:    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
cannam@154:    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
cannam@154:    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
cannam@154:    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
cannam@154:    OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
cannam@154:    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
cannam@154:    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
cannam@154:    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
cannam@154:    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
cannam@154:    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
cannam@154:    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
cannam@154: */
cannam@154: 
cannam@154: #ifndef MATHOPS_H
cannam@154: #define MATHOPS_H
cannam@154: 
cannam@154: #include "arch.h"
cannam@154: #include "entcode.h"
cannam@154: #include "os_support.h"
cannam@154: 
cannam@154: #define PI 3.141592653f
cannam@154: 
cannam@154: /* Multiplies two 16-bit fractional values. Bit-exactness of this macro is important */
cannam@154: #define FRAC_MUL16(a,b) ((16384+((opus_int32)(opus_int16)(a)*(opus_int16)(b)))>>15)
cannam@154: 
cannam@154: unsigned isqrt32(opus_uint32 _val);
cannam@154: 
cannam@154: /* CELT doesn't need it for fixed-point, by analysis.c does. */
cannam@154: #if !defined(FIXED_POINT) || defined(ANALYSIS_C)
cannam@154: #define cA 0.43157974f
cannam@154: #define cB 0.67848403f
cannam@154: #define cC 0.08595542f
cannam@154: #define cE ((float)PI/2)
cannam@154: static OPUS_INLINE float fast_atan2f(float y, float x) {
cannam@154:    float x2, y2;
cannam@154:    x2 = x*x;
cannam@154:    y2 = y*y;
cannam@154:    /* For very small values, we don't care about the answer, so
cannam@154:       we can just return 0. */
cannam@154:    if (x2 + y2 < 1e-18f)
cannam@154:    {
cannam@154:       return 0;
cannam@154:    }
cannam@154:    if(x2<y2){
cannam@154:       float den = (y2 + cB*x2) * (y2 + cC*x2);
cannam@154:       return -x*y*(y2 + cA*x2) / den + (y<0 ? -cE : cE);
cannam@154:    }else{
cannam@154:       float den = (x2 + cB*y2) * (x2 + cC*y2);
cannam@154:       return  x*y*(x2 + cA*y2) / den + (y<0 ? -cE : cE) - (x*y<0 ? -cE : cE);
cannam@154:    }
cannam@154: }
cannam@154: #undef cA
cannam@154: #undef cB
cannam@154: #undef cC
cannam@154: #undef cE
cannam@154: #endif
cannam@154: 
cannam@154: 
cannam@154: #ifndef OVERRIDE_CELT_MAXABS16
cannam@154: static OPUS_INLINE opus_val32 celt_maxabs16(const opus_val16 *x, int len)
cannam@154: {
cannam@154:    int i;
cannam@154:    opus_val16 maxval = 0;
cannam@154:    opus_val16 minval = 0;
cannam@154:    for (i=0;i<len;i++)
cannam@154:    {
cannam@154:       maxval = MAX16(maxval, x[i]);
cannam@154:       minval = MIN16(minval, x[i]);
cannam@154:    }
cannam@154:    return MAX32(EXTEND32(maxval),-EXTEND32(minval));
cannam@154: }
cannam@154: #endif
cannam@154: 
cannam@154: #ifndef OVERRIDE_CELT_MAXABS32
cannam@154: #ifdef FIXED_POINT
cannam@154: static OPUS_INLINE opus_val32 celt_maxabs32(const opus_val32 *x, int len)
cannam@154: {
cannam@154:    int i;
cannam@154:    opus_val32 maxval = 0;
cannam@154:    opus_val32 minval = 0;
cannam@154:    for (i=0;i<len;i++)
cannam@154:    {
cannam@154:       maxval = MAX32(maxval, x[i]);
cannam@154:       minval = MIN32(minval, x[i]);
cannam@154:    }
cannam@154:    return MAX32(maxval, -minval);
cannam@154: }
cannam@154: #else
cannam@154: #define celt_maxabs32(x,len) celt_maxabs16(x,len)
cannam@154: #endif
cannam@154: #endif
cannam@154: 
cannam@154: 
cannam@154: #ifndef FIXED_POINT
cannam@154: 
cannam@154: #define celt_sqrt(x) ((float)sqrt(x))
cannam@154: #define celt_rsqrt(x) (1.f/celt_sqrt(x))
cannam@154: #define celt_rsqrt_norm(x) (celt_rsqrt(x))
cannam@154: #define celt_cos_norm(x) ((float)cos((.5f*PI)*(x)))
cannam@154: #define celt_rcp(x) (1.f/(x))
cannam@154: #define celt_div(a,b) ((a)/(b))
cannam@154: #define frac_div32(a,b) ((float)(a)/(b))
cannam@154: 
cannam@154: #ifdef FLOAT_APPROX
cannam@154: 
cannam@154: /* Note: This assumes radix-2 floating point with the exponent at bits 23..30 and an offset of 127
cannam@154:          denorm, +/- inf and NaN are *not* handled */
cannam@154: 
cannam@154: /** Base-2 log approximation (log2(x)). */
cannam@154: static OPUS_INLINE float celt_log2(float x)
cannam@154: {
cannam@154:    int integer;
cannam@154:    float frac;
cannam@154:    union {
cannam@154:       float f;
cannam@154:       opus_uint32 i;
cannam@154:    } in;
cannam@154:    in.f = x;
cannam@154:    integer = (in.i>>23)-127;
cannam@154:    in.i -= integer<<23;
cannam@154:    frac = in.f - 1.5f;
cannam@154:    frac = -0.41445418f + frac*(0.95909232f
cannam@154:           + frac*(-0.33951290f + frac*0.16541097f));
cannam@154:    return 1+integer+frac;
cannam@154: }
cannam@154: 
cannam@154: /** Base-2 exponential approximation (2^x). */
cannam@154: static OPUS_INLINE float celt_exp2(float x)
cannam@154: {
cannam@154:    int integer;
cannam@154:    float frac;
cannam@154:    union {
cannam@154:       float f;
cannam@154:       opus_uint32 i;
cannam@154:    } res;
cannam@154:    integer = floor(x);
cannam@154:    if (integer < -50)
cannam@154:       return 0;
cannam@154:    frac = x-integer;
cannam@154:    /* K0 = 1, K1 = log(2), K2 = 3-4*log(2), K3 = 3*log(2) - 2 */
cannam@154:    res.f = 0.99992522f + frac * (0.69583354f
cannam@154:            + frac * (0.22606716f + 0.078024523f*frac));
cannam@154:    res.i = (res.i + (integer<<23)) & 0x7fffffff;
cannam@154:    return res.f;
cannam@154: }
cannam@154: 
cannam@154: #else
cannam@154: #define celt_log2(x) ((float)(1.442695040888963387*log(x)))
cannam@154: #define celt_exp2(x) ((float)exp(0.6931471805599453094*(x)))
cannam@154: #endif
cannam@154: 
cannam@154: #endif
cannam@154: 
cannam@154: #ifdef FIXED_POINT
cannam@154: 
cannam@154: #include "os_support.h"
cannam@154: 
cannam@154: #ifndef OVERRIDE_CELT_ILOG2
cannam@154: /** Integer log in base2. Undefined for zero and negative numbers */
cannam@154: static OPUS_INLINE opus_int16 celt_ilog2(opus_int32 x)
cannam@154: {
cannam@154:    celt_sig_assert(x>0);
cannam@154:    return EC_ILOG(x)-1;
cannam@154: }
cannam@154: #endif
cannam@154: 
cannam@154: 
cannam@154: /** Integer log in base2. Defined for zero, but not for negative numbers */
cannam@154: static OPUS_INLINE opus_int16 celt_zlog2(opus_val32 x)
cannam@154: {
cannam@154:    return x <= 0 ? 0 : celt_ilog2(x);
cannam@154: }
cannam@154: 
cannam@154: opus_val16 celt_rsqrt_norm(opus_val32 x);
cannam@154: 
cannam@154: opus_val32 celt_sqrt(opus_val32 x);
cannam@154: 
cannam@154: opus_val16 celt_cos_norm(opus_val32 x);
cannam@154: 
cannam@154: /** Base-2 logarithm approximation (log2(x)). (Q14 input, Q10 output) */
cannam@154: static OPUS_INLINE opus_val16 celt_log2(opus_val32 x)
cannam@154: {
cannam@154:    int i;
cannam@154:    opus_val16 n, frac;
cannam@154:    /* -0.41509302963303146, 0.9609890551383969, -0.31836011537636605,
cannam@154:        0.15530808010959576, -0.08556153059057618 */
cannam@154:    static const opus_val16 C[5] = {-6801+(1<<(13-DB_SHIFT)), 15746, -5217, 2545, -1401};
cannam@154:    if (x==0)
cannam@154:       return -32767;
cannam@154:    i = celt_ilog2(x);
cannam@154:    n = VSHR32(x,i-15)-32768-16384;
cannam@154:    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]))))))));
cannam@154:    return SHL16(i-13,DB_SHIFT)+SHR16(frac,14-DB_SHIFT);
cannam@154: }
cannam@154: 
cannam@154: /*
cannam@154:  K0 = 1
cannam@154:  K1 = log(2)
cannam@154:  K2 = 3-4*log(2)
cannam@154:  K3 = 3*log(2) - 2
cannam@154: */
cannam@154: #define D0 16383
cannam@154: #define D1 22804
cannam@154: #define D2 14819
cannam@154: #define D3 10204
cannam@154: 
cannam@154: static OPUS_INLINE opus_val32 celt_exp2_frac(opus_val16 x)
cannam@154: {
cannam@154:    opus_val16 frac;
cannam@154:    frac = SHL16(x, 4);
cannam@154:    return ADD16(D0, MULT16_16_Q15(frac, ADD16(D1, MULT16_16_Q15(frac, ADD16(D2 , MULT16_16_Q15(D3,frac))))));
cannam@154: }
cannam@154: /** Base-2 exponential approximation (2^x). (Q10 input, Q16 output) */
cannam@154: static OPUS_INLINE opus_val32 celt_exp2(opus_val16 x)
cannam@154: {
cannam@154:    int integer;
cannam@154:    opus_val16 frac;
cannam@154:    integer = SHR16(x,10);
cannam@154:    if (integer>14)
cannam@154:       return 0x7f000000;
cannam@154:    else if (integer < -15)
cannam@154:       return 0;
cannam@154:    frac = celt_exp2_frac(x-SHL16(integer,10));
cannam@154:    return VSHR32(EXTEND32(frac), -integer-2);
cannam@154: }
cannam@154: 
cannam@154: opus_val32 celt_rcp(opus_val32 x);
cannam@154: 
cannam@154: #define celt_div(a,b) MULT32_32_Q31((opus_val32)(a),celt_rcp(b))
cannam@154: 
cannam@154: opus_val32 frac_div32(opus_val32 a, opus_val32 b);
cannam@154: 
cannam@154: #define M1 32767
cannam@154: #define M2 -21
cannam@154: #define M3 -11943
cannam@154: #define M4 4936
cannam@154: 
cannam@154: /* Atan approximation using a 4th order polynomial. Input is in Q15 format
cannam@154:    and normalized by pi/4. Output is in Q15 format */
cannam@154: static OPUS_INLINE opus_val16 celt_atan01(opus_val16 x)
cannam@154: {
cannam@154:    return MULT16_16_P15(x, ADD32(M1, MULT16_16_P15(x, ADD32(M2, MULT16_16_P15(x, ADD32(M3, MULT16_16_P15(M4, x)))))));
cannam@154: }
cannam@154: 
cannam@154: #undef M1
cannam@154: #undef M2
cannam@154: #undef M3
cannam@154: #undef M4
cannam@154: 
cannam@154: /* atan2() approximation valid for positive input values */
cannam@154: static OPUS_INLINE opus_val16 celt_atan2p(opus_val16 y, opus_val16 x)
cannam@154: {
cannam@154:    if (y < x)
cannam@154:    {
cannam@154:       opus_val32 arg;
cannam@154:       arg = celt_div(SHL32(EXTEND32(y),15),x);
cannam@154:       if (arg >= 32767)
cannam@154:          arg = 32767;
cannam@154:       return SHR16(celt_atan01(EXTRACT16(arg)),1);
cannam@154:    } else {
cannam@154:       opus_val32 arg;
cannam@154:       arg = celt_div(SHL32(EXTEND32(x),15),y);
cannam@154:       if (arg >= 32767)
cannam@154:          arg = 32767;
cannam@154:       return 25736-SHR16(celt_atan01(EXTRACT16(arg)),1);
cannam@154:    }
cannam@154: }
cannam@154: 
cannam@154: #endif /* FIXED_POINT */
cannam@154: #endif /* MATHOPS_H */