• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 
2 /* -----------------------------------------------------------------------------------------------------------
3 Software License for The Fraunhofer FDK AAC Codec Library for Android
4 
5 � Copyright  1995 - 2012 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
CalcLdData(FIXP_DBL op)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
LdDataVector(FIXP_DBL * srcVector,FIXP_DBL * destVector,INT n)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
mul_dbl_sgl_rnd(const FIXP_DBL op1,const FIXP_SGL op2)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
CalcInvLdData(FIXP_DBL op)245 FIXP_DBL CalcInvLdData(FIXP_DBL op)
246 {
247   FIXP_DBL result_m;
248 
249   if ( op == FL2FXCONST_DBL(0.0f) ) {
250     result_m = (FIXP_DBL)MAXVAL_DBL;
251   }
252   else if ( op < FL2FXCONST_DBL(0.0f) ) {
253     result_m = f2Pow(op, LD_DATA_SHIFT);
254   }
255   else {
256     int result_e;
257 
258     result_m = f2Pow(op, LD_DATA_SHIFT, &result_e);
259     result_e = fixMin(fixMax(result_e+1-(DFRACT_BITS-1), -(DFRACT_BITS-1)), (DFRACT_BITS-1)); /* rounding and saturation */
260 
261     if ( (result_e>0) && ( result_m > (((FIXP_DBL)MAXVAL_DBL)>>result_e) ) ) {
262       result_m = (FIXP_DBL)MAXVAL_DBL; /* saturate to max representable value */
263     }
264     else {
265       result_m = (scaleValue(result_m, result_e)+(FIXP_DBL)1)>>1; /* descale result + rounding */
266     }
267   }
268   return result_m;
269 }
270 
271 
272 
273 
274 
275 /*****************************************************************************
276     functionname: InitLdInt and CalcLdInt
277     description:  Create and access table with integer LdData (0 to 193)
278 *****************************************************************************/
279 
280 
281     LNK_SECTION_CONSTDATA_L1
282     static const FIXP_DBL ldIntCoeff[] = {
283       0x80000001, 0x00000000, 0x02000000, 0x032b8034, 0x04000000, 0x04a4d3c2, 0x052b8034, 0x059d5da0,
284       0x06000000, 0x06570069, 0x06a4d3c2, 0x06eb3a9f, 0x072b8034, 0x0766a009, 0x079d5da0, 0x07d053f7,
285       0x08000000, 0x082cc7ee, 0x08570069, 0x087ef05b, 0x08a4d3c2, 0x08c8ddd4, 0x08eb3a9f, 0x090c1050,
286       0x092b8034, 0x0949a785, 0x0966a009, 0x0982809d, 0x099d5da0, 0x09b74949, 0x09d053f7, 0x09e88c6b,
287       0x0a000000, 0x0a16bad3, 0x0a2cc7ee, 0x0a423162, 0x0a570069, 0x0a6b3d79, 0x0a7ef05b, 0x0a92203d,
288       0x0aa4d3c2, 0x0ab7110e, 0x0ac8ddd4, 0x0ada3f60, 0x0aeb3a9f, 0x0afbd42b, 0x0b0c1050, 0x0b1bf312,
289       0x0b2b8034, 0x0b3abb40, 0x0b49a785, 0x0b584822, 0x0b66a009, 0x0b74b1fd, 0x0b82809d, 0x0b900e61,
290       0x0b9d5da0, 0x0baa708f, 0x0bb74949, 0x0bc3e9ca, 0x0bd053f7, 0x0bdc899b, 0x0be88c6b, 0x0bf45e09,
291       0x0c000000, 0x0c0b73cb, 0x0c16bad3, 0x0c21d671, 0x0c2cc7ee, 0x0c379085, 0x0c423162, 0x0c4caba8,
292       0x0c570069, 0x0c6130af, 0x0c6b3d79, 0x0c7527b9, 0x0c7ef05b, 0x0c88983f, 0x0c92203d, 0x0c9b8926,
293       0x0ca4d3c2, 0x0cae00d2, 0x0cb7110e, 0x0cc0052b, 0x0cc8ddd4, 0x0cd19bb0, 0x0cda3f60, 0x0ce2c97d,
294       0x0ceb3a9f, 0x0cf39355, 0x0cfbd42b, 0x0d03fda9, 0x0d0c1050, 0x0d140ca0, 0x0d1bf312, 0x0d23c41d,
295       0x0d2b8034, 0x0d3327c7, 0x0d3abb40, 0x0d423b08, 0x0d49a785, 0x0d510118, 0x0d584822, 0x0d5f7cff,
296       0x0d66a009, 0x0d6db197, 0x0d74b1fd, 0x0d7ba190, 0x0d82809d, 0x0d894f75, 0x0d900e61, 0x0d96bdad,
297       0x0d9d5da0, 0x0da3ee7f, 0x0daa708f, 0x0db0e412, 0x0db74949, 0x0dbda072, 0x0dc3e9ca, 0x0dca258e,
298       0x0dd053f7, 0x0dd6753e, 0x0ddc899b, 0x0de29143, 0x0de88c6b, 0x0dee7b47, 0x0df45e09, 0x0dfa34e1,
299       0x0e000000, 0x0e05bf94, 0x0e0b73cb, 0x0e111cd2, 0x0e16bad3, 0x0e1c4dfb, 0x0e21d671, 0x0e275460,
300       0x0e2cc7ee, 0x0e323143, 0x0e379085, 0x0e3ce5d8, 0x0e423162, 0x0e477346, 0x0e4caba8, 0x0e51daa8,
301       0x0e570069, 0x0e5c1d0b, 0x0e6130af, 0x0e663b74, 0x0e6b3d79, 0x0e7036db, 0x0e7527b9, 0x0e7a1030,
302       0x0e7ef05b, 0x0e83c857, 0x0e88983f, 0x0e8d602e, 0x0e92203d, 0x0e96d888, 0x0e9b8926, 0x0ea03232,
303       0x0ea4d3c2, 0x0ea96df0, 0x0eae00d2, 0x0eb28c7f, 0x0eb7110e, 0x0ebb8e96, 0x0ec0052b, 0x0ec474e4,
304       0x0ec8ddd4, 0x0ecd4012, 0x0ed19bb0, 0x0ed5f0c4, 0x0eda3f60, 0x0ede8797, 0x0ee2c97d, 0x0ee70525,
305       0x0eeb3a9f, 0x0eef69ff, 0x0ef39355, 0x0ef7b6b4, 0x0efbd42b, 0x0effebcd, 0x0f03fda9, 0x0f0809cf,
306       0x0f0c1050, 0x0f10113b, 0x0f140ca0, 0x0f18028d, 0x0f1bf312, 0x0f1fde3d, 0x0f23c41d, 0x0f27a4c0,
307       0x0f2b8034
308     };
309 
310 
311   LNK_SECTION_INITCODE
InitLdInt()312   void InitLdInt()
313   {
314     /* nothing to do! Use preinitialized logarithm table */
315   }
316 
317 
318 
319 LNK_SECTION_CODE_L1
CalcLdInt(INT i)320 FIXP_DBL CalcLdInt(INT i)
321 {
322   /* calculates ld(op)/LD_DATA_SCALING */
323   /* op is assumed to be an integer value between 1 and 193 */
324 
325   FDK_ASSERT((193>0) && ((FIXP_DBL)ldIntCoeff[0]==(FIXP_DBL)0x80000001)); /* tab has to be initialized */
326 
327   if ((i>0)&&(i<193))
328     return ldIntCoeff[i];
329   else
330   {
331     return (0);
332   }
333 }
334 
335 
336 /*****************************************************************************
337 
338     functionname: invSqrtNorm2
339     description:  delivers 1/sqrt(op) normalized to .5...1 and the shift value of the OUTPUT
340 
341 *****************************************************************************/
342 #define SQRT_BITS           7
343 #define SQRT_VALUES       128
344 #define SQRT_BITS_MASK   0x7f
345 
346 LNK_SECTION_CONSTDATA_L1
347 static const FIXP_DBL invSqrtTab[SQRT_VALUES] = {
348   0x5a827999, 0x5a287e03, 0x59cf8cbb, 0x5977a0ab, 0x5920b4de, 0x58cac480, 0x5875cade, 0x5821c364,
349   0x57cea99c, 0x577c792f, 0x572b2ddf, 0x56dac38d, 0x568b3631, 0x563c81df, 0x55eea2c3, 0x55a19521,
350   0x55555555, 0x5509dfd0, 0x54bf311a, 0x547545d0, 0x542c1aa3, 0x53e3ac5a, 0x539bf7cc, 0x5354f9e6,
351   0x530eafa4, 0x52c91617, 0x52842a5e, 0x523fe9ab, 0x51fc513f, 0x51b95e6b, 0x51770e8e, 0x51355f19,
352   0x50f44d88, 0x50b3d768, 0x5073fa4f, 0x5034b3e6, 0x4ff601df, 0x4fb7e1f9, 0x4f7a5201, 0x4f3d4fce,
353   0x4f00d943, 0x4ec4ec4e, 0x4e8986e9, 0x4e4ea718, 0x4e144ae8, 0x4dda7072, 0x4da115d9, 0x4d683948,
354   0x4d2fd8f4, 0x4cf7f31b, 0x4cc08604, 0x4c898fff, 0x4c530f64, 0x4c1d0293, 0x4be767f5, 0x4bb23df9,
355   0x4b7d8317, 0x4b4935ce, 0x4b1554a6, 0x4ae1de2a, 0x4aaed0f0, 0x4a7c2b92, 0x4a49ecb3, 0x4a1812fa,
356   0x49e69d16, 0x49b589bb, 0x4984d7a4, 0x49548591, 0x49249249, 0x48f4fc96, 0x48c5c34a, 0x4896e53c,
357   0x48686147, 0x483a364c, 0x480c6331, 0x47dee6e0, 0x47b1c049, 0x4784ee5f, 0x4758701c, 0x472c447c,
358   0x47006a80, 0x46d4e130, 0x46a9a793, 0x467ebcb9, 0x46541fb3, 0x4629cf98, 0x45ffcb80, 0x45d61289,
359   0x45aca3d5, 0x45837e88, 0x455aa1ca, 0x45320cc8, 0x4509beb0, 0x44e1b6b4, 0x44b9f40b, 0x449275ec,
360   0x446b3b95, 0x44444444, 0x441d8f3b, 0x43f71bbe, 0x43d0e917, 0x43aaf68e, 0x43854373, 0x435fcf14,
361   0x433a98c5, 0x43159fdb, 0x42f0e3ae, 0x42cc6397, 0x42a81ef5, 0x42841527, 0x4260458d, 0x423caf8c,
362   0x4219528b, 0x41f62df1, 0x41d3412a, 0x41b08ba1, 0x418e0cc7, 0x416bc40d, 0x4149b0e4, 0x4127d2c3,
363   0x41062920, 0x40e4b374, 0x40c3713a, 0x40a261ef, 0x40818511, 0x4060da21, 0x404060a1, 0x40201814
364 };
365 
366 LNK_SECTION_INITCODE
InitInvSqrtTab()367 void InitInvSqrtTab()
368 {
369   /* nothing to do !
370      use preinitialized square root table
371   */
372 }
373 
374 
375 
376 #if !defined(FUNCTION_invSqrtNorm2)
377 /*****************************************************************************
378   delivers 1/sqrt(op) normalized to .5...1 and the shift value of the OUTPUT,
379   i.e. the denormalized result is 1/sqrt(op) = invSqrtNorm(op) * 2^(shift)
380   uses Newton-iteration for approximation
381       Q(n+1) = Q(n) + Q(n) * (0.5 - 2 * V * Q(n)^2)
382       with Q = 0.5* V ^-0.5; 0.5 <= V < 1.0
383 *****************************************************************************/
invSqrtNorm2(FIXP_DBL op,INT * shift)384 FIXP_DBL invSqrtNorm2(FIXP_DBL op, INT *shift)
385 {
386 
387   FIXP_DBL val = op ;
388   FIXP_DBL reg1, reg2, regtmp ;
389 
390   if (val == FL2FXCONST_DBL(0.0)) {
391     *shift = 1 ;
392     return((LONG)1);  /* minimum positive value */
393   }
394 
395 
396   /* normalize input, calculate shift value */
397   FDK_ASSERT(val > FL2FXCONST_DBL(0.0));
398   *shift = fNormz(val) - 1;  /* CountLeadingBits() is not necessary here since test value is always > 0 */
399   val <<=*shift ;  /* normalized input V */
400   *shift+=2 ;      /* bias for exponent */
401 
402   /* Newton iteration of 1/sqrt(V) */
403   reg1 = invSqrtTab[ (INT)(val>>(DFRACT_BITS-1-(SQRT_BITS+1))) & SQRT_BITS_MASK ];
404   reg2 = FL2FXCONST_DBL(0.0625f);   /* 0.5 >> 3 */
405 
406   regtmp= fPow2Div2(reg1);              /* a = Q^2 */
407   regtmp= reg2 - fMultDiv2(regtmp, val);      /* b = 0.5 - 2 * V * Q^2 */
408   reg1 += (fMultDiv2(regtmp, reg1)<<4);       /* Q = Q + Q*b */
409 
410   /* calculate the output exponent = input exp/2 */
411   if (*shift & 0x00000001) { /* odd shift values ? */
412     reg2 = FL2FXCONST_DBL(0.707106781186547524400844362104849f); /* 1/sqrt(2); */
413     reg1 = fMultDiv2(reg1, reg2) << 2;
414   }
415 
416   *shift = *shift>>1;
417 
418   return(reg1);
419 }
420 #endif /* !defined(FUNCTION_invSqrtNorm2) */
421 
422 /*****************************************************************************
423 
424     functionname: sqrtFixp
425     description:  delivers sqrt(op)
426 
427 *****************************************************************************/
sqrtFixp(FIXP_DBL op)428 FIXP_DBL sqrtFixp(FIXP_DBL op)
429 {
430   INT tmp_exp = 0;
431   FIXP_DBL tmp_inv = invSqrtNorm2(op, &tmp_exp);
432 
433   FDK_ASSERT(tmp_exp > 0) ;
434   return( (FIXP_DBL) ( fMultDiv2( (op<<(tmp_exp-1)), tmp_inv ) << 2 ));
435 }
436 
437 
438 #if !defined(FUNCTION_schur_div)
439 /*****************************************************************************
440 
441     functionname: schur_div
442     description:  delivers op1/op2 with op3-bit accuracy
443 
444 *****************************************************************************/
445 
446 
schur_div(FIXP_DBL num,FIXP_DBL denum,INT count)447 FIXP_DBL schur_div(FIXP_DBL num, FIXP_DBL denum, INT count)
448 {
449     INT L_num   = (LONG)num>>1;
450     INT L_denum = (LONG)denum>>1;
451     INT div     = 0;
452     INT k       = count;
453 
454     FDK_ASSERT (num>=(FIXP_DBL)0);
455     FDK_ASSERT (denum>(FIXP_DBL)0);
456     FDK_ASSERT (num <= denum);
457 
458     if (L_num != 0)
459         while (--k)
460         {
461             div   <<= 1;
462             L_num <<= 1;
463             if (L_num >= L_denum)
464             {
465                 L_num -= L_denum;
466                 div++;
467             }
468         }
469     return (FIXP_DBL)(div << (DFRACT_BITS - count));
470 }
471 
472 
473 #endif /* !defined(FUNCTION_schur_div) */
474 
475 
476 #ifndef FUNCTION_fMultNorm
fMultNorm(FIXP_DBL f1,FIXP_DBL f2,INT * result_e)477 FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2, INT *result_e)
478 {
479     INT    product = 0;
480     INT    norm_f1, norm_f2;
481 
482     if (  (f1 == (FIXP_DBL)0) || (f2 == (FIXP_DBL)0) ) {
483         *result_e = 0;
484         return (FIXP_DBL)0;
485     }
486     norm_f1 = CountLeadingBits(f1);
487     f1 = f1 << norm_f1;
488     norm_f2 = CountLeadingBits(f2);
489     f2 = f2 << norm_f2;
490 
491     product = fMult(f1, f2);
492     *result_e  = - (norm_f1 + norm_f2);
493 
494     return (FIXP_DBL)product;
495 }
496 #endif
497 
498 #ifndef FUNCTION_fDivNorm
fDivNorm(FIXP_DBL L_num,FIXP_DBL L_denum,INT * result_e)499 FIXP_DBL fDivNorm(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e)
500 {
501     FIXP_DBL div;
502     INT norm_num, norm_den;
503 
504     FDK_ASSERT (L_num   >= (FIXP_DBL)0);
505     FDK_ASSERT (L_denum >  (FIXP_DBL)0);
506 
507     if(L_num == (FIXP_DBL)0)
508     {
509         *result_e = 0;
510         return ((FIXP_DBL)0);
511     }
512 
513     norm_num = CountLeadingBits(L_num);
514     L_num = L_num << norm_num;
515     L_num = L_num >> 1;
516     *result_e = - norm_num + 1;
517 
518     norm_den = CountLeadingBits(L_denum);
519     L_denum = L_denum << norm_den;
520     *result_e -= - norm_den;
521 
522     div = schur_div(L_num, L_denum, FRACT_BITS);
523 
524     return div;
525 }
526 #endif /* !FUNCTION_fDivNorm */
527 
528 #ifndef FUNCTION_fDivNorm
fDivNorm(FIXP_DBL num,FIXP_DBL denom)529 FIXP_DBL fDivNorm(FIXP_DBL num, FIXP_DBL denom)
530 {
531     INT e;
532     FIXP_DBL res;
533 
534     FDK_ASSERT (denom >= num);
535 
536     res = fDivNorm(num, denom, &e);
537 
538     /* Avoid overflow since we must output a value with exponent 0
539        there is no other choice than saturating to almost 1.0f */
540     if(res == (FIXP_DBL)(1<<(DFRACT_BITS-2)) && e == 1)
541     {
542         res = (FIXP_DBL)MAXVAL_DBL;
543     }
544     else
545     {
546         res = scaleValue(res, e);
547     }
548 
549     return res;
550 }
551 #endif /* !FUNCTION_fDivNorm */
552 
553 #ifndef FUNCTION_fDivNormHighPrec
fDivNormHighPrec(FIXP_DBL num,FIXP_DBL denom,INT * result_e)554 FIXP_DBL fDivNormHighPrec(FIXP_DBL num, FIXP_DBL denom, INT *result_e)
555 {
556     FIXP_DBL div;
557     INT norm_num, norm_den;
558 
559     FDK_ASSERT (num   >= (FIXP_DBL)0);
560     FDK_ASSERT (denom >  (FIXP_DBL)0);
561 
562     if(num == (FIXP_DBL)0)
563     {
564         *result_e = 0;
565         return ((FIXP_DBL)0);
566     }
567 
568     norm_num = CountLeadingBits(num);
569     num = num << norm_num;
570     num = num >> 1;
571     *result_e = - norm_num + 1;
572 
573     norm_den = CountLeadingBits(denom);
574     denom = denom << norm_den;
575     *result_e -=  - norm_den;
576 
577     div = schur_div(num, denom, 31);
578     return div;
579 }
580 #endif /* !FUNCTION_fDivNormHighPrec */
581 
582 
583 
CalcLog2(FIXP_DBL base_m,INT base_e,INT * result_e)584 FIXP_DBL CalcLog2(FIXP_DBL base_m, INT base_e, INT *result_e)
585 {
586   return fLog2(base_m, base_e, result_e);
587 }
588 
f2Pow(const FIXP_DBL exp_m,const INT exp_e,INT * result_e)589 FIXP_DBL f2Pow(
590         const FIXP_DBL exp_m, const INT exp_e,
591         INT *result_e
592         )
593 {
594   FIXP_DBL frac_part, result_m;
595   INT int_part;
596 
597   if (exp_e > 0)
598   {
599       INT exp_bits = DFRACT_BITS-1 - exp_e;
600       int_part = exp_m >> exp_bits;
601       frac_part = exp_m - (FIXP_DBL)(int_part << exp_bits);
602       frac_part = frac_part << exp_e;
603   }
604   else
605   {
606       int_part = 0;
607       frac_part = exp_m >> -exp_e;
608   }
609 
610   /* Best accuracy is around 0, so try to get there with the fractional part. */
611   if( frac_part > FL2FXCONST_DBL(0.5f) )
612   {
613       int_part = int_part + 1;
614       frac_part = frac_part + FL2FXCONST_DBL(-1.0f);
615   }
616   if( frac_part < FL2FXCONST_DBL(-0.5f) )
617   {
618       int_part = int_part - 1;
619       frac_part = -(FL2FXCONST_DBL(-1.0f) - frac_part);
620   }
621 
622   /* Evaluate taylor polynomial which approximates 2^x */
623   {
624     FIXP_DBL p;
625 
626     /* result_m ~= 2^frac_part */
627     p = frac_part;
628     /* First taylor series coefficient a_0 = 1.0, scaled by 0.5 due to fMultDiv2(). */
629     result_m = FL2FXCONST_DBL(1.0f/2.0f);
630     for (INT i = 0; i < POW2_PRECISION; i++) {
631       /* next taylor series term: a_i * x^i, x=0 */
632       result_m = fMultAddDiv2(result_m, pow2Coeff[i], p);
633       p  = fMult(p, frac_part);
634     }
635   }
636 
637   /* "+ 1" compensates fMultAddDiv2() of the polynomial evaluation above. */
638   *result_e = int_part + 1;
639 
640   return result_m;
641 }
642 
f2Pow(const FIXP_DBL exp_m,const INT exp_e)643 FIXP_DBL f2Pow(
644         const FIXP_DBL exp_m, const INT exp_e
645         )
646 {
647   FIXP_DBL result_m;
648   INT result_e;
649 
650   result_m = f2Pow(exp_m, exp_e, &result_e);
651   result_e = fixMin(DFRACT_BITS-1,fixMax(-(DFRACT_BITS-1),result_e));
652 
653   return scaleValue(result_m, result_e);
654 }
655 
fPow(FIXP_DBL base_m,INT base_e,FIXP_DBL exp_m,INT exp_e,INT * result_e)656 FIXP_DBL fPow(
657         FIXP_DBL base_m, INT base_e,
658         FIXP_DBL exp_m, INT exp_e,
659         INT *result_e
660         )
661 {
662     INT ans_lg2_e, baselg2_e;
663     FIXP_DBL base_lg2, ans_lg2, result;
664 
665     /* Calc log2 of base */
666     base_lg2 = fLog2(base_m, base_e, &baselg2_e);
667 
668     /* Prepare exp */
669     {
670       INT leadingBits;
671 
672       leadingBits = CountLeadingBits(fAbs(exp_m));
673       exp_m = exp_m << leadingBits;
674       exp_e -= leadingBits;
675     }
676 
677     /* Calc base pow exp */
678     ans_lg2 = fMult(base_lg2, exp_m);
679     ans_lg2_e = exp_e + baselg2_e;
680 
681     /* Calc antilog */
682     result = f2Pow(ans_lg2, ans_lg2_e, result_e);
683 
684     return result;
685 }
686 
fLdPow(FIXP_DBL baseLd_m,INT baseLd_e,FIXP_DBL exp_m,INT exp_e,INT * result_e)687 FIXP_DBL fLdPow(
688         FIXP_DBL baseLd_m,
689         INT baseLd_e,
690         FIXP_DBL exp_m, INT exp_e,
691         INT *result_e
692         )
693 {
694     INT ans_lg2_e;
695     FIXP_DBL ans_lg2, result;
696 
697     /* Prepare exp */
698     {
699       INT leadingBits;
700 
701       leadingBits = CountLeadingBits(fAbs(exp_m));
702       exp_m = exp_m << leadingBits;
703       exp_e -= leadingBits;
704     }
705 
706     /* Calc base pow exp */
707     ans_lg2 = fMult(baseLd_m, exp_m);
708     ans_lg2_e = exp_e + baseLd_e;
709 
710     /* Calc antilog */
711     result = f2Pow(ans_lg2, ans_lg2_e, result_e);
712 
713     return result;
714 }
715 
fLdPow(FIXP_DBL baseLd_m,INT baseLd_e,FIXP_DBL exp_m,INT exp_e)716 FIXP_DBL fLdPow(
717         FIXP_DBL baseLd_m, INT baseLd_e,
718         FIXP_DBL exp_m, INT exp_e
719         )
720 {
721   FIXP_DBL result_m;
722   int result_e;
723 
724   result_m = fLdPow(baseLd_m, baseLd_e, exp_m, exp_e, &result_e);
725 
726   return SATURATE_SHIFT(result_m, -result_e, DFRACT_BITS);
727 }
728 
fPowInt(FIXP_DBL base_m,INT base_e,INT exp,INT * pResult_e)729 FIXP_DBL fPowInt(
730         FIXP_DBL base_m, INT base_e,
731         INT exp,
732         INT *pResult_e
733         )
734 {
735   FIXP_DBL result;
736 
737   if (exp != 0) {
738     INT result_e = 0;
739 
740     if (base_m != (FIXP_DBL)0) {
741       {
742         INT leadingBits;
743         leadingBits = CountLeadingBits( base_m );
744         base_m <<= leadingBits;
745         base_e -= leadingBits;
746       }
747 
748       result = base_m;
749 
750       {
751         int i;
752         for (i = 1; i < fAbs(exp); i++) {
753           result = fMult(result, base_m);
754         }
755       }
756 
757       if (exp < 0) {
758         /* 1.0 / ans */
759         result = fDivNorm( FL2FXCONST_DBL(0.5f), result, &result_e );
760         result_e++;
761       } else {
762         int ansScale = CountLeadingBits( result );
763         result <<= ansScale;
764         result_e -= ansScale;
765       }
766 
767       result_e += exp * base_e;
768 
769     } else {
770       result = (FIXP_DBL)0;
771     }
772     *pResult_e = result_e;
773   }
774   else {
775     result =  FL2FXCONST_DBL(0.5f);
776     *pResult_e = 1;
777   }
778 
779   return result;
780 }
781 
fLog2(FIXP_DBL x_m,INT x_e,INT * result_e)782 FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e, INT *result_e)
783 {
784   FIXP_DBL result_m;
785 
786   /* Short cut for zero and negative numbers. */
787   if ( x_m <= FL2FXCONST_DBL(0.0f) ) {
788     *result_e = DFRACT_BITS-1;
789     return FL2FXCONST_DBL(-1.0f);
790   }
791 
792   /* Calculate log2() */
793   {
794     FIXP_DBL px2_m, x2_m;
795 
796     /* Move input value x_m * 2^x_e toward 1.0, where the taylor approximation
797        of the function log(1-x) centered at 0 is most accurate. */
798     {
799       INT b_norm;
800 
801       b_norm = fNormz(x_m)-1;
802       x2_m = x_m << b_norm;
803       x_e = x_e - b_norm;
804     }
805 
806     /* map x from log(x) domain to log(1-x) domain. */
807     x2_m = - (x2_m + FL2FXCONST_DBL(-1.0) );
808 
809     /* Taylor polinomial approximation of ln(1-x) */
810     result_m  = FL2FXCONST_DBL(0.0);
811     px2_m = x2_m;
812     for (int i=0; i<LD_PRECISION; i++) {
813       result_m = fMultAddDiv2(result_m, ldCoeff[i], px2_m);
814       px2_m = fMult(px2_m, x2_m);
815     }
816     /* Multiply result with 1/ln(2) = 1.0 + 0.442695040888 (get log2(x) from ln(x) result). */
817     result_m = fMultAddDiv2(result_m, result_m, FL2FXCONST_DBL(2.0*0.4426950408889634073599246810019));
818 
819     /* Add exponent part. log2(x_m * 2^x_e) = log2(x_m) + x_e */
820     if (x_e != 0)
821     {
822       int enorm;
823 
824       enorm = DFRACT_BITS - fNorm((FIXP_DBL)x_e);
825       /* The -1 in the right shift of result_m compensates the fMultDiv2() above in the taylor polinomial evaluation loop.*/
826       result_m = (result_m >> (enorm-1)) + ((FIXP_DBL)x_e << (DFRACT_BITS-1-enorm));
827 
828       *result_e = enorm;
829     } else {
830       /* 1 compensates the fMultDiv2() above in the taylor polinomial evaluation loop.*/
831       *result_e = 1;
832     }
833   }
834 
835   return result_m;
836 }
837 
fLog2(FIXP_DBL x_m,INT x_e)838 FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e)
839 {
840   if ( x_m <= FL2FXCONST_DBL(0.0f) ) {
841     x_m = FL2FXCONST_DBL(-1.0f);
842   }
843   else {
844     INT result_e;
845     x_m = fLog2(x_m, x_e, &result_e);
846     x_m = scaleValue(x_m, result_e-LD_DATA_SHIFT);
847   }
848   return  x_m;
849 }
850 
851 
852 
853 
854