libFDK/src/fixpoint_math.cpp

   1
   2 /* -----------------------------------------------------------------------------------------------------------
   3 Software License for The Fraunhofer FDK AAC Codec Library for Android
   4
   5 © Copyright  1995 - 2013 Fraunhofer-Gesellschaft zur Förderung 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@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