1 /* -----------------------------------------------------------------------------
2 Software License for The Fraunhofer FDK AAC Codec Library for Android
3
4 © Copyright 1995 - 2019 Fraunhofer-Gesellschaft zur Förderung der angewandten
5 Forschung e.V. All rights reserved.
6
7 1. INTRODUCTION
8 The Fraunhofer FDK AAC Codec Library for Android ("FDK AAC Codec") is software
9 that implements the MPEG Advanced Audio Coding ("AAC") encoding and decoding
10 scheme for digital audio. This FDK AAC Codec software is intended to be used on
11 a wide variety of Android devices.
12
13 AAC's HE-AAC and HE-AAC v2 versions are regarded as today's most efficient
14 general perceptual audio codecs. AAC-ELD is considered the best-performing
15 full-bandwidth communications codec by independent studies and is widely
16 deployed. AAC has been standardized by ISO and IEC as part of the MPEG
17 specifications.
18
19 Patent licenses for necessary patent claims for the FDK AAC Codec (including
20 those of Fraunhofer) may be obtained through Via Licensing
21 (www.vialicensing.com) or through the respective patent owners individually for
22 the purpose of encoding or decoding bit streams in products that are compliant
23 with the ISO/IEC MPEG audio standards. Please note that most manufacturers of
24 Android devices already license these patent claims through Via Licensing or
25 directly from the patent owners, and therefore FDK AAC Codec software may
26 already be covered under those patent licenses when it is used for those
27 licensed purposes only.
28
29 Commercially-licensed AAC software libraries, including floating-point versions
30 with enhanced sound quality, are also available from Fraunhofer. Users are
31 encouraged to check the Fraunhofer website for additional applications
32 information and documentation.
33
34 2. COPYRIGHT LICENSE
35
36 Redistribution and use in source and binary forms, with or without modification,
37 are permitted without payment of copyright license fees provided that you
38 satisfy the following conditions:
39
40 You must retain the complete text of this software license in redistributions of
41 the FDK AAC Codec or your modifications thereto in source code form.
42
43 You must retain the complete text of this software license in the documentation
44 and/or other materials provided with redistributions of the FDK AAC Codec or
45 your modifications thereto in binary form. You must make available free of
46 charge copies of the complete source code of the FDK AAC Codec and your
47 modifications thereto to recipients of copies in binary form.
48
49 The name of Fraunhofer may not be used to endorse or promote products derived
50 from this library without prior written permission.
51
52 You may not charge copyright license fees for anyone to use, copy or distribute
53 the FDK AAC Codec software or your modifications thereto.
54
55 Your modified versions of the FDK AAC Codec must carry prominent notices stating
56 that you changed the software and the date of any change. For modified versions
57 of the FDK AAC Codec, the term "Fraunhofer FDK AAC Codec Library for Android"
58 must be replaced by the term "Third-Party Modified Version of the Fraunhofer FDK
59 AAC Codec Library for Android."
60
61 3. NO PATENT LICENSE
62
63 NO EXPRESS OR IMPLIED LICENSES TO ANY PATENT CLAIMS, including without
64 limitation the patents of Fraunhofer, ARE GRANTED BY THIS SOFTWARE LICENSE.
65 Fraunhofer provides no warranty of patent non-infringement with respect to this
66 software.
67
68 You may use this FDK AAC Codec software or modifications thereto only for
69 purposes that are authorized by appropriate patent licenses.
70
71 4. DISCLAIMER
72
73 This FDK AAC Codec software is provided by Fraunhofer on behalf of the copyright
74 holders and contributors "AS IS" and WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES,
75 including but not limited to the implied warranties of merchantability and
76 fitness for a particular purpose. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
77 CONTRIBUTORS BE LIABLE for any direct, indirect, incidental, special, exemplary,
78 or consequential damages, including but not limited to procurement of substitute
79 goods or services; loss of use, data, or profits, or business interruption,
80 however caused and on any theory of liability, whether in contract, strict
81 liability, or tort (including negligence), arising in any way out of the use of
82 this software, even if advised of the possibility of such damage.
83
84 5. CONTACT INFORMATION
85
86 Fraunhofer Institute for Integrated Circuits IIS
87 Attention: Audio and Multimedia Departments - FDK AAC LL
88 Am Wolfsmantel 33
89 91058 Erlangen, Germany
90
91 www.iis.fraunhofer.de/amm
92 amm-info@iis.fraunhofer.de
93 ----------------------------------------------------------------------------- */
94
95 /******************* Library for basic calculation routines ********************
96
97 Author(s): M. Gayer
98
99 Description: Fixed point specific mathematical functions
100
101 *******************************************************************************/
102
103 #ifndef FIXPOINT_MATH_H
104 #define FIXPOINT_MATH_H
105
106 #include "common_fix.h"
107 #include "scale.h"
108
109 /*
110 * Data definitions
111 */
112
113 #define LD_DATA_SCALING (64.0f)
114 #define LD_DATA_SHIFT 6 /* pow(2, LD_DATA_SHIFT) = LD_DATA_SCALING */
115
116 #define MAX_LD_PRECISION 10
117 #define LD_PRECISION 10
118
119 /* Taylor series coefficients for ln(1-x), centered at 0 (MacLaurin polynomial).
120 */
121 #ifndef LDCOEFF_16BIT
122 LNK_SECTION_CONSTDATA_L1
123 static const FIXP_DBL ldCoeff[MAX_LD_PRECISION] = {
124 FL2FXCONST_DBL(-1.0), FL2FXCONST_DBL(-1.0 / 2.0),
125 FL2FXCONST_DBL(-1.0 / 3.0), FL2FXCONST_DBL(-1.0 / 4.0),
126 FL2FXCONST_DBL(-1.0 / 5.0), FL2FXCONST_DBL(-1.0 / 6.0),
127 FL2FXCONST_DBL(-1.0 / 7.0), FL2FXCONST_DBL(-1.0 / 8.0),
128 FL2FXCONST_DBL(-1.0 / 9.0), FL2FXCONST_DBL(-1.0 / 10.0)};
129 #else /* LDCOEFF_16BIT */
130 LNK_SECTION_CONSTDATA_L1
131 static const FIXP_SGL ldCoeff[MAX_LD_PRECISION] = {
132 FL2FXCONST_SGL(-1.0), FL2FXCONST_SGL(-1.0 / 2.0),
133 FL2FXCONST_SGL(-1.0 / 3.0), FL2FXCONST_SGL(-1.0 / 4.0),
134 FL2FXCONST_SGL(-1.0 / 5.0), FL2FXCONST_SGL(-1.0 / 6.0),
135 FL2FXCONST_SGL(-1.0 / 7.0), FL2FXCONST_SGL(-1.0 / 8.0),
136 FL2FXCONST_SGL(-1.0 / 9.0), FL2FXCONST_SGL(-1.0 / 10.0)};
137 #endif /* LDCOEFF_16BIT */
138
139 /*****************************************************************************
140
141 functionname: invSqrtNorm2
142 description: delivers 1/sqrt(op) normalized to .5...1 and the shift value
143 of the OUTPUT
144
145 *****************************************************************************/
146 #define SQRT_BITS 7
147 #define SQRT_VALUES (128 + 2)
148 #define SQRT_BITS_MASK 0x7f
149 #define SQRT_FRACT_BITS_MASK 0x007FFFFF
150
151 extern const FIXP_DBL invSqrtTab[SQRT_VALUES];
152
153 /*
154 * Hardware specific implementations
155 */
156
157 #if defined(__x86__)
158 #include "x86/fixpoint_math_x86.h"
159 #endif /* target architecture selector */
160
161 /*
162 * Fallback implementations
163 */
164 #if !defined(FUNCTION_fIsLessThan)
165 /**
166 * \brief Compares two fixpoint values incl. scaling.
167 * \param a_m mantissa of the first input value.
168 * \param a_e exponent of the first input value.
169 * \param b_m mantissa of the second input value.
170 * \param b_e exponent of the second input value.
171 * \return non-zero if (a_m*2^a_e) < (b_m*2^b_e), 0 otherwise
172 */
fIsLessThan(FIXP_DBL a_m,INT a_e,FIXP_DBL b_m,INT b_e)173 FDK_INLINE INT fIsLessThan(FIXP_DBL a_m, INT a_e, FIXP_DBL b_m, INT b_e) {
174 INT n;
175
176 n = fixnorm_D(a_m);
177 a_m <<= n;
178 a_e -= n;
179
180 n = fixnorm_D(b_m);
181 b_m <<= n;
182 b_e -= n;
183
184 if (a_m == (FIXP_DBL)0) a_e = b_e;
185 if (b_m == (FIXP_DBL)0) b_e = a_e;
186
187 if (a_e > b_e) {
188 return ((b_m >> fMin(a_e - b_e, DFRACT_BITS - 1)) > a_m);
189 } else {
190 return ((a_m >> fMin(b_e - a_e, DFRACT_BITS - 1)) < b_m);
191 }
192 }
193
fIsLessThan(FIXP_SGL a_m,INT a_e,FIXP_SGL b_m,INT b_e)194 FDK_INLINE INT fIsLessThan(FIXP_SGL a_m, INT a_e, FIXP_SGL b_m, INT b_e) {
195 INT n;
196
197 n = fixnorm_S(a_m);
198 a_m <<= n;
199 a_e -= n;
200
201 n = fixnorm_S(b_m);
202 b_m <<= n;
203 b_e -= n;
204
205 if (a_m == (FIXP_SGL)0) a_e = b_e;
206 if (b_m == (FIXP_SGL)0) b_e = a_e;
207
208 if (a_e > b_e) {
209 return ((b_m >> fMin(a_e - b_e, FRACT_BITS - 1)) > a_m);
210 } else {
211 return ((a_m >> fMin(b_e - a_e, FRACT_BITS - 1)) < b_m);
212 }
213 }
214 #endif
215
216 /**
217 * \brief deprecated. Use fLog2() instead.
218 */
219 #define CalcLdData(op) fLog2(op, 0)
220
221 void LdDataVector(FIXP_DBL *srcVector, FIXP_DBL *destVector, INT number);
222
223 extern const UINT exp2_tab_long[32];
224 extern const UINT exp2w_tab_long[32];
225 extern const UINT exp2x_tab_long[32];
226
227 LNK_SECTION_CODE_L1
CalcInvLdData(const FIXP_DBL x)228 FDK_INLINE FIXP_DBL CalcInvLdData(const FIXP_DBL x) {
229 int set_zero = (x < FL2FXCONST_DBL(-31.0 / 64.0)) ? 0 : 1;
230 int set_max = (x >= FL2FXCONST_DBL(31.0 / 64.0)) | (x == FL2FXCONST_DBL(0.0));
231
232 FIXP_SGL frac = (FIXP_SGL)((LONG)x & 0x3FF);
233 UINT index3 = (UINT)(LONG)(x >> 10) & 0x1F;
234 UINT index2 = (UINT)(LONG)(x >> 15) & 0x1F;
235 UINT index1 = (UINT)(LONG)(x >> 20) & 0x1F;
236 int exp = fMin(31, ((x > FL2FXCONST_DBL(0.0f)) ? (31 - (int)(x >> 25))
237 : (int)(-(x >> 25))));
238
239 UINT lookup1 = exp2_tab_long[index1] * set_zero;
240 UINT lookup2 = exp2w_tab_long[index2];
241 UINT lookup3 = exp2x_tab_long[index3];
242 UINT lookup3f =
243 lookup3 + (UINT)(LONG)fMultDiv2((FIXP_DBL)(0x0016302F), (FIXP_SGL)frac);
244
245 UINT lookup12 = (UINT)(LONG)fMult((FIXP_DBL)lookup1, (FIXP_DBL)lookup2);
246 UINT lookup = (UINT)(LONG)fMult((FIXP_DBL)lookup12, (FIXP_DBL)lookup3f);
247
248 FIXP_DBL retVal = (lookup << 3) >> exp;
249
250 if (set_max) {
251 retVal = (FIXP_DBL)MAXVAL_DBL;
252 }
253
254 return retVal;
255 }
256
257 void InitLdInt();
258 FIXP_DBL CalcLdInt(INT i);
259
260 extern const USHORT sqrt_tab[49];
261
sqrtFixp_lookup(FIXP_DBL x)262 inline FIXP_DBL sqrtFixp_lookup(FIXP_DBL x) {
263 UINT y = (INT)x;
264 UCHAR is_zero = (y == 0);
265 INT zeros = fixnormz_D(y) & 0x1e;
266 y <<= zeros;
267 UINT idx = (y >> 26) - 16;
268 USHORT frac = (y >> 10) & 0xffff;
269 USHORT nfrac = 0xffff ^ frac;
270 UINT t = (UINT)nfrac * sqrt_tab[idx] + (UINT)frac * sqrt_tab[idx + 1];
271 t = t >> (zeros >> 1);
272 return (is_zero ? 0 : t);
273 }
274
sqrtFixp_lookup(FIXP_DBL x,INT * x_e)275 inline FIXP_DBL sqrtFixp_lookup(FIXP_DBL x, INT *x_e) {
276 UINT y = (INT)x;
277 INT e;
278
279 if (x == (FIXP_DBL)0) {
280 return x;
281 }
282
283 /* Normalize */
284 e = fixnormz_D(y);
285 y <<= e;
286 e = *x_e - e + 2;
287
288 /* Correct odd exponent. */
289 if (e & 1) {
290 y >>= 1;
291 e++;
292 }
293 /* Get square root */
294 UINT idx = (y >> 26) - 16;
295 USHORT frac = (y >> 10) & 0xffff;
296 USHORT nfrac = 0xffff ^ frac;
297 UINT t = (UINT)nfrac * sqrt_tab[idx] + (UINT)frac * sqrt_tab[idx + 1];
298
299 /* Write back exponent */
300 *x_e = e >> 1;
301 return (FIXP_DBL)(LONG)(t >> 1);
302 }
303
304 void InitInvSqrtTab();
305
306 #ifndef FUNCTION_invSqrtNorm2
307 /**
308 * \brief calculate 1.0/sqrt(op)
309 * \param op_m mantissa of input value.
310 * \param result_e pointer to return the exponent of the result
311 * \return mantissa of the result
312 */
313 /*****************************************************************************
314 delivers 1/sqrt(op) normalized to .5...1 and the shift value of the OUTPUT,
315 i.e. the denormalized result is 1/sqrt(op) = invSqrtNorm(op) * 2^(shift)
316 uses Newton-iteration for approximation
317 Q(n+1) = Q(n) + Q(n) * (0.5 - 2 * V * Q(n)^2)
318 with Q = 0.5* V ^-0.5; 0.5 <= V < 1.0
319 *****************************************************************************/
invSqrtNorm2(FIXP_DBL op,INT * shift)320 static FDK_FORCEINLINE FIXP_DBL invSqrtNorm2(FIXP_DBL op, INT *shift) {
321 FIXP_DBL val = op;
322 FIXP_DBL reg1, reg2;
323
324 if (val == FL2FXCONST_DBL(0.0)) {
325 *shift = 16;
326 return ((LONG)MAXVAL_DBL); /* maximum positive value */
327 }
328
329 #define INVSQRTNORM2_LINEAR_INTERPOLATE
330 #define INVSQRTNORM2_LINEAR_INTERPOLATE_HQ
331
332 /* normalize input, calculate shift value */
333 FDK_ASSERT(val > FL2FXCONST_DBL(0.0));
334 *shift = fNormz(val) - 1; /* CountLeadingBits() is not necessary here since
335 test value is always > 0 */
336 val <<= *shift; /* normalized input V */
337 *shift += 2; /* bias for exponent */
338
339 #if defined(INVSQRTNORM2_LINEAR_INTERPOLATE)
340 INT index =
341 (INT)(val >> (DFRACT_BITS - 1 - (SQRT_BITS + 1))) & SQRT_BITS_MASK;
342 FIXP_DBL Fract =
343 (FIXP_DBL)(((INT)val & SQRT_FRACT_BITS_MASK) << (SQRT_BITS + 1));
344 FIXP_DBL diff = invSqrtTab[index + 1] - invSqrtTab[index];
345 reg1 = invSqrtTab[index] + (fMultDiv2(diff, Fract) << 1);
346 #if defined(INVSQRTNORM2_LINEAR_INTERPOLATE_HQ)
347 /* reg1 = t[i] + (t[i+1]-t[i])*fract ... already computed ...
348 + (1-fract)fract*(t[i+2]-t[i+1])/2 */
349 if (Fract != (FIXP_DBL)0) {
350 /* fract = fract * (1 - fract) */
351 Fract = fMultDiv2(Fract, (FIXP_DBL)((ULONG)0x80000000 - (ULONG)Fract)) << 1;
352 diff = diff - (invSqrtTab[index + 2] - invSqrtTab[index + 1]);
353 reg1 = fMultAddDiv2(reg1, Fract, diff);
354 }
355 #endif /* INVSQRTNORM2_LINEAR_INTERPOLATE_HQ */
356 #else
357 #error \
358 "Either define INVSQRTNORM2_NEWTON_ITERATE or INVSQRTNORM2_LINEAR_INTERPOLATE"
359 #endif
360 /* calculate the output exponent = input exp/2 */
361 if (*shift & 0x00000001) { /* odd shift values ? */
362 /* Note: Do not use rounded value 0x5A82799A to avoid overflow with
363 * shift-by-2 */
364 reg2 = (FIXP_DBL)0x5A827999;
365 /* FL2FXCONST_DBL(0.707106781186547524400844362104849f);*/ /* 1/sqrt(2);
366 */
367 reg1 = fMultDiv2(reg1, reg2) << 2;
368 }
369
370 *shift = *shift >> 1;
371
372 return (reg1);
373 }
374 #endif /* FUNCTION_invSqrtNorm2 */
375
376 #ifndef FUNCTION_sqrtFixp
sqrtFixp(FIXP_DBL op)377 static FDK_FORCEINLINE FIXP_DBL sqrtFixp(FIXP_DBL op) {
378 INT tmp_exp = 0;
379 FIXP_DBL tmp_inv = invSqrtNorm2(op, &tmp_exp);
380
381 FDK_ASSERT(tmp_exp > 0);
382 return ((FIXP_DBL)(fMultDiv2((op << (tmp_exp - 1)), tmp_inv) << 2));
383 }
384 #endif /* FUNCTION_sqrtFixp */
385
386 #ifndef FUNCTION_invFixp
387 /**
388 * \brief calculate 1.0/op
389 * \param op mantissa of the input value.
390 * \return mantissa of the result with implicit exponent of 31
391 * \exceptions are provided for op=0,1 setting max. positive value
392 */
invFixp(FIXP_DBL op)393 static inline FIXP_DBL invFixp(FIXP_DBL op) {
394 if ((op == (FIXP_DBL)0x00000000) || (op == (FIXP_DBL)0x00000001)) {
395 return ((LONG)MAXVAL_DBL);
396 }
397 INT tmp_exp;
398 FIXP_DBL tmp_inv = invSqrtNorm2(op, &tmp_exp);
399 FDK_ASSERT((31 - (2 * tmp_exp + 1)) >= 0);
400 int shift = 31 - (2 * tmp_exp + 1);
401 tmp_inv = fPow2Div2(tmp_inv);
402 if (shift) {
403 tmp_inv = ((tmp_inv >> (shift - 1)) + (FIXP_DBL)1) >> 1;
404 }
405 return tmp_inv;
406 }
407
408 /**
409 * \brief calculate 1.0/(op_m * 2^op_e)
410 * \param op_m mantissa of the input value.
411 * \param op_e pointer into were the exponent of the input value is stored, and
412 * the result will be stored into.
413 * \return mantissa of the result
414 */
invFixp(FIXP_DBL op_m,int * op_e)415 static inline FIXP_DBL invFixp(FIXP_DBL op_m, int *op_e) {
416 if ((op_m == (FIXP_DBL)0x00000000) || (op_m == (FIXP_DBL)0x00000001)) {
417 *op_e = 31 - *op_e;
418 return ((LONG)MAXVAL_DBL);
419 }
420
421 INT tmp_exp;
422 FIXP_DBL tmp_inv = invSqrtNorm2(op_m, &tmp_exp);
423
424 *op_e = (tmp_exp << 1) - *op_e + 1;
425 return fPow2Div2(tmp_inv);
426 }
427 #endif /* FUNCTION_invFixp */
428
429 #ifndef FUNCTION_schur_div
430
431 /**
432 * \brief Divide two FIXP_DBL values with given precision.
433 * \param num dividend
434 * \param denum divisor
435 * \param count amount of significant bits of the result (starting to the MSB)
436 * \return num/divisor
437 */
438
439 FIXP_DBL schur_div(FIXP_DBL num, FIXP_DBL denum, INT count);
440
441 #endif /* FUNCTION_schur_div */
442
443 FIXP_DBL mul_dbl_sgl_rnd(const FIXP_DBL op1, const FIXP_SGL op2);
444
445 #ifndef FUNCTION_fMultNorm
446 /**
447 * \brief multiply two values with normalization, thus max precision.
448 * Author: Robert Weidner
449 *
450 * \param f1 first factor
451 * \param f2 second factor
452 * \param result_e pointer to an INT where the exponent of the result is stored
453 * into
454 * \return mantissa of the product f1*f2
455 */
456 FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2, INT *result_e);
457
458 /**
459 * \brief Multiply 2 values using maximum precision. The exponent of the result
460 * is 0.
461 * \param f1_m mantissa of factor 1
462 * \param f2_m mantissa of factor 2
463 * \return mantissa of the result with exponent equal to 0
464 */
fMultNorm(FIXP_DBL f1,FIXP_DBL f2)465 inline FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2) {
466 FIXP_DBL m;
467 INT e;
468
469 m = fMultNorm(f1, f2, &e);
470
471 m = scaleValueSaturate(m, e);
472
473 return m;
474 }
475
476 /**
477 * \brief Multiply 2 values with exponent and use given exponent for the
478 * mantissa of the result.
479 * \param f1_m mantissa of factor 1
480 * \param f1_e exponent of factor 1
481 * \param f2_m mantissa of factor 2
482 * \param f2_e exponent of factor 2
483 * \param result_e exponent for the returned mantissa of the result
484 * \return mantissa of the result with exponent equal to result_e
485 */
fMultNorm(FIXP_DBL f1_m,INT f1_e,FIXP_DBL f2_m,INT f2_e,INT result_e)486 inline FIXP_DBL fMultNorm(FIXP_DBL f1_m, INT f1_e, FIXP_DBL f2_m, INT f2_e,
487 INT result_e) {
488 FIXP_DBL m;
489 INT e;
490
491 m = fMultNorm(f1_m, f2_m, &e);
492
493 m = scaleValueSaturate(m, e + f1_e + f2_e - result_e);
494
495 return m;
496 }
497 #endif /* FUNCTION_fMultNorm */
498
499 #ifndef FUNCTION_fMultI
500 /**
501 * \brief Multiplies a fractional value and a integer value and performs
502 * rounding to nearest
503 * \param a fractional value
504 * \param b integer value
505 * \return integer value
506 */
fMultI(FIXP_DBL a,INT b)507 inline INT fMultI(FIXP_DBL a, INT b) {
508 FIXP_DBL m, mi;
509 INT m_e;
510
511 m = fMultNorm(a, (FIXP_DBL)b, &m_e);
512
513 if (m_e < (INT)0) {
514 if (m_e > (INT)-DFRACT_BITS) {
515 m = m >> ((-m_e) - 1);
516 mi = (m + (FIXP_DBL)1) >> 1;
517 } else {
518 mi = (FIXP_DBL)0;
519 }
520 } else {
521 mi = scaleValueSaturate(m, m_e);
522 }
523
524 return ((INT)mi);
525 }
526 #endif /* FUNCTION_fMultI */
527
528 #ifndef FUNCTION_fMultIfloor
529 /**
530 * \brief Multiplies a fractional value and a integer value and performs floor
531 * rounding
532 * \param a fractional value
533 * \param b integer value
534 * \return integer value
535 */
fMultIfloor(FIXP_DBL a,INT b)536 inline INT fMultIfloor(FIXP_DBL a, INT b) {
537 FIXP_DBL m, mi;
538 INT m_e;
539
540 m = fMultNorm(a, (FIXP_DBL)b, &m_e);
541
542 if (m_e < (INT)0) {
543 if (m_e > (INT)-DFRACT_BITS) {
544 mi = m >> (-m_e);
545 } else {
546 mi = (FIXP_DBL)0;
547 if (m < (FIXP_DBL)0) {
548 mi = (FIXP_DBL)-1;
549 }
550 }
551 } else {
552 mi = scaleValueSaturate(m, m_e);
553 }
554
555 return ((INT)mi);
556 }
557 #endif /* FUNCTION_fMultIfloor */
558
559 #ifndef FUNCTION_fMultIceil
560 /**
561 * \brief Multiplies a fractional value and a integer value and performs ceil
562 * rounding
563 * \param a fractional value
564 * \param b integer value
565 * \return integer value
566 */
fMultIceil(FIXP_DBL a,INT b)567 inline INT fMultIceil(FIXP_DBL a, INT b) {
568 FIXP_DBL m, mi;
569 INT m_e;
570
571 m = fMultNorm(a, (FIXP_DBL)b, &m_e);
572
573 if (m_e < (INT)0) {
574 if (m_e > (INT) - (DFRACT_BITS - 1)) {
575 mi = (m >> (-m_e));
576 if ((LONG)m & ((1 << (-m_e)) - 1)) {
577 mi = mi + (FIXP_DBL)1;
578 }
579 } else {
580 if (m > (FIXP_DBL)0) {
581 mi = (FIXP_DBL)1;
582 } else {
583 if ((m_e == -(DFRACT_BITS - 1)) && (m == (FIXP_DBL)MINVAL_DBL)) {
584 mi = (FIXP_DBL)-1;
585 } else {
586 mi = (FIXP_DBL)0;
587 }
588 }
589 }
590 } else {
591 mi = scaleValueSaturate(m, m_e);
592 }
593
594 return ((INT)mi);
595 }
596 #endif /* FUNCTION_fMultIceil */
597
598 #ifndef FUNCTION_fDivNorm
599 /**
600 * \brief Divide 2 FIXP_DBL values with normalization of input values.
601 * \param num numerator
602 * \param denum denominator
603 * \param result_e pointer to an INT where the exponent of the result is stored
604 * into
605 * \return num/denum with exponent = *result_e
606 */
607 FIXP_DBL fDivNorm(FIXP_DBL num, FIXP_DBL denom, INT *result_e);
608
609 /**
610 * \brief Divide 2 positive FIXP_DBL values with normalization of input values.
611 * \param num numerator
612 * \param denum denominator
613 * \return num/denum with exponent = 0
614 */
615 FIXP_DBL fDivNorm(FIXP_DBL num, FIXP_DBL denom);
616
617 /**
618 * \brief Divide 2 signed FIXP_DBL values with normalization of input values.
619 * \param num numerator
620 * \param denum denominator
621 * \param result_e pointer to an INT where the exponent of the result is stored
622 * into
623 * \return num/denum with exponent = *result_e
624 */
625 FIXP_DBL fDivNormSigned(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e);
626
627 /**
628 * \brief Divide 2 signed FIXP_DBL values with normalization of input values.
629 * \param num numerator
630 * \param denum denominator
631 * \return num/denum with exponent = 0
632 */
633 FIXP_DBL fDivNormSigned(FIXP_DBL num, FIXP_DBL denom);
634 #endif /* FUNCTION_fDivNorm */
635
636 /**
637 * \brief Adjust mantissa to exponent -1
638 * \param a_m mantissa of value to be adjusted
639 * \param pA_e pointer to the exponen of a_m
640 * \return adjusted mantissa
641 */
fAdjust(FIXP_DBL a_m,INT * pA_e)642 inline FIXP_DBL fAdjust(FIXP_DBL a_m, INT *pA_e) {
643 INT shift;
644
645 shift = fNorm(a_m) - 1;
646 *pA_e -= shift;
647
648 return scaleValue(a_m, shift);
649 }
650
651 #ifndef FUNCTION_fAddNorm
652 /**
653 * \brief Add two values with normalization
654 * \param a_m mantissa of first summand
655 * \param a_e exponent of first summand
656 * \param a_m mantissa of second summand
657 * \param a_e exponent of second summand
658 * \param pResult_e pointer to where the exponent of the result will be stored
659 * to.
660 * \return mantissa of result
661 */
fAddNorm(FIXP_DBL a_m,INT a_e,FIXP_DBL b_m,INT b_e,INT * pResult_e)662 inline FIXP_DBL fAddNorm(FIXP_DBL a_m, INT a_e, FIXP_DBL b_m, INT b_e,
663 INT *pResult_e) {
664 INT result_e;
665 FIXP_DBL result_m;
666
667 /* If one of the summands is zero, return the other.
668 This is necessary for the summation of a very small number to zero */
669 if (a_m == (FIXP_DBL)0) {
670 *pResult_e = b_e;
671 return b_m;
672 }
673 if (b_m == (FIXP_DBL)0) {
674 *pResult_e = a_e;
675 return a_m;
676 }
677
678 a_m = fAdjust(a_m, &a_e);
679 b_m = fAdjust(b_m, &b_e);
680
681 if (a_e > b_e) {
682 result_m = a_m + (b_m >> fMin(a_e - b_e, DFRACT_BITS - 1));
683 result_e = a_e;
684 } else {
685 result_m = (a_m >> fMin(b_e - a_e, DFRACT_BITS - 1)) + b_m;
686 result_e = b_e;
687 }
688
689 *pResult_e = result_e;
690 return result_m;
691 }
692
fAddNorm(FIXP_DBL a_m,INT a_e,FIXP_DBL b_m,INT b_e,INT result_e)693 inline FIXP_DBL fAddNorm(FIXP_DBL a_m, INT a_e, FIXP_DBL b_m, INT b_e,
694 INT result_e) {
695 FIXP_DBL result_m;
696
697 a_m = scaleValue(a_m, a_e - result_e);
698 b_m = scaleValue(b_m, b_e - result_e);
699
700 result_m = a_m + b_m;
701
702 return result_m;
703 }
704 #endif /* FUNCTION_fAddNorm */
705
706 /**
707 * \brief Divide 2 FIXP_DBL values with normalization of input values.
708 * \param num numerator
709 * \param denum denomintator
710 * \return num/denum with exponent = 0
711 */
712 FIXP_DBL fDivNormHighPrec(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e);
713
714 #ifndef FUNCTION_fPow
715 /**
716 * \brief return 2 ^ (exp_m * 2^exp_e)
717 * \param exp_m mantissa of the exponent to 2.0f
718 * \param exp_e exponent of the exponent to 2.0f
719 * \param result_e pointer to a INT where the exponent of the result will be
720 * stored into
721 * \return mantissa of the result
722 */
723 FIXP_DBL f2Pow(const FIXP_DBL exp_m, const INT exp_e, INT *result_e);
724
725 /**
726 * \brief return 2 ^ (exp_m * 2^exp_e). This version returns only the mantissa
727 * with implicit exponent of zero.
728 * \param exp_m mantissa of the exponent to 2.0f
729 * \param exp_e exponent of the exponent to 2.0f
730 * \return mantissa of the result
731 */
732 FIXP_DBL f2Pow(const FIXP_DBL exp_m, const INT exp_e);
733
734 /**
735 * \brief return x ^ (exp_m * 2^exp_e), where log2(x) = baseLd_m * 2^(baseLd_e).
736 * This saves the need to compute log2() of constant values (when x is a
737 * constant).
738 * \param baseLd_m mantissa of log2() of x.
739 * \param baseLd_e exponent of log2() of x.
740 * \param exp_m mantissa of the exponent to 2.0f
741 * \param exp_e exponent of the exponent to 2.0f
742 * \param result_e pointer to a INT where the exponent of the result will be
743 * stored into
744 * \return mantissa of the result
745 */
746 FIXP_DBL fLdPow(FIXP_DBL baseLd_m, INT baseLd_e, FIXP_DBL exp_m, INT exp_e,
747 INT *result_e);
748
749 /**
750 * \brief return x ^ (exp_m * 2^exp_e), where log2(x) = baseLd_m * 2^(baseLd_e).
751 * This saves the need to compute log2() of constant values (when x is a
752 * constant). This version does not return an exponent, which is
753 * implicitly 0.
754 * \param baseLd_m mantissa of log2() of x.
755 * \param baseLd_e exponent of log2() of x.
756 * \param exp_m mantissa of the exponent to 2.0f
757 * \param exp_e exponent of the exponent to 2.0f
758 * \return mantissa of the result
759 */
760 FIXP_DBL fLdPow(FIXP_DBL baseLd_m, INT baseLd_e, FIXP_DBL exp_m, INT exp_e);
761
762 /**
763 * \brief return (base_m * 2^base_e) ^ (exp * 2^exp_e). Use fLdPow() instead
764 * whenever possible.
765 * \param base_m mantissa of the base.
766 * \param base_e exponent of the base.
767 * \param exp_m mantissa of power to be calculated of the base.
768 * \param exp_e exponent of power to be calculated of the base.
769 * \param result_e pointer to a INT where the exponent of the result will be
770 * stored into.
771 * \return mantissa of the result.
772 */
773 FIXP_DBL fPow(FIXP_DBL base_m, INT base_e, FIXP_DBL exp_m, INT exp_e,
774 INT *result_e);
775
776 /**
777 * \brief return (base_m * 2^base_e) ^ N
778 * \param base_m mantissa of the base. Must not be negative.
779 * \param base_e exponent of the base
780 * \param N power to be calculated of the base
781 * \param result_e pointer to a INT where the exponent of the result will be
782 * stored into
783 * \return mantissa of the result
784 */
785 FIXP_DBL fPowInt(FIXP_DBL base_m, INT base_e, INT N, INT *result_e);
786 #endif /* #ifndef FUNCTION_fPow */
787
788 #ifndef FUNCTION_fLog2
789 /**
790 * \brief Calculate log(argument)/log(2) (logarithm with base 2). deprecated.
791 * Use fLog2() instead.
792 * \param arg mantissa of the argument
793 * \param arg_e exponent of the argument
794 * \param result_e pointer to an INT to store the exponent of the result
795 * \return the mantissa of the result.
796 * \param
797 */
798 FIXP_DBL CalcLog2(FIXP_DBL arg, INT arg_e, INT *result_e);
799
800 /**
801 * \brief calculate logarithm of base 2 of x_m * 2^(x_e)
802 * \param x_m mantissa of the input value.
803 * \param x_e exponent of the input value.
804 * \param pointer to an INT where the exponent of the result is returned into.
805 * \return mantissa of the result.
806 */
fLog2(FIXP_DBL x_m,INT x_e,INT * result_e)807 FDK_INLINE FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e, INT *result_e) {
808 FIXP_DBL result_m;
809
810 /* Short cut for zero and negative numbers. */
811 if (x_m <= FL2FXCONST_DBL(0.0f)) {
812 *result_e = DFRACT_BITS - 1;
813 return FL2FXCONST_DBL(-1.0f);
814 }
815
816 /* Calculate log2() */
817 {
818 FIXP_DBL x2_m;
819
820 /* Move input value x_m * 2^x_e toward 1.0, where the taylor approximation
821 of the function log(1-x) centered at 0 is most accurate. */
822 {
823 INT b_norm;
824
825 b_norm = fNormz(x_m) - 1;
826 x2_m = x_m << b_norm;
827 x_e = x_e - b_norm;
828 }
829
830 /* map x from log(x) domain to log(1-x) domain. */
831 x2_m = -(x2_m + FL2FXCONST_DBL(-1.0));
832
833 /* Taylor polynomial approximation of ln(1-x) */
834 {
835 FIXP_DBL px2_m;
836 result_m = FL2FXCONST_DBL(0.0);
837 px2_m = x2_m;
838 for (int i = 0; i < LD_PRECISION; i++) {
839 result_m = fMultAddDiv2(result_m, ldCoeff[i], px2_m);
840 px2_m = fMult(px2_m, x2_m);
841 }
842 }
843 /* Multiply result with 1/ln(2) = 1.0 + 0.442695040888 (get log2(x) from
844 * ln(x) result). */
845 result_m =
846 fMultAddDiv2(result_m, result_m,
847 FL2FXCONST_DBL(2.0 * 0.4426950408889634073599246810019));
848
849 /* Add exponent part. log2(x_m * 2^x_e) = log2(x_m) + x_e */
850 if (x_e != 0) {
851 int enorm;
852
853 enorm = DFRACT_BITS - fNorm((FIXP_DBL)x_e);
854 /* The -1 in the right shift of result_m compensates the fMultDiv2() above
855 * in the taylor polynomial evaluation loop.*/
856 result_m = (result_m >> (enorm - 1)) +
857 ((FIXP_DBL)x_e << (DFRACT_BITS - 1 - enorm));
858
859 *result_e = enorm;
860 } else {
861 /* 1 compensates the fMultDiv2() above in the taylor polynomial evaluation
862 * loop.*/
863 *result_e = 1;
864 }
865 }
866
867 return result_m;
868 }
869
870 /**
871 * \brief calculate logarithm of base 2 of x_m * 2^(x_e)
872 * \param x_m mantissa of the input value.
873 * \param x_e exponent of the input value.
874 * \return mantissa of the result with implicit exponent of LD_DATA_SHIFT.
875 */
fLog2(FIXP_DBL x_m,INT x_e)876 FDK_INLINE FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e) {
877 if (x_m <= FL2FXCONST_DBL(0.0f)) {
878 x_m = FL2FXCONST_DBL(-1.0f);
879 } else {
880 INT result_e;
881 x_m = fLog2(x_m, x_e, &result_e);
882 x_m = scaleValue(x_m, result_e - LD_DATA_SHIFT);
883 }
884 return x_m;
885 }
886
887 #endif /* FUNCTION_fLog2 */
888
889 #ifndef FUNCTION_fAddSaturate
890 /**
891 * \brief Add with saturation of the result.
892 * \param a first summand
893 * \param b second summand
894 * \return saturated sum of a and b.
895 */
fAddSaturate(const FIXP_SGL a,const FIXP_SGL b)896 inline FIXP_SGL fAddSaturate(const FIXP_SGL a, const FIXP_SGL b) {
897 LONG sum;
898
899 sum = (LONG)(SHORT)a + (LONG)(SHORT)b;
900 sum = fMax(fMin((INT)sum, (INT)MAXVAL_SGL), (INT)MINVAL_SGL);
901 return (FIXP_SGL)(SHORT)sum;
902 }
903
904 /**
905 * \brief Add with saturation of the result.
906 * \param a first summand
907 * \param b second summand
908 * \return saturated sum of a and b.
909 */
fAddSaturate(const FIXP_DBL a,const FIXP_DBL b)910 inline FIXP_DBL fAddSaturate(const FIXP_DBL a, const FIXP_DBL b) {
911 LONG sum;
912
913 sum = (LONG)(a >> 1) + (LONG)(b >> 1);
914 sum = fMax(fMin((INT)sum, (INT)(MAXVAL_DBL >> 1)), (INT)(MINVAL_DBL >> 1));
915 return (FIXP_DBL)(LONG)(sum << 1);
916 }
917 #endif /* FUNCTION_fAddSaturate */
918
919 INT fixp_floorToInt(FIXP_DBL f_inp, INT sf);
920 FIXP_DBL fixp_floor(FIXP_DBL f_inp, INT sf);
921
922 INT fixp_ceilToInt(FIXP_DBL f_inp, INT sf);
923 FIXP_DBL fixp_ceil(FIXP_DBL f_inp, INT sf);
924
925 INT fixp_truncateToInt(FIXP_DBL f_inp, INT sf);
926 FIXP_DBL fixp_truncate(FIXP_DBL f_inp, INT sf);
927
928 INT fixp_roundToInt(FIXP_DBL f_inp, INT sf);
929 FIXP_DBL fixp_round(FIXP_DBL f_inp, INT sf);
930
931 /*****************************************************************************
932
933 array for 1/n, n=1..80
934
935 ****************************************************************************/
936
937 extern const FIXP_DBL invCount[80];
938
939 LNK_SECTION_INITCODE
InitInvInt(void)940 inline void InitInvInt(void) {}
941
942 /**
943 * \brief Calculate the value of 1/i where i is a integer value. It supports
944 * input values from 1 upto (80-1).
945 * \param intValue Integer input value.
946 * \param FIXP_DBL representation of 1/intValue
947 */
GetInvInt(int intValue)948 inline FIXP_DBL GetInvInt(int intValue) {
949 return invCount[fMin(fMax(intValue, 0), 80 - 1)];
950 }
951
952 #endif /* FIXPOINT_MATH_H */
953