Home | History | Annotate | Download | only in src
      1 
      2 /* -----------------------------------------------------------------------------------------------------------
      3 Software License for The Fraunhofer FDK AAC Codec Library for Android
      4 
      5  Copyright  1995 - 2013 Fraunhofer-Gesellschaft zur Frderung der angewandten Forschung e.V.
      6   All rights reserved.
      7 
      8  1.    INTRODUCTION
      9 The Fraunhofer FDK AAC Codec Library for Android ("FDK AAC Codec") is software that implements
     10 the MPEG Advanced Audio Coding ("AAC") encoding and decoding scheme for digital audio.
     11 This FDK AAC Codec software is intended to be used on a wide variety of Android devices.
     12 
     13 AAC's HE-AAC and HE-AAC v2 versions are regarded as today's most efficient general perceptual
     14 audio codecs. AAC-ELD is considered the best-performing full-bandwidth communications codec by
     15 independent studies and is widely deployed. AAC has been standardized by ISO and IEC as part
     16 of the MPEG specifications.
     17 
     18 Patent licenses for necessary patent claims for the FDK AAC Codec (including those of Fraunhofer)
     19 may be obtained through Via Licensing (www.vialicensing.com) or through the respective patent owners
     20 individually for the purpose of encoding or decoding bit streams in products that are compliant with
     21 the ISO/IEC MPEG audio standards. Please note that most manufacturers of Android devices already license
     22 these patent claims through Via Licensing or directly from the patent owners, and therefore FDK AAC Codec
     23 software may already be covered under those patent licenses when it is used for those licensed purposes only.
     24 
     25 Commercially-licensed AAC software libraries, including floating-point versions with enhanced sound quality,
     26 are also available from Fraunhofer. Users are encouraged to check the Fraunhofer website for additional
     27 applications information and documentation.
     28 
     29 2.    COPYRIGHT LICENSE
     30 
     31 Redistribution and use in source and binary forms, with or without modification, are permitted without
     32 payment of copyright license fees provided that you satisfy the following conditions:
     33 
     34 You must retain the complete text of this software license in redistributions of the FDK AAC Codec or
     35 your modifications thereto in source code form.
     36 
     37 You must retain the complete text of this software license in the documentation and/or other materials
     38 provided with redistributions of the FDK AAC Codec or your modifications thereto in binary form.
     39 You must make available free of charge copies of the complete source code of the FDK AAC Codec and your
     40 modifications thereto to recipients of copies in binary form.
     41 
     42 The name of Fraunhofer may not be used to endorse or promote products derived from this library without
     43 prior written permission.
     44 
     45 You may not charge copyright license fees for anyone to use, copy or distribute the FDK AAC Codec
     46 software or your modifications thereto.
     47 
     48 Your modified versions of the FDK AAC Codec must carry prominent notices stating that you changed the software
     49 and the date of any change. For modified versions of the FDK AAC Codec, the term
     50 "Fraunhofer FDK AAC Codec Library for Android" must be replaced by the term
     51 "Third-Party Modified Version of the Fraunhofer FDK AAC Codec Library for Android."
     52 
     53 3.    NO PATENT LICENSE
     54 
     55 NO EXPRESS OR IMPLIED LICENSES TO ANY PATENT CLAIMS, including without limitation the patents of Fraunhofer,
     56 ARE GRANTED BY THIS SOFTWARE LICENSE. Fraunhofer provides no warranty of patent non-infringement with
     57 respect to this software.
     58 
     59 You may use this FDK AAC Codec software or modifications thereto only for purposes that are authorized
     60 by appropriate patent licenses.
     61 
     62 4.    DISCLAIMER
     63 
     64 This FDK AAC Codec software is provided by Fraunhofer on behalf of the copyright holders and contributors
     65 "AS IS" and WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES, including but not limited to the implied warranties
     66 of merchantability and fitness for a particular purpose. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
     67 CONTRIBUTORS BE LIABLE for any direct, indirect, incidental, special, exemplary, or consequential damages,
     68 including but not limited to procurement of substitute goods or services; loss of use, data, or profits,
     69 or business interruption, however caused and on any theory of liability, whether in contract, strict
     70 liability, or tort (including negligence), arising in any way out of the use of this software, even if
     71 advised of the possibility of such damage.
     72 
     73 5.    CONTACT INFORMATION
     74 
     75 Fraunhofer Institute for Integrated Circuits IIS
     76 Attention: Audio and Multimedia Departments - FDK AAC LL
     77 Am Wolfsmantel 33
     78 91058 Erlangen, Germany
     79 
     80 www.iis.fraunhofer.de/amm
     81 amm-info (at) iis.fraunhofer.de
     82 ----------------------------------------------------------------------------------------------------------- */
     83 
     84 /***************************  Fraunhofer IIS FDK Tools  **********************
     85 
     86    Author(s):   M. Gayer
     87    Description: Fixed point specific mathematical functions
     88 
     89 ******************************************************************************/
     90 
     91 #include "fixpoint_math.h"
     92 
     93 
     94 #define MAX_LD_PRECISION 10
     95 #define LD_PRECISION     10
     96 
     97 /* Taylor series coeffcients for ln(1-x), centered at 0 (MacLaurin polinomial). */
     98 #ifndef LDCOEFF_16BIT
     99 LNK_SECTION_CONSTDATA_L1
    100 static const FIXP_DBL ldCoeff[MAX_LD_PRECISION] = {
    101     FL2FXCONST_DBL(-1.0),
    102     FL2FXCONST_DBL(-1.0/2.0),
    103     FL2FXCONST_DBL(-1.0/3.0),
    104     FL2FXCONST_DBL(-1.0/4.0),
    105     FL2FXCONST_DBL(-1.0/5.0),
    106     FL2FXCONST_DBL(-1.0/6.0),
    107     FL2FXCONST_DBL(-1.0/7.0),
    108     FL2FXCONST_DBL(-1.0/8.0),
    109     FL2FXCONST_DBL(-1.0/9.0),
    110     FL2FXCONST_DBL(-1.0/10.0)
    111 };
    112 #else
    113 LNK_SECTION_CONSTDATA_L1
    114 static const FIXP_SGL ldCoeff[MAX_LD_PRECISION] = {
    115     FL2FXCONST_SGL(-1.0),
    116     FL2FXCONST_SGL(-1.0/2.0),
    117     FL2FXCONST_SGL(-1.0/3.0),
    118     FL2FXCONST_SGL(-1.0/4.0),
    119     FL2FXCONST_SGL(-1.0/5.0),
    120     FL2FXCONST_SGL(-1.0/6.0),
    121     FL2FXCONST_SGL(-1.0/7.0),
    122     FL2FXCONST_SGL(-1.0/8.0),
    123     FL2FXCONST_SGL(-1.0/9.0),
    124     FL2FXCONST_SGL(-1.0/10.0)
    125 };
    126 #endif
    127 
    128 /*****************************************************************************
    129 
    130   functionname: CalcLdData
    131   description:  Delivers the Logarithm Dualis ld(op)/LD_DATA_SCALING with polynomial approximation.
    132   input:        Input op is assumed to be double precision fractional 0 < op < 1.0
    133                 This function does not accept negative values.
    134   output:       For op == 0, the result is saturated to -1.0
    135                 This function does not return positive values since input values are treated as fractional values.
    136                 It does not make sense to input an integer value into this function (and expect a positive output value)
    137                 since input values are treated as fractional values.
    138 
    139 *****************************************************************************/
    140 
    141 LNK_SECTION_CODE_L1
    142 FIXP_DBL CalcLdData(FIXP_DBL op)
    143 {
    144   return fLog2(op, 0);
    145 }
    146 
    147 
    148 /*****************************************************************************
    149   functionname: LdDataVector
    150 *****************************************************************************/
    151 LNK_SECTION_CODE_L1
    152 void LdDataVector(  FIXP_DBL    *srcVector,
    153                     FIXP_DBL    *destVector,
    154                     INT         n)
    155 {
    156     INT i;
    157     for ( i=0; i<n; i++) {
    158         destVector[i] = CalcLdData(srcVector[i]);
    159     }
    160 }
    161 
    162 
    163 
    164 #define MAX_POW2_PRECISION 8
    165 #ifndef SINETABLE_16BIT
    166 	#define POW2_PRECISION MAX_POW2_PRECISION
    167 #else
    168 	#define POW2_PRECISION 5
    169 #endif
    170 
    171 /*
    172   Taylor series coefficients of the function x^2. The first coefficient is
    173   ommited (equal to 1.0).
    174 
    175   pow2Coeff[i-1] = (1/i!) d^i(2^x)/dx^i, i=1..MAX_POW2_PRECISION
    176   To evaluate the taylor series around x = 0, the coefficients are: 1/!i * ln(2)^i
    177  */
    178 #ifndef POW2COEFF_16BIT
    179 LNK_SECTION_CONSTDATA_L1
    180 static const FIXP_DBL pow2Coeff[MAX_POW2_PRECISION] = {
    181     FL2FXCONST_DBL(0.693147180559945309417232121458177),    /* ln(2)^1 /1! */
    182     FL2FXCONST_DBL(0.240226506959100712333551263163332),    /* ln(2)^2 /2! */
    183     FL2FXCONST_DBL(0.0555041086648215799531422637686218),   /* ln(2)^3 /3! */
    184     FL2FXCONST_DBL(0.00961812910762847716197907157365887),  /* ln(2)^4 /4! */
    185     FL2FXCONST_DBL(0.00133335581464284434234122219879962),  /* ln(2)^5 /5! */
    186     FL2FXCONST_DBL(1.54035303933816099544370973327423e-4),  /* ln(2)^6 /6! */
    187     FL2FXCONST_DBL(1.52527338040598402800254390120096e-5),  /* ln(2)^7 /7! */
    188     FL2FXCONST_DBL(1.32154867901443094884037582282884e-6)   /* ln(2)^8 /8! */
    189 };
    190 #else
    191 LNK_SECTION_CONSTDATA_L1
    192 static const FIXP_SGL pow2Coeff[MAX_POW2_PRECISION] = {
    193     FL2FXCONST_SGL(0.693147180559945309417232121458177),    /* ln(2)^1 /1! */
    194     FL2FXCONST_SGL(0.240226506959100712333551263163332),    /* ln(2)^2 /2! */
    195     FL2FXCONST_SGL(0.0555041086648215799531422637686218),   /* ln(2)^3 /3! */
    196     FL2FXCONST_SGL(0.00961812910762847716197907157365887),  /* ln(2)^4 /4! */
    197     FL2FXCONST_SGL(0.00133335581464284434234122219879962),  /* ln(2)^5 /5! */
    198     FL2FXCONST_SGL(1.54035303933816099544370973327423e-4),  /* ln(2)^6 /6! */
    199     FL2FXCONST_SGL(1.52527338040598402800254390120096e-5),  /* ln(2)^7 /7! */
    200     FL2FXCONST_SGL(1.32154867901443094884037582282884e-6)   /* ln(2)^8 /8! */
    201 };
    202 #endif
    203 
    204 
    205 
    206 /*****************************************************************************
    207 
    208     functionname: mul_dbl_sgl_rnd
    209     description:  Multiply with round.
    210 *****************************************************************************/
    211 
    212 /* for rounding a dfract to fract */
    213 #define ACCU_R (LONG) 0x00008000
    214 
    215 LNK_SECTION_CODE_L1
    216 FIXP_DBL mul_dbl_sgl_rnd (const FIXP_DBL op1, const FIXP_SGL op2)
    217 {
    218   FIXP_DBL prod;
    219   LONG v = (LONG)(op1);
    220   SHORT u = (SHORT)(op2);
    221 
    222   LONG low = u*(v&SGL_MASK);
    223   low = (low+(ACCU_R>>1)) >> (FRACT_BITS-1);  /* round */
    224   LONG high = u * ((v>>FRACT_BITS)<<1);
    225 
    226   prod = (LONG)(high+low);
    227 
    228   return((FIXP_DBL)prod);
    229 }
    230 
    231 
    232 /*****************************************************************************
    233 
    234   functionname: CalcInvLdData
    235   description:  Delivers the inverse of function CalcLdData().
    236                 Delivers 2^(op*LD_DATA_SCALING)
    237   input:        Input op is assumed to be fractional -1.0 < op < 1.0
    238   output:       For op == 0, the result is MAXVAL_DBL (almost 1.0).
    239                 For negative input values the output should be treated as a positive fractional value.
    240                 For positive input values the output should be treated as a positive integer value.
    241                 This function does not output negative values.
    242 
    243 *****************************************************************************/
    244 LNK_SECTION_CODE_L1
    245 /* This table is used for lookup 2^x with   */
    246 /* x in range [0...1.0[ in steps of 1/32 */
    247 LNK_SECTION_DATA_L1 static const UINT exp2_tab_long[32]={
    248 0x40000000,0x4166C34C,0x42D561B4,0x444C0740,
    249 0x45CAE0F2,0x47521CC6,0x48E1E9BA,0x4A7A77D4,
    250 0x4C1BF829,0x4DC69CDD,0x4F7A9930,0x51382182,
    251 0x52FF6B55,0x54D0AD5A,0x56AC1F75,0x5891FAC1,
    252 0x5A82799A,0x5C7DD7A4,0x5E8451D0,0x60962665,
    253 0x62B39509,0x64DCDEC3,0x6712460B,0x69540EC9,
    254 0x6BA27E65,0x6DFDDBCC,0x70666F76,0x72DC8374,
    255 0x75606374,0x77F25CCE,0x7A92BE8B,0x7D41D96E
    256 // 0x80000000
    257 };
    258 
    259 /* This table is used for lookup 2^x with   */
    260 /* x in range [0...1/32[ in steps of 1/1024 */
    261 LNK_SECTION_DATA_L1 static const UINT exp2w_tab_long[32]={
    262 0x40000000,0x400B1818,0x4016321B,0x40214E0C,
    263 0x402C6BE9,0x40378BB4,0x4042AD6D,0x404DD113,
    264 0x4058F6A8,0x40641E2B,0x406F479E,0x407A7300,
    265 0x4085A051,0x4090CF92,0x409C00C4,0x40A733E6,
    266 0x40B268FA,0x40BD9FFF,0x40C8D8F5,0x40D413DD,
    267 0x40DF50B8,0x40EA8F86,0x40F5D046,0x410112FA,
    268 0x410C57A2,0x41179E3D,0x4122E6CD,0x412E3152,
    269 0x41397DCC,0x4144CC3B,0x41501CA0,0x415B6EFB,
    270 // 0x4166C34C,
    271 };
    272 /* This table is used for lookup 2^x with   */
    273 /* x in range [0...1/1024[ in steps of 1/32768 */
    274 LNK_SECTION_DATA_L1 static const UINT exp2x_tab_long[32]={
    275 0x40000000,0x400058B9,0x4000B173,0x40010A2D,
    276 0x400162E8,0x4001BBA3,0x4002145F,0x40026D1B,
    277 0x4002C5D8,0x40031E95,0x40037752,0x4003D011,
    278 0x400428CF,0x4004818E,0x4004DA4E,0x4005330E,
    279 0x40058BCE,0x4005E48F,0x40063D51,0x40069613,
    280 0x4006EED5,0x40074798,0x4007A05B,0x4007F91F,
    281 0x400851E4,0x4008AAA8,0x4009036E,0x40095C33,
    282 0x4009B4FA,0x400A0DC0,0x400A6688,0x400ABF4F,
    283 //0x400B1818
    284 };
    285 
    286 LNK_SECTION_CODE_L1 FIXP_DBL CalcInvLdData(FIXP_DBL x)
    287 {
    288   int set_zero = (x <  FL2FXCONST_DBL(-31.0/64.0))? 0 : 1;
    289   int set_max  = (x >= FL2FXCONST_DBL( 31.0/64.0)) | (x == FL2FXCONST_DBL(0.0));
    290 
    291   FIXP_SGL frac = (FIXP_SGL)(LONG)(x & 0x3FF);
    292   UINT index3 = (UINT)(LONG)(x >> 10) & 0x1F;
    293   UINT index2 = (UINT)(LONG)(x >> 15) & 0x1F;
    294   UINT index1 = (UINT)(LONG)(x >> 20) & 0x1F;
    295   int exp  = (x >  FL2FXCONST_DBL(0.0f)) ? (31 - (int)(x>>25)) : (int)(-(x>>25));
    296 
    297   UINT lookup1 = exp2_tab_long[index1]*set_zero;
    298   UINT lookup2 = exp2w_tab_long[index2];
    299   UINT lookup3 = exp2x_tab_long[index3];
    300   UINT lookup3f = lookup3 + (UINT)(LONG)fMultDiv2((FIXP_DBL)(0x0016302F),(FIXP_SGL)frac);
    301 
    302   UINT lookup12 = (UINT)(LONG)fMult((FIXP_DBL)lookup1,  (FIXP_DBL) lookup2);
    303   UINT lookup   = (UINT)(LONG)fMult((FIXP_DBL)lookup12, (FIXP_DBL) lookup3f);
    304 
    305   FIXP_DBL retVal = (lookup<<3) >> exp;
    306 
    307   if (set_max)
    308     retVal=FL2FXCONST_DBL(1.0f);
    309 
    310   return retVal;
    311 }
    312 
    313 
    314 
    315 
    316 
    317 /*****************************************************************************
    318     functionname: InitLdInt and CalcLdInt
    319     description:  Create and access table with integer LdData (0 to 193)
    320 *****************************************************************************/
    321 
    322 
    323     LNK_SECTION_CONSTDATA_L1
    324     static const FIXP_DBL ldIntCoeff[] = {
    325       0x80000001, 0x00000000, 0x02000000, 0x032b8034, 0x04000000, 0x04a4d3c2, 0x052b8034, 0x059d5da0,
    326       0x06000000, 0x06570069, 0x06a4d3c2, 0x06eb3a9f, 0x072b8034, 0x0766a009, 0x079d5da0, 0x07d053f7,
    327       0x08000000, 0x082cc7ee, 0x08570069, 0x087ef05b, 0x08a4d3c2, 0x08c8ddd4, 0x08eb3a9f, 0x090c1050,
    328       0x092b8034, 0x0949a785, 0x0966a009, 0x0982809d, 0x099d5da0, 0x09b74949, 0x09d053f7, 0x09e88c6b,
    329       0x0a000000, 0x0a16bad3, 0x0a2cc7ee, 0x0a423162, 0x0a570069, 0x0a6b3d79, 0x0a7ef05b, 0x0a92203d,
    330       0x0aa4d3c2, 0x0ab7110e, 0x0ac8ddd4, 0x0ada3f60, 0x0aeb3a9f, 0x0afbd42b, 0x0b0c1050, 0x0b1bf312,
    331       0x0b2b8034, 0x0b3abb40, 0x0b49a785, 0x0b584822, 0x0b66a009, 0x0b74b1fd, 0x0b82809d, 0x0b900e61,
    332       0x0b9d5da0, 0x0baa708f, 0x0bb74949, 0x0bc3e9ca, 0x0bd053f7, 0x0bdc899b, 0x0be88c6b, 0x0bf45e09,
    333       0x0c000000, 0x0c0b73cb, 0x0c16bad3, 0x0c21d671, 0x0c2cc7ee, 0x0c379085, 0x0c423162, 0x0c4caba8,
    334       0x0c570069, 0x0c6130af, 0x0c6b3d79, 0x0c7527b9, 0x0c7ef05b, 0x0c88983f, 0x0c92203d, 0x0c9b8926,
    335       0x0ca4d3c2, 0x0cae00d2, 0x0cb7110e, 0x0cc0052b, 0x0cc8ddd4, 0x0cd19bb0, 0x0cda3f60, 0x0ce2c97d,
    336       0x0ceb3a9f, 0x0cf39355, 0x0cfbd42b, 0x0d03fda9, 0x0d0c1050, 0x0d140ca0, 0x0d1bf312, 0x0d23c41d,
    337       0x0d2b8034, 0x0d3327c7, 0x0d3abb40, 0x0d423b08, 0x0d49a785, 0x0d510118, 0x0d584822, 0x0d5f7cff,
    338       0x0d66a009, 0x0d6db197, 0x0d74b1fd, 0x0d7ba190, 0x0d82809d, 0x0d894f75, 0x0d900e61, 0x0d96bdad,
    339       0x0d9d5da0, 0x0da3ee7f, 0x0daa708f, 0x0db0e412, 0x0db74949, 0x0dbda072, 0x0dc3e9ca, 0x0dca258e,
    340       0x0dd053f7, 0x0dd6753e, 0x0ddc899b, 0x0de29143, 0x0de88c6b, 0x0dee7b47, 0x0df45e09, 0x0dfa34e1,
    341       0x0e000000, 0x0e05bf94, 0x0e0b73cb, 0x0e111cd2, 0x0e16bad3, 0x0e1c4dfb, 0x0e21d671, 0x0e275460,
    342       0x0e2cc7ee, 0x0e323143, 0x0e379085, 0x0e3ce5d8, 0x0e423162, 0x0e477346, 0x0e4caba8, 0x0e51daa8,
    343       0x0e570069, 0x0e5c1d0b, 0x0e6130af, 0x0e663b74, 0x0e6b3d79, 0x0e7036db, 0x0e7527b9, 0x0e7a1030,
    344       0x0e7ef05b, 0x0e83c857, 0x0e88983f, 0x0e8d602e, 0x0e92203d, 0x0e96d888, 0x0e9b8926, 0x0ea03232,
    345       0x0ea4d3c2, 0x0ea96df0, 0x0eae00d2, 0x0eb28c7f, 0x0eb7110e, 0x0ebb8e96, 0x0ec0052b, 0x0ec474e4,
    346       0x0ec8ddd4, 0x0ecd4012, 0x0ed19bb0, 0x0ed5f0c4, 0x0eda3f60, 0x0ede8797, 0x0ee2c97d, 0x0ee70525,
    347       0x0eeb3a9f, 0x0eef69ff, 0x0ef39355, 0x0ef7b6b4, 0x0efbd42b, 0x0effebcd, 0x0f03fda9, 0x0f0809cf,
    348       0x0f0c1050, 0x0f10113b, 0x0f140ca0, 0x0f18028d, 0x0f1bf312, 0x0f1fde3d, 0x0f23c41d, 0x0f27a4c0,
    349       0x0f2b8034
    350     };
    351 
    352 
    353   LNK_SECTION_INITCODE
    354   void InitLdInt()
    355   {
    356     /* nothing to do! Use preinitialized logarithm table */
    357   }
    358 
    359 
    360 
    361 LNK_SECTION_CODE_L1
    362 FIXP_DBL CalcLdInt(INT i)
    363 {
    364   /* calculates ld(op)/LD_DATA_SCALING */
    365   /* op is assumed to be an integer value between 1 and 193 */
    366 
    367   FDK_ASSERT((193>0) && ((FIXP_DBL)ldIntCoeff[0]==(FIXP_DBL)0x80000001)); /* tab has to be initialized */
    368 
    369   if ((i>0)&&(i<193))
    370     return ldIntCoeff[i];
    371   else
    372   {
    373     return (0);
    374   }
    375 }
    376 
    377 
    378 /*****************************************************************************
    379 
    380     functionname: invSqrtNorm2
    381     description:  delivers 1/sqrt(op) normalized to .5...1 and the shift value of the OUTPUT
    382 
    383 *****************************************************************************/
    384 #define SQRT_BITS           7
    385 #define SQRT_VALUES       128
    386 #define SQRT_BITS_MASK   0x7f
    387 
    388 LNK_SECTION_CONSTDATA_L1
    389 static const FIXP_DBL invSqrtTab[SQRT_VALUES] = {
    390   0x5a827999, 0x5a287e03, 0x59cf8cbb, 0x5977a0ab, 0x5920b4de, 0x58cac480, 0x5875cade, 0x5821c364,
    391   0x57cea99c, 0x577c792f, 0x572b2ddf, 0x56dac38d, 0x568b3631, 0x563c81df, 0x55eea2c3, 0x55a19521,
    392   0x55555555, 0x5509dfd0, 0x54bf311a, 0x547545d0, 0x542c1aa3, 0x53e3ac5a, 0x539bf7cc, 0x5354f9e6,
    393   0x530eafa4, 0x52c91617, 0x52842a5e, 0x523fe9ab, 0x51fc513f, 0x51b95e6b, 0x51770e8e, 0x51355f19,
    394   0x50f44d88, 0x50b3d768, 0x5073fa4f, 0x5034b3e6, 0x4ff601df, 0x4fb7e1f9, 0x4f7a5201, 0x4f3d4fce,
    395   0x4f00d943, 0x4ec4ec4e, 0x4e8986e9, 0x4e4ea718, 0x4e144ae8, 0x4dda7072, 0x4da115d9, 0x4d683948,
    396   0x4d2fd8f4, 0x4cf7f31b, 0x4cc08604, 0x4c898fff, 0x4c530f64, 0x4c1d0293, 0x4be767f5, 0x4bb23df9,
    397   0x4b7d8317, 0x4b4935ce, 0x4b1554a6, 0x4ae1de2a, 0x4aaed0f0, 0x4a7c2b92, 0x4a49ecb3, 0x4a1812fa,
    398   0x49e69d16, 0x49b589bb, 0x4984d7a4, 0x49548591, 0x49249249, 0x48f4fc96, 0x48c5c34a, 0x4896e53c,
    399   0x48686147, 0x483a364c, 0x480c6331, 0x47dee6e0, 0x47b1c049, 0x4784ee5f, 0x4758701c, 0x472c447c,
    400   0x47006a80, 0x46d4e130, 0x46a9a793, 0x467ebcb9, 0x46541fb3, 0x4629cf98, 0x45ffcb80, 0x45d61289,
    401   0x45aca3d5, 0x45837e88, 0x455aa1ca, 0x45320cc8, 0x4509beb0, 0x44e1b6b4, 0x44b9f40b, 0x449275ec,
    402   0x446b3b95, 0x44444444, 0x441d8f3b, 0x43f71bbe, 0x43d0e917, 0x43aaf68e, 0x43854373, 0x435fcf14,
    403   0x433a98c5, 0x43159fdb, 0x42f0e3ae, 0x42cc6397, 0x42a81ef5, 0x42841527, 0x4260458d, 0x423caf8c,
    404   0x4219528b, 0x41f62df1, 0x41d3412a, 0x41b08ba1, 0x418e0cc7, 0x416bc40d, 0x4149b0e4, 0x4127d2c3,
    405   0x41062920, 0x40e4b374, 0x40c3713a, 0x40a261ef, 0x40818511, 0x4060da21, 0x404060a1, 0x40201814
    406 };
    407 
    408 LNK_SECTION_INITCODE
    409 void InitInvSqrtTab()
    410 {
    411   /* nothing to do !
    412      use preinitialized square root table
    413   */
    414 }
    415 
    416 
    417 
    418 #if !defined(FUNCTION_invSqrtNorm2)
    419 /*****************************************************************************
    420   delivers 1/sqrt(op) normalized to .5...1 and the shift value of the OUTPUT,
    421   i.e. the denormalized result is 1/sqrt(op) = invSqrtNorm(op) * 2^(shift)
    422   uses Newton-iteration for approximation
    423       Q(n+1) = Q(n) + Q(n) * (0.5 - 2 * V * Q(n)^2)
    424       with Q = 0.5* V ^-0.5; 0.5 <= V < 1.0
    425 *****************************************************************************/
    426 FIXP_DBL invSqrtNorm2(FIXP_DBL op, INT *shift)
    427 {
    428 
    429   FIXP_DBL val = op ;
    430   FIXP_DBL reg1, reg2, regtmp ;
    431 
    432   if (val == FL2FXCONST_DBL(0.0)) {
    433     *shift = 1 ;
    434     return((LONG)1);  /* minimum positive value */
    435   }
    436 
    437 
    438   /* normalize input, calculate shift value */
    439   FDK_ASSERT(val > FL2FXCONST_DBL(0.0));
    440   *shift = fNormz(val) - 1;  /* CountLeadingBits() is not necessary here since test value is always > 0 */
    441   val <<=*shift ;  /* normalized input V */
    442   *shift+=2 ;      /* bias for exponent */
    443 
    444   /* Newton iteration of 1/sqrt(V) */
    445   reg1 = invSqrtTab[ (INT)(val>>(DFRACT_BITS-1-(SQRT_BITS+1))) & SQRT_BITS_MASK ];
    446   reg2 = FL2FXCONST_DBL(0.0625f);   /* 0.5 >> 3 */
    447 
    448   regtmp= fPow2Div2(reg1);              /* a = Q^2 */
    449   regtmp= reg2 - fMultDiv2(regtmp, val);      /* b = 0.5 - 2 * V * Q^2 */
    450   reg1 += (fMultDiv2(regtmp, reg1)<<4);       /* Q = Q + Q*b */
    451 
    452   /* calculate the output exponent = input exp/2 */
    453   if (*shift & 0x00000001) { /* odd shift values ? */
    454     reg2 = FL2FXCONST_DBL(0.707106781186547524400844362104849f); /* 1/sqrt(2); */
    455     reg1 = fMultDiv2(reg1, reg2) << 2;
    456   }
    457 
    458   *shift = *shift>>1;
    459 
    460   return(reg1);
    461 }
    462 #endif /* !defined(FUNCTION_invSqrtNorm2) */
    463 
    464 /*****************************************************************************
    465 
    466     functionname: sqrtFixp
    467     description:  delivers sqrt(op)
    468 
    469 *****************************************************************************/
    470 FIXP_DBL sqrtFixp(FIXP_DBL op)
    471 {
    472   INT tmp_exp = 0;
    473   FIXP_DBL tmp_inv = invSqrtNorm2(op, &tmp_exp);
    474 
    475   FDK_ASSERT(tmp_exp > 0) ;
    476   return( (FIXP_DBL) ( fMultDiv2( (op<<(tmp_exp-1)), tmp_inv ) << 2 ));
    477 }
    478 
    479 
    480 #if !defined(FUNCTION_schur_div)
    481 /*****************************************************************************
    482 
    483     functionname: schur_div
    484     description:  delivers op1/op2 with op3-bit accuracy
    485 
    486 *****************************************************************************/
    487 
    488 
    489 FIXP_DBL schur_div(FIXP_DBL num, FIXP_DBL denum, INT count)
    490 {
    491     INT L_num   = (LONG)num>>1;
    492     INT L_denum = (LONG)denum>>1;
    493     INT div     = 0;
    494     INT k       = count;
    495 
    496     FDK_ASSERT (num>=(FIXP_DBL)0);
    497     FDK_ASSERT (denum>(FIXP_DBL)0);
    498     FDK_ASSERT (num <= denum);
    499 
    500     if (L_num != 0)
    501         while (--k)
    502         {
    503             div   <<= 1;
    504             L_num <<= 1;
    505             if (L_num >= L_denum)
    506             {
    507                 L_num -= L_denum;
    508                 div++;
    509             }
    510         }
    511     return (FIXP_DBL)(div << (DFRACT_BITS - count));
    512 }
    513 
    514 
    515 #endif /* !defined(FUNCTION_schur_div) */
    516 
    517 
    518 #ifndef FUNCTION_fMultNorm
    519 FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2, INT *result_e)
    520 {
    521     INT    product = 0;
    522     INT    norm_f1, norm_f2;
    523 
    524     if (  (f1 == (FIXP_DBL)0) || (f2 == (FIXP_DBL)0) ) {
    525         *result_e = 0;
    526         return (FIXP_DBL)0;
    527     }
    528     norm_f1 = CountLeadingBits(f1);
    529     f1 = f1 << norm_f1;
    530     norm_f2 = CountLeadingBits(f2);
    531     f2 = f2 << norm_f2;
    532 
    533     product = fMult(f1, f2);
    534     *result_e  = - (norm_f1 + norm_f2);
    535 
    536     return (FIXP_DBL)product;
    537 }
    538 #endif
    539 
    540 #ifndef FUNCTION_fDivNorm
    541 FIXP_DBL fDivNorm(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e)
    542 {
    543     FIXP_DBL div;
    544     INT norm_num, norm_den;
    545 
    546     FDK_ASSERT (L_num   >= (FIXP_DBL)0);
    547     FDK_ASSERT (L_denum >  (FIXP_DBL)0);
    548 
    549     if(L_num == (FIXP_DBL)0)
    550     {
    551         *result_e = 0;
    552         return ((FIXP_DBL)0);
    553     }
    554 
    555     norm_num = CountLeadingBits(L_num);
    556     L_num = L_num << norm_num;
    557     L_num = L_num >> 1;
    558     *result_e = - norm_num + 1;
    559 
    560     norm_den = CountLeadingBits(L_denum);
    561     L_denum = L_denum << norm_den;
    562     *result_e -= - norm_den;
    563 
    564     div = schur_div(L_num, L_denum, FRACT_BITS);
    565 
    566     return div;
    567 }
    568 #endif /* !FUNCTION_fDivNorm */
    569 
    570 #ifndef FUNCTION_fDivNorm
    571 FIXP_DBL fDivNorm(FIXP_DBL num, FIXP_DBL denom)
    572 {
    573     INT e;
    574     FIXP_DBL res;
    575 
    576     FDK_ASSERT (denom >= num);
    577 
    578     res = fDivNorm(num, denom, &e);
    579 
    580     /* Avoid overflow since we must output a value with exponent 0
    581        there is no other choice than saturating to almost 1.0f */
    582     if(res == (FIXP_DBL)(1<<(DFRACT_BITS-2)) && e == 1)
    583     {
    584         res = (FIXP_DBL)MAXVAL_DBL;
    585     }
    586     else
    587     {
    588         res = scaleValue(res, e);
    589     }
    590 
    591     return res;
    592 }
    593 #endif /* !FUNCTION_fDivNorm */
    594 
    595 #ifndef FUNCTION_fDivNormHighPrec
    596 FIXP_DBL fDivNormHighPrec(FIXP_DBL num, FIXP_DBL denom, INT *result_e)
    597 {
    598     FIXP_DBL div;
    599     INT norm_num, norm_den;
    600 
    601     FDK_ASSERT (num   >= (FIXP_DBL)0);
    602     FDK_ASSERT (denom >  (FIXP_DBL)0);
    603 
    604     if(num == (FIXP_DBL)0)
    605     {
    606         *result_e = 0;
    607         return ((FIXP_DBL)0);
    608     }
    609 
    610     norm_num = CountLeadingBits(num);
    611     num = num << norm_num;
    612     num = num >> 1;
    613     *result_e = - norm_num + 1;
    614 
    615     norm_den = CountLeadingBits(denom);
    616     denom = denom << norm_den;
    617     *result_e -=  - norm_den;
    618 
    619     div = schur_div(num, denom, 31);
    620     return div;
    621 }
    622 #endif /* !FUNCTION_fDivNormHighPrec */
    623 
    624 
    625 
    626 FIXP_DBL CalcLog2(FIXP_DBL base_m, INT base_e, INT *result_e)
    627 {
    628   return fLog2(base_m, base_e, result_e);
    629 }
    630 
    631 FIXP_DBL f2Pow(
    632         const FIXP_DBL exp_m, const INT exp_e,
    633         INT *result_e
    634         )
    635 {
    636   FIXP_DBL frac_part, result_m;
    637   INT int_part;
    638 
    639   if (exp_e > 0)
    640   {
    641       INT exp_bits = DFRACT_BITS-1 - exp_e;
    642       int_part = exp_m >> exp_bits;
    643       frac_part = exp_m - (FIXP_DBL)(int_part << exp_bits);
    644       frac_part = frac_part << exp_e;
    645   }
    646   else
    647   {
    648       int_part = 0;
    649       frac_part = exp_m >> -exp_e;
    650   }
    651 
    652   /* Best accuracy is around 0, so try to get there with the fractional part. */
    653   if( frac_part > FL2FXCONST_DBL(0.5f) )
    654   {
    655       int_part = int_part + 1;
    656       frac_part = frac_part + FL2FXCONST_DBL(-1.0f);
    657   }
    658   if( frac_part < FL2FXCONST_DBL(-0.5f) )
    659   {
    660       int_part = int_part - 1;
    661       frac_part = -(FL2FXCONST_DBL(-1.0f) - frac_part);
    662   }
    663 
    664   /* Evaluate taylor polynomial which approximates 2^x */
    665   {
    666     FIXP_DBL p;
    667 
    668     /* result_m ~= 2^frac_part */
    669     p = frac_part;
    670     /* First taylor series coefficient a_0 = 1.0, scaled by 0.5 due to fMultDiv2(). */
    671     result_m = FL2FXCONST_DBL(1.0f/2.0f);
    672     for (INT i = 0; i < POW2_PRECISION; i++) {
    673       /* next taylor series term: a_i * x^i, x=0 */
    674       result_m = fMultAddDiv2(result_m, pow2Coeff[i], p);
    675       p  = fMult(p, frac_part);
    676     }
    677   }
    678 
    679   /* "+ 1" compensates fMultAddDiv2() of the polynomial evaluation above. */
    680   *result_e = int_part + 1;
    681 
    682   return result_m;
    683 }
    684 
    685 FIXP_DBL f2Pow(
    686         const FIXP_DBL exp_m, const INT exp_e
    687         )
    688 {
    689   FIXP_DBL result_m;
    690   INT result_e;
    691 
    692   result_m = f2Pow(exp_m, exp_e, &result_e);
    693   result_e = fixMin(DFRACT_BITS-1,fixMax(-(DFRACT_BITS-1),result_e));
    694 
    695   return scaleValue(result_m, result_e);
    696 }
    697 
    698 FIXP_DBL fPow(
    699         FIXP_DBL base_m, INT base_e,
    700         FIXP_DBL exp_m, INT exp_e,
    701         INT *result_e
    702         )
    703 {
    704     INT ans_lg2_e, baselg2_e;
    705     FIXP_DBL base_lg2, ans_lg2, result;
    706 
    707     /* Calc log2 of base */
    708     base_lg2 = fLog2(base_m, base_e, &baselg2_e);
    709 
    710     /* Prepare exp */
    711     {
    712       INT leadingBits;
    713 
    714       leadingBits = CountLeadingBits(fAbs(exp_m));
    715       exp_m = exp_m << leadingBits;
    716       exp_e -= leadingBits;
    717     }
    718 
    719     /* Calc base pow exp */
    720     ans_lg2 = fMult(base_lg2, exp_m);
    721     ans_lg2_e = exp_e + baselg2_e;
    722 
    723     /* Calc antilog */
    724     result = f2Pow(ans_lg2, ans_lg2_e, result_e);
    725 
    726     return result;
    727 }
    728 
    729 FIXP_DBL fLdPow(
    730         FIXP_DBL baseLd_m,
    731         INT baseLd_e,
    732         FIXP_DBL exp_m, INT exp_e,
    733         INT *result_e
    734         )
    735 {
    736     INT ans_lg2_e;
    737     FIXP_DBL ans_lg2, result;
    738 
    739     /* Prepare exp */
    740     {
    741       INT leadingBits;
    742 
    743       leadingBits = CountLeadingBits(fAbs(exp_m));
    744       exp_m = exp_m << leadingBits;
    745       exp_e -= leadingBits;
    746     }
    747 
    748     /* Calc base pow exp */
    749     ans_lg2 = fMult(baseLd_m, exp_m);
    750     ans_lg2_e = exp_e + baseLd_e;
    751 
    752     /* Calc antilog */
    753     result = f2Pow(ans_lg2, ans_lg2_e, result_e);
    754 
    755     return result;
    756 }
    757 
    758 FIXP_DBL fLdPow(
    759         FIXP_DBL baseLd_m, INT baseLd_e,
    760         FIXP_DBL exp_m, INT exp_e
    761         )
    762 {
    763   FIXP_DBL result_m;
    764   int result_e;
    765 
    766   result_m = fLdPow(baseLd_m, baseLd_e, exp_m, exp_e, &result_e);
    767 
    768   return SATURATE_SHIFT(result_m, -result_e, DFRACT_BITS);
    769 }
    770 
    771 FIXP_DBL fPowInt(
    772         FIXP_DBL base_m, INT base_e,
    773         INT exp,
    774         INT *pResult_e
    775         )
    776 {
    777   FIXP_DBL result;
    778 
    779   if (exp != 0) {
    780     INT result_e = 0;
    781 
    782     if (base_m != (FIXP_DBL)0) {
    783       {
    784         INT leadingBits;
    785         leadingBits = CountLeadingBits( base_m );
    786         base_m <<= leadingBits;
    787         base_e -= leadingBits;
    788       }
    789 
    790       result = base_m;
    791 
    792       {
    793         int i;
    794         for (i = 1; i < fAbs(exp); i++) {
    795           result = fMult(result, base_m);
    796         }
    797       }
    798 
    799       if (exp < 0) {
    800         /* 1.0 / ans */
    801         result = fDivNorm( FL2FXCONST_DBL(0.5f), result, &result_e );
    802         result_e++;
    803       } else {
    804         int ansScale = CountLeadingBits( result );
    805         result <<= ansScale;
    806         result_e -= ansScale;
    807       }
    808 
    809       result_e += exp * base_e;
    810 
    811     } else {
    812       result = (FIXP_DBL)0;
    813     }
    814     *pResult_e = result_e;
    815   }
    816   else {
    817     result =  FL2FXCONST_DBL(0.5f);
    818     *pResult_e = 1;
    819   }
    820 
    821   return result;
    822 }
    823 
    824 FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e, INT *result_e)
    825 {
    826   FIXP_DBL result_m;
    827 
    828   /* Short cut for zero and negative numbers. */
    829   if ( x_m <= FL2FXCONST_DBL(0.0f) ) {
    830     *result_e = DFRACT_BITS-1;
    831     return FL2FXCONST_DBL(-1.0f);
    832   }
    833 
    834   /* Calculate log2() */
    835   {
    836     FIXP_DBL px2_m, x2_m;
    837 
    838     /* Move input value x_m * 2^x_e toward 1.0, where the taylor approximation
    839        of the function log(1-x) centered at 0 is most accurate. */
    840     {
    841       INT b_norm;
    842 
    843       b_norm = fNormz(x_m)-1;
    844       x2_m = x_m << b_norm;
    845       x_e = x_e - b_norm;
    846     }
    847 
    848     /* map x from log(x) domain to log(1-x) domain. */
    849     x2_m = - (x2_m + FL2FXCONST_DBL(-1.0) );
    850 
    851     /* Taylor polinomial approximation of ln(1-x) */
    852     result_m  = FL2FXCONST_DBL(0.0);
    853     px2_m = x2_m;
    854     for (int i=0; i<LD_PRECISION; i++) {
    855       result_m = fMultAddDiv2(result_m, ldCoeff[i], px2_m);
    856       px2_m = fMult(px2_m, x2_m);
    857     }
    858     /* Multiply result with 1/ln(2) = 1.0 + 0.442695040888 (get log2(x) from ln(x) result). */
    859     result_m = fMultAddDiv2(result_m, result_m, FL2FXCONST_DBL(2.0*0.4426950408889634073599246810019));
    860 
    861     /* Add exponent part. log2(x_m * 2^x_e) = log2(x_m) + x_e */
    862     if (x_e != 0)
    863     {
    864       int enorm;
    865 
    866       enorm = DFRACT_BITS - fNorm((FIXP_DBL)x_e);
    867       /* The -1 in the right shift of result_m compensates the fMultDiv2() above in the taylor polinomial evaluation loop.*/
    868       result_m = (result_m >> (enorm-1)) + ((FIXP_DBL)x_e << (DFRACT_BITS-1-enorm));
    869 
    870       *result_e = enorm;
    871     } else {
    872       /* 1 compensates the fMultDiv2() above in the taylor polinomial evaluation loop.*/
    873       *result_e = 1;
    874     }
    875   }
    876 
    877   return result_m;
    878 }
    879 
    880 FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e)
    881 {
    882   if ( x_m <= FL2FXCONST_DBL(0.0f) ) {
    883     x_m = FL2FXCONST_DBL(-1.0f);
    884   }
    885   else {
    886     INT result_e;
    887     x_m = fLog2(x_m, x_e, &result_e);
    888     x_m = scaleValue(x_m, result_e-LD_DATA_SHIFT);
    889   }
    890   return  x_m;
    891 }
    892 
    893 
    894 
    895 
    896