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
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
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
CalcInvLdData(FIXP_DBL x)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
InitLdInt()354 void InitLdInt()
355 {
356 /* nothing to do! Use preinitialized logarithm table */
357 }
358
359
360
361 LNK_SECTION_CODE_L1
CalcLdInt(INT i)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
InitInvSqrtTab()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 *****************************************************************************/
invSqrtNorm2(FIXP_DBL op,INT * shift)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 *****************************************************************************/
sqrtFixp(FIXP_DBL op)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
schur_div(FIXP_DBL num,FIXP_DBL denum,INT count)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
fMultNorm(FIXP_DBL f1,FIXP_DBL f2,INT * result_e)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
fDivNorm(FIXP_DBL L_num,FIXP_DBL L_denum,INT * result_e)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
fDivNorm(FIXP_DBL num,FIXP_DBL denom)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
fDivNormHighPrec(FIXP_DBL num,FIXP_DBL denom,INT * result_e)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
CalcLog2(FIXP_DBL base_m,INT base_e,INT * result_e)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
f2Pow(const FIXP_DBL exp_m,const INT exp_e,INT * result_e)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
f2Pow(const FIXP_DBL exp_m,const INT exp_e)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
fPow(FIXP_DBL base_m,INT base_e,FIXP_DBL exp_m,INT exp_e,INT * result_e)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
fLdPow(FIXP_DBL baseLd_m,INT baseLd_e,FIXP_DBL exp_m,INT exp_e,INT * result_e)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
fLdPow(FIXP_DBL baseLd_m,INT baseLd_e,FIXP_DBL exp_m,INT exp_e)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
fPowInt(FIXP_DBL base_m,INT base_e,INT exp,INT * pResult_e)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
fLog2(FIXP_DBL x_m,INT x_e,INT * result_e)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
fLog2(FIXP_DBL x_m,INT x_e)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