1 /* -----------------------------------------------------------------------------
2 Software License for The Fraunhofer FDK AAC Codec Library for Android
3 
4 © Copyright  1995 - 2018 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   if (a_e > b_e) {
175     return ((b_m >> fMin(a_e - b_e, DFRACT_BITS - 1)) > a_m);
176   } else {
177     return ((a_m >> fMin(b_e - a_e, DFRACT_BITS - 1)) < b_m);
178   }
179 }
180 
fIsLessThan(FIXP_SGL a_m,INT a_e,FIXP_SGL b_m,INT b_e)181 FDK_INLINE INT fIsLessThan(FIXP_SGL a_m, INT a_e, FIXP_SGL b_m, INT b_e) {
182   if (a_e > b_e) {
183     return ((b_m >> fMin(a_e - b_e, FRACT_BITS - 1)) > a_m);
184   } else {
185     return ((a_m >> fMin(b_e - a_e, FRACT_BITS - 1)) < b_m);
186   }
187 }
188 #endif
189 
190 /**
191  * \brief deprecated. Use fLog2() instead.
192  */
193 #define CalcLdData(op) fLog2(op, 0)
194 
195 void LdDataVector(FIXP_DBL *srcVector, FIXP_DBL *destVector, INT number);
196 
197 extern const UINT exp2_tab_long[32];
198 extern const UINT exp2w_tab_long[32];
199 extern const UINT exp2x_tab_long[32];
200 
201 LNK_SECTION_CODE_L1
CalcInvLdData(const FIXP_DBL x)202 FDK_INLINE FIXP_DBL CalcInvLdData(const FIXP_DBL x) {
203   int set_zero = (x < FL2FXCONST_DBL(-31.0 / 64.0)) ? 0 : 1;
204   int set_max = (x >= FL2FXCONST_DBL(31.0 / 64.0)) | (x == FL2FXCONST_DBL(0.0));
205 
206   FIXP_SGL frac = (FIXP_SGL)((LONG)x & 0x3FF);
207   UINT index3 = (UINT)(LONG)(x >> 10) & 0x1F;
208   UINT index2 = (UINT)(LONG)(x >> 15) & 0x1F;
209   UINT index1 = (UINT)(LONG)(x >> 20) & 0x1F;
210   int exp = fMin(31, ((x > FL2FXCONST_DBL(0.0f)) ? (31 - (int)(x >> 25))
211                                                  : (int)(-(x >> 25))));
212 
213   UINT lookup1 = exp2_tab_long[index1] * set_zero;
214   UINT lookup2 = exp2w_tab_long[index2];
215   UINT lookup3 = exp2x_tab_long[index3];
216   UINT lookup3f =
217       lookup3 + (UINT)(LONG)fMultDiv2((FIXP_DBL)(0x0016302F), (FIXP_SGL)frac);
218 
219   UINT lookup12 = (UINT)(LONG)fMult((FIXP_DBL)lookup1, (FIXP_DBL)lookup2);
220   UINT lookup = (UINT)(LONG)fMult((FIXP_DBL)lookup12, (FIXP_DBL)lookup3f);
221 
222   FIXP_DBL retVal = (lookup << 3) >> exp;
223 
224   if (set_max) {
225     retVal = (FIXP_DBL)MAXVAL_DBL;
226   }
227 
228   return retVal;
229 }
230 
231 void InitLdInt();
232 FIXP_DBL CalcLdInt(INT i);
233 
234 extern const USHORT sqrt_tab[49];
235 
sqrtFixp_lookup(FIXP_DBL x)236 inline FIXP_DBL sqrtFixp_lookup(FIXP_DBL x) {
237   UINT y = (INT)x;
238   UCHAR is_zero = (y == 0);
239   INT zeros = fixnormz_D(y) & 0x1e;
240   y <<= zeros;
241   UINT idx = (y >> 26) - 16;
242   USHORT frac = (y >> 10) & 0xffff;
243   USHORT nfrac = 0xffff ^ frac;
244   UINT t = (UINT)nfrac * sqrt_tab[idx] + (UINT)frac * sqrt_tab[idx + 1];
245   t = t >> (zeros >> 1);
246   return (is_zero ? 0 : t);
247 }
248 
sqrtFixp_lookup(FIXP_DBL x,INT * x_e)249 inline FIXP_DBL sqrtFixp_lookup(FIXP_DBL x, INT *x_e) {
250   UINT y = (INT)x;
251   INT e;
252 
253   if (x == (FIXP_DBL)0) {
254     return x;
255   }
256 
257   /* Normalize */
258   e = fixnormz_D(y);
259   y <<= e;
260   e = *x_e - e + 2;
261 
262   /* Correct odd exponent. */
263   if (e & 1) {
264     y >>= 1;
265     e++;
266   }
267   /* Get square root */
268   UINT idx = (y >> 26) - 16;
269   USHORT frac = (y >> 10) & 0xffff;
270   USHORT nfrac = 0xffff ^ frac;
271   UINT t = (UINT)nfrac * sqrt_tab[idx] + (UINT)frac * sqrt_tab[idx + 1];
272 
273   /* Write back exponent */
274   *x_e = e >> 1;
275   return (FIXP_DBL)(LONG)(t >> 1);
276 }
277 
278 void InitInvSqrtTab();
279 
280 #ifndef FUNCTION_invSqrtNorm2
281 /**
282  * \brief calculate 1.0/sqrt(op)
283  * \param op_m mantissa of input value.
284  * \param result_e pointer to return the exponent of the result
285  * \return mantissa of the result
286  */
287 /*****************************************************************************
288   delivers 1/sqrt(op) normalized to .5...1 and the shift value of the OUTPUT,
289   i.e. the denormalized result is 1/sqrt(op) = invSqrtNorm(op) * 2^(shift)
290   uses Newton-iteration for approximation
291       Q(n+1) = Q(n) + Q(n) * (0.5 - 2 * V * Q(n)^2)
292       with Q = 0.5* V ^-0.5; 0.5 <= V < 1.0
293 *****************************************************************************/
invSqrtNorm2(FIXP_DBL op,INT * shift)294 static FDK_FORCEINLINE FIXP_DBL invSqrtNorm2(FIXP_DBL op, INT *shift) {
295   FIXP_DBL val = op;
296   FIXP_DBL reg1, reg2;
297 
298   if (val == FL2FXCONST_DBL(0.0)) {
299     *shift = 16;
300     return ((LONG)MAXVAL_DBL); /* maximum positive value */
301   }
302 
303 #define INVSQRTNORM2_LINEAR_INTERPOLATE
304 #define INVSQRTNORM2_LINEAR_INTERPOLATE_HQ
305 
306   /* normalize input, calculate shift value */
307   FDK_ASSERT(val > FL2FXCONST_DBL(0.0));
308   *shift = fNormz(val) - 1; /* CountLeadingBits() is not necessary here since
309                                test value is always > 0 */
310   val <<= *shift;           /* normalized input V */
311   *shift += 2;              /* bias for exponent */
312 
313 #if defined(INVSQRTNORM2_LINEAR_INTERPOLATE)
314   INT index =
315       (INT)(val >> (DFRACT_BITS - 1 - (SQRT_BITS + 1))) & SQRT_BITS_MASK;
316   FIXP_DBL Fract =
317       (FIXP_DBL)(((INT)val & SQRT_FRACT_BITS_MASK) << (SQRT_BITS + 1));
318   FIXP_DBL diff = invSqrtTab[index + 1] - invSqrtTab[index];
319   reg1 = invSqrtTab[index] + (fMultDiv2(diff, Fract) << 1);
320 #if defined(INVSQRTNORM2_LINEAR_INTERPOLATE_HQ)
321   /* reg1 = t[i] + (t[i+1]-t[i])*fract ... already computed ...
322                                        + (1-fract)fract*(t[i+2]-t[i+1])/2 */
323   if (Fract != (FIXP_DBL)0) {
324     /* fract = fract * (1 - fract) */
325     Fract = fMultDiv2(Fract, (FIXP_DBL)((ULONG)0x80000000 - (ULONG)Fract)) << 1;
326     diff = diff - (invSqrtTab[index + 2] - invSqrtTab[index + 1]);
327     reg1 = fMultAddDiv2(reg1, Fract, diff);
328   }
329 #endif /* INVSQRTNORM2_LINEAR_INTERPOLATE_HQ */
330 #else
331 #error \
332     "Either define INVSQRTNORM2_NEWTON_ITERATE or INVSQRTNORM2_LINEAR_INTERPOLATE"
333 #endif
334   /* calculate the output exponent = input exp/2 */
335   if (*shift & 0x00000001) { /* odd shift values ? */
336     /* Note: Do not use rounded value 0x5A82799A to avoid overflow with
337      * shift-by-2 */
338     reg2 = (FIXP_DBL)0x5A827999;
339     /* FL2FXCONST_DBL(0.707106781186547524400844362104849f);*/ /* 1/sqrt(2);
340                                                                 */
341     reg1 = fMultDiv2(reg1, reg2) << 2;
342   }
343 
344   *shift = *shift >> 1;
345 
346   return (reg1);
347 }
348 #endif /* FUNCTION_invSqrtNorm2 */
349 
350 #ifndef FUNCTION_sqrtFixp
sqrtFixp(FIXP_DBL op)351 static FDK_FORCEINLINE FIXP_DBL sqrtFixp(FIXP_DBL op) {
352   INT tmp_exp = 0;
353   FIXP_DBL tmp_inv = invSqrtNorm2(op, &tmp_exp);
354 
355   FDK_ASSERT(tmp_exp > 0);
356   return ((FIXP_DBL)(fMultDiv2((op << (tmp_exp - 1)), tmp_inv) << 2));
357 }
358 #endif /* FUNCTION_sqrtFixp */
359 
360 #ifndef FUNCTION_invFixp
361 /**
362  * \brief calculate 1.0/op
363  * \param op mantissa of the input value.
364  * \return mantissa of the result with implicit exponent of 31
365  * \exceptions are provided for op=0,1 setting max. positive value
366  */
invFixp(FIXP_DBL op)367 static inline FIXP_DBL invFixp(FIXP_DBL op) {
368   if ((op == (FIXP_DBL)0x00000000) || (op == (FIXP_DBL)0x00000001)) {
369     return ((LONG)MAXVAL_DBL);
370   }
371   INT tmp_exp;
372   FIXP_DBL tmp_inv = invSqrtNorm2(op, &tmp_exp);
373   FDK_ASSERT((31 - (2 * tmp_exp + 1)) >= 0);
374   int shift = 31 - (2 * tmp_exp + 1);
375   tmp_inv = fPow2Div2(tmp_inv);
376   if (shift) {
377     tmp_inv = ((tmp_inv >> (shift - 1)) + (FIXP_DBL)1) >> 1;
378   }
379   return tmp_inv;
380 }
381 
382 /**
383  * \brief calculate 1.0/(op_m * 2^op_e)
384  * \param op_m mantissa of the input value.
385  * \param op_e pointer into were the exponent of the input value is stored, and
386  * the result will be stored into.
387  * \return mantissa of the result
388  */
invFixp(FIXP_DBL op_m,int * op_e)389 static inline FIXP_DBL invFixp(FIXP_DBL op_m, int *op_e) {
390   if ((op_m == (FIXP_DBL)0x00000000) || (op_m == (FIXP_DBL)0x00000001)) {
391     *op_e = 31 - *op_e;
392     return ((LONG)MAXVAL_DBL);
393   }
394 
395   INT tmp_exp;
396   FIXP_DBL tmp_inv = invSqrtNorm2(op_m, &tmp_exp);
397 
398   *op_e = (tmp_exp << 1) - *op_e + 1;
399   return fPow2Div2(tmp_inv);
400 }
401 #endif /* FUNCTION_invFixp */
402 
403 #ifndef FUNCTION_schur_div
404 
405 /**
406  * \brief Divide two FIXP_DBL values with given precision.
407  * \param num dividend
408  * \param denum divisor
409  * \param count amount of significant bits of the result (starting to the MSB)
410  * \return num/divisor
411  */
412 
413 FIXP_DBL schur_div(FIXP_DBL num, FIXP_DBL denum, INT count);
414 
415 #endif /* FUNCTION_schur_div */
416 
417 FIXP_DBL mul_dbl_sgl_rnd(const FIXP_DBL op1, const FIXP_SGL op2);
418 
419 #ifndef FUNCTION_fMultNorm
420 /**
421  * \brief multiply two values with normalization, thus max precision.
422  * Author: Robert Weidner
423  *
424  * \param f1 first factor
425  * \param f2 second factor
426  * \param result_e pointer to an INT where the exponent of the result is stored
427  * into
428  * \return mantissa of the product f1*f2
429  */
430 FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2, INT *result_e);
431 
432 /**
433  * \brief Multiply 2 values using maximum precision. The exponent of the result
434  * is 0.
435  * \param f1_m mantissa of factor 1
436  * \param f2_m mantissa of factor 2
437  * \return mantissa of the result with exponent equal to 0
438  */
fMultNorm(FIXP_DBL f1,FIXP_DBL f2)439 inline FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2) {
440   FIXP_DBL m;
441   INT e;
442 
443   m = fMultNorm(f1, f2, &e);
444 
445   m = scaleValueSaturate(m, e);
446 
447   return m;
448 }
449 
450 /**
451  * \brief Multiply 2 values with exponent and use given exponent for the
452  * mantissa of the result.
453  * \param f1_m mantissa of factor 1
454  * \param f1_e exponent of factor 1
455  * \param f2_m mantissa of factor 2
456  * \param f2_e exponent of factor 2
457  * \param result_e exponent for the returned mantissa of the result
458  * \return mantissa of the result with exponent equal to result_e
459  */
fMultNorm(FIXP_DBL f1_m,INT f1_e,FIXP_DBL f2_m,INT f2_e,INT result_e)460 inline FIXP_DBL fMultNorm(FIXP_DBL f1_m, INT f1_e, FIXP_DBL f2_m, INT f2_e,
461                           INT result_e) {
462   FIXP_DBL m;
463   INT e;
464 
465   m = fMultNorm(f1_m, f2_m, &e);
466 
467   m = scaleValueSaturate(m, e + f1_e + f2_e - result_e);
468 
469   return m;
470 }
471 #endif /* FUNCTION_fMultNorm */
472 
473 #ifndef FUNCTION_fMultI
474 /**
475  * \brief Multiplies a fractional value and a integer value and performs
476  * rounding to nearest
477  * \param a fractional value
478  * \param b integer value
479  * \return integer value
480  */
fMultI(FIXP_DBL a,INT b)481 inline INT fMultI(FIXP_DBL a, INT b) {
482   FIXP_DBL m, mi;
483   INT m_e;
484 
485   m = fMultNorm(a, (FIXP_DBL)b, &m_e);
486 
487   if (m_e < (INT)0) {
488     if (m_e > (INT)-DFRACT_BITS) {
489       m = m >> ((-m_e) - 1);
490       mi = (m + (FIXP_DBL)1) >> 1;
491     } else {
492       mi = (FIXP_DBL)0;
493     }
494   } else {
495     mi = scaleValueSaturate(m, m_e);
496   }
497 
498   return ((INT)mi);
499 }
500 #endif /* FUNCTION_fMultI */
501 
502 #ifndef FUNCTION_fMultIfloor
503 /**
504  * \brief Multiplies a fractional value and a integer value and performs floor
505  * rounding
506  * \param a fractional value
507  * \param b integer value
508  * \return integer value
509  */
fMultIfloor(FIXP_DBL a,INT b)510 inline INT fMultIfloor(FIXP_DBL a, INT b) {
511   FIXP_DBL m, mi;
512   INT m_e;
513 
514   m = fMultNorm(a, (FIXP_DBL)b, &m_e);
515 
516   if (m_e < (INT)0) {
517     if (m_e > (INT)-DFRACT_BITS) {
518       mi = m >> (-m_e);
519     } else {
520       mi = (FIXP_DBL)0;
521       if (m < (FIXP_DBL)0) {
522         mi = (FIXP_DBL)-1;
523       }
524     }
525   } else {
526     mi = scaleValueSaturate(m, m_e);
527   }
528 
529   return ((INT)mi);
530 }
531 #endif /* FUNCTION_fMultIfloor */
532 
533 #ifndef FUNCTION_fMultIceil
534 /**
535  * \brief Multiplies a fractional value and a integer value and performs ceil
536  * rounding
537  * \param a fractional value
538  * \param b integer value
539  * \return integer value
540  */
fMultIceil(FIXP_DBL a,INT b)541 inline INT fMultIceil(FIXP_DBL a, INT b) {
542   FIXP_DBL m, mi;
543   INT m_e;
544 
545   m = fMultNorm(a, (FIXP_DBL)b, &m_e);
546 
547   if (m_e < (INT)0) {
548     if (m_e > (INT)-DFRACT_BITS) {
549       mi = (m >> (-m_e));
550       if ((LONG)m & ((1 << (-m_e)) - 1)) {
551         mi = mi + (FIXP_DBL)1;
552       }
553     } else {
554       mi = (FIXP_DBL)1;
555       if (m < (FIXP_DBL)0) {
556         mi = (FIXP_DBL)0;
557       }
558     }
559   } else {
560     mi = scaleValueSaturate(m, m_e);
561   }
562 
563   return ((INT)mi);
564 }
565 #endif /* FUNCTION_fMultIceil */
566 
567 #ifndef FUNCTION_fDivNorm
568 /**
569  * \brief Divide 2 FIXP_DBL values with normalization of input values.
570  * \param num numerator
571  * \param denum denominator
572  * \param result_e pointer to an INT where the exponent of the result is stored
573  * into
574  * \return num/denum with exponent = *result_e
575  */
576 FIXP_DBL fDivNorm(FIXP_DBL num, FIXP_DBL denom, INT *result_e);
577 
578 /**
579  * \brief Divide 2 positive FIXP_DBL values with normalization of input values.
580  * \param num numerator
581  * \param denum denominator
582  * \return num/denum with exponent = 0
583  */
584 FIXP_DBL fDivNorm(FIXP_DBL num, FIXP_DBL denom);
585 
586 /**
587  * \brief Divide 2 signed FIXP_DBL values with normalization of input values.
588  * \param num numerator
589  * \param denum denominator
590  * \param result_e pointer to an INT where the exponent of the result is stored
591  * into
592  * \return num/denum with exponent = *result_e
593  */
594 FIXP_DBL fDivNormSigned(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e);
595 
596 /**
597  * \brief Divide 2 signed FIXP_DBL values with normalization of input values.
598  * \param num numerator
599  * \param denum denominator
600  * \return num/denum with exponent = 0
601  */
602 FIXP_DBL fDivNormSigned(FIXP_DBL num, FIXP_DBL denom);
603 #endif /* FUNCTION_fDivNorm */
604 
605 /**
606  * \brief Adjust mantissa to exponent -1
607  * \param a_m mantissa of value to be adjusted
608  * \param pA_e pointer to the exponen of a_m
609  * \return adjusted mantissa
610  */
fAdjust(FIXP_DBL a_m,INT * pA_e)611 inline FIXP_DBL fAdjust(FIXP_DBL a_m, INT *pA_e) {
612   INT shift;
613 
614   shift = fNorm(a_m) - 1;
615   *pA_e -= shift;
616 
617   return scaleValue(a_m, shift);
618 }
619 
620 #ifndef FUNCTION_fAddNorm
621 /**
622  * \brief Add two values with normalization
623  * \param a_m mantissa of first summand
624  * \param a_e exponent of first summand
625  * \param a_m mantissa of second summand
626  * \param a_e exponent of second summand
627  * \param pResult_e pointer to where the exponent of the result will be stored
628  * to.
629  * \return mantissa of result
630  */
fAddNorm(FIXP_DBL a_m,INT a_e,FIXP_DBL b_m,INT b_e,INT * pResult_e)631 inline FIXP_DBL fAddNorm(FIXP_DBL a_m, INT a_e, FIXP_DBL b_m, INT b_e,
632                          INT *pResult_e) {
633   INT result_e;
634   FIXP_DBL result_m;
635 
636   /* If one of the summands is zero, return the other.
637      This is necessary for the summation of a very small number to zero */
638   if (a_m == (FIXP_DBL)0) {
639     *pResult_e = b_e;
640     return b_m;
641   }
642   if (b_m == (FIXP_DBL)0) {
643     *pResult_e = a_e;
644     return a_m;
645   }
646 
647   a_m = fAdjust(a_m, &a_e);
648   b_m = fAdjust(b_m, &b_e);
649 
650   if (a_e > b_e) {
651     result_m = a_m + (b_m >> fMin(a_e - b_e, DFRACT_BITS - 1));
652     result_e = a_e;
653   } else {
654     result_m = (a_m >> fMin(b_e - a_e, DFRACT_BITS - 1)) + b_m;
655     result_e = b_e;
656   }
657 
658   *pResult_e = result_e;
659   return result_m;
660 }
661 
fAddNorm(FIXP_DBL a_m,INT a_e,FIXP_DBL b_m,INT b_e,INT result_e)662 inline FIXP_DBL fAddNorm(FIXP_DBL a_m, INT a_e, FIXP_DBL b_m, INT b_e,
663                          INT result_e) {
664   FIXP_DBL result_m;
665 
666   a_m = scaleValue(a_m, a_e - result_e);
667   b_m = scaleValue(b_m, b_e - result_e);
668 
669   result_m = a_m + b_m;
670 
671   return result_m;
672 }
673 #endif /* FUNCTION_fAddNorm */
674 
675 /**
676  * \brief Divide 2 FIXP_DBL values with normalization of input values.
677  * \param num numerator
678  * \param denum denomintator
679  * \return num/denum with exponent = 0
680  */
681 FIXP_DBL fDivNormHighPrec(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e);
682 
683 #ifndef FUNCTION_fPow
684 /**
685  * \brief return 2 ^ (exp_m * 2^exp_e)
686  * \param exp_m mantissa of the exponent to 2.0f
687  * \param exp_e exponent of the exponent to 2.0f
688  * \param result_e pointer to a INT where the exponent of the result will be
689  * stored into
690  * \return mantissa of the result
691  */
692 FIXP_DBL f2Pow(const FIXP_DBL exp_m, const INT exp_e, INT *result_e);
693 
694 /**
695  * \brief return 2 ^ (exp_m * 2^exp_e). This version returns only the mantissa
696  * with implicit exponent of zero.
697  * \param exp_m mantissa of the exponent to 2.0f
698  * \param exp_e exponent of the exponent to 2.0f
699  * \return mantissa of the result
700  */
701 FIXP_DBL f2Pow(const FIXP_DBL exp_m, const INT exp_e);
702 
703 /**
704  * \brief return x ^ (exp_m * 2^exp_e), where log2(x) = baseLd_m * 2^(baseLd_e).
705  * This saves the need to compute log2() of constant values (when x is a
706  * constant).
707  * \param baseLd_m mantissa of log2() of x.
708  * \param baseLd_e exponent of log2() of x.
709  * \param exp_m mantissa of the exponent to 2.0f
710  * \param exp_e exponent of the exponent to 2.0f
711  * \param result_e pointer to a INT where the exponent of the result will be
712  * stored into
713  * \return mantissa of the result
714  */
715 FIXP_DBL fLdPow(FIXP_DBL baseLd_m, INT baseLd_e, FIXP_DBL exp_m, INT exp_e,
716                 INT *result_e);
717 
718 /**
719  * \brief return x ^ (exp_m * 2^exp_e), where log2(x) = baseLd_m * 2^(baseLd_e).
720  * This saves the need to compute log2() of constant values (when x is a
721  * constant). This version does not return an exponent, which is
722  * implicitly 0.
723  * \param baseLd_m mantissa of log2() of x.
724  * \param baseLd_e exponent of log2() of x.
725  * \param exp_m mantissa of the exponent to 2.0f
726  * \param exp_e exponent of the exponent to 2.0f
727  * \return mantissa of the result
728  */
729 FIXP_DBL fLdPow(FIXP_DBL baseLd_m, INT baseLd_e, FIXP_DBL exp_m, INT exp_e);
730 
731 /**
732  * \brief return (base_m * 2^base_e) ^ (exp * 2^exp_e). Use fLdPow() instead
733  * whenever possible.
734  * \param base_m mantissa of the base.
735  * \param base_e exponent of the base.
736  * \param exp_m mantissa of power to be calculated of the base.
737  * \param exp_e exponent of power to be calculated of the base.
738  * \param result_e pointer to a INT where the exponent of the result will be
739  * stored into.
740  * \return mantissa of the result.
741  */
742 FIXP_DBL fPow(FIXP_DBL base_m, INT base_e, FIXP_DBL exp_m, INT exp_e,
743               INT *result_e);
744 
745 /**
746  * \brief return (base_m * 2^base_e) ^ N
747  * \param base_m mantissa of the base
748  * \param base_e exponent of the base
749  * \param N power to be calculated of the base
750  * \param result_e pointer to a INT where the exponent of the result will be
751  * stored into
752  * \return mantissa of the result
753  */
754 FIXP_DBL fPowInt(FIXP_DBL base_m, INT base_e, INT N, INT *result_e);
755 #endif /* #ifndef FUNCTION_fPow */
756 
757 #ifndef FUNCTION_fLog2
758 /**
759  * \brief Calculate log(argument)/log(2) (logarithm with base 2). deprecated.
760  * Use fLog2() instead.
761  * \param arg mantissa of the argument
762  * \param arg_e exponent of the argument
763  * \param result_e pointer to an INT to store the exponent of the result
764  * \return the mantissa of the result.
765  * \param
766  */
767 FIXP_DBL CalcLog2(FIXP_DBL arg, INT arg_e, INT *result_e);
768 
769 /**
770  * \brief calculate logarithm of base 2 of x_m * 2^(x_e)
771  * \param x_m mantissa of the input value.
772  * \param x_e exponent of the input value.
773  * \param pointer to an INT where the exponent of the result is returned into.
774  * \return mantissa of the result.
775  */
fLog2(FIXP_DBL x_m,INT x_e,INT * result_e)776 FDK_INLINE FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e, INT *result_e) {
777   FIXP_DBL result_m;
778 
779   /* Short cut for zero and negative numbers. */
780   if (x_m <= FL2FXCONST_DBL(0.0f)) {
781     *result_e = DFRACT_BITS - 1;
782     return FL2FXCONST_DBL(-1.0f);
783   }
784 
785   /* Calculate log2() */
786   {
787     FIXP_DBL x2_m;
788 
789     /* Move input value x_m * 2^x_e toward 1.0, where the taylor approximation
790        of the function log(1-x) centered at 0 is most accurate. */
791     {
792       INT b_norm;
793 
794       b_norm = fNormz(x_m) - 1;
795       x2_m = x_m << b_norm;
796       x_e = x_e - b_norm;
797     }
798 
799     /* map x from log(x) domain to log(1-x) domain. */
800     x2_m = -(x2_m + FL2FXCONST_DBL(-1.0));
801 
802     /* Taylor polynomial approximation of ln(1-x) */
803     {
804       FIXP_DBL px2_m;
805       result_m = FL2FXCONST_DBL(0.0);
806       px2_m = x2_m;
807       for (int i = 0; i < LD_PRECISION; i++) {
808         result_m = fMultAddDiv2(result_m, ldCoeff[i], px2_m);
809         px2_m = fMult(px2_m, x2_m);
810       }
811     }
812     /* Multiply result with 1/ln(2) = 1.0 + 0.442695040888 (get log2(x) from
813      * ln(x) result). */
814     result_m =
815         fMultAddDiv2(result_m, result_m,
816                      FL2FXCONST_DBL(2.0 * 0.4426950408889634073599246810019));
817 
818     /* Add exponent part. log2(x_m * 2^x_e) = log2(x_m) + x_e */
819     if (x_e != 0) {
820       int enorm;
821 
822       enorm = DFRACT_BITS - fNorm((FIXP_DBL)x_e);
823       /* The -1 in the right shift of result_m compensates the fMultDiv2() above
824        * in the taylor polynomial evaluation loop.*/
825       result_m = (result_m >> (enorm - 1)) +
826                  ((FIXP_DBL)x_e << (DFRACT_BITS - 1 - enorm));
827 
828       *result_e = enorm;
829     } else {
830       /* 1 compensates the fMultDiv2() above in the taylor polynomial evaluation
831        * loop.*/
832       *result_e = 1;
833     }
834   }
835 
836   return result_m;
837 }
838 
839 /**
840  * \brief calculate logarithm of base 2 of x_m * 2^(x_e)
841  * \param x_m mantissa of the input value.
842  * \param x_e exponent of the input value.
843  * \return mantissa of the result with implicit exponent of LD_DATA_SHIFT.
844  */
fLog2(FIXP_DBL x_m,INT x_e)845 FDK_INLINE FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e) {
846   if (x_m <= FL2FXCONST_DBL(0.0f)) {
847     x_m = FL2FXCONST_DBL(-1.0f);
848   } else {
849     INT result_e;
850     x_m = fLog2(x_m, x_e, &result_e);
851     x_m = scaleValue(x_m, result_e - LD_DATA_SHIFT);
852   }
853   return x_m;
854 }
855 
856 #endif /* FUNCTION_fLog2 */
857 
858 #ifndef FUNCTION_fAddSaturate
859 /**
860  * \brief Add with saturation of the result.
861  * \param a first summand
862  * \param b second summand
863  * \return saturated sum of a and b.
864  */
fAddSaturate(const FIXP_SGL a,const FIXP_SGL b)865 inline FIXP_SGL fAddSaturate(const FIXP_SGL a, const FIXP_SGL b) {
866   LONG sum;
867 
868   sum = (LONG)(SHORT)a + (LONG)(SHORT)b;
869   sum = fMax(fMin((INT)sum, (INT)MAXVAL_SGL), (INT)MINVAL_SGL);
870   return (FIXP_SGL)(SHORT)sum;
871 }
872 
873 /**
874  * \brief Add with saturation of the result.
875  * \param a first summand
876  * \param b second summand
877  * \return saturated sum of a and b.
878  */
fAddSaturate(const FIXP_DBL a,const FIXP_DBL b)879 inline FIXP_DBL fAddSaturate(const FIXP_DBL a, const FIXP_DBL b) {
880   LONG sum;
881 
882   sum = (LONG)(a >> 1) + (LONG)(b >> 1);
883   sum = fMax(fMin((INT)sum, (INT)(MAXVAL_DBL >> 1)), (INT)(MINVAL_DBL >> 1));
884   return (FIXP_DBL)(LONG)(sum << 1);
885 }
886 #endif /* FUNCTION_fAddSaturate */
887 
888 INT fixp_floorToInt(FIXP_DBL f_inp, INT sf);
889 FIXP_DBL fixp_floor(FIXP_DBL f_inp, INT sf);
890 
891 INT fixp_ceilToInt(FIXP_DBL f_inp, INT sf);
892 FIXP_DBL fixp_ceil(FIXP_DBL f_inp, INT sf);
893 
894 INT fixp_truncateToInt(FIXP_DBL f_inp, INT sf);
895 FIXP_DBL fixp_truncate(FIXP_DBL f_inp, INT sf);
896 
897 INT fixp_roundToInt(FIXP_DBL f_inp, INT sf);
898 FIXP_DBL fixp_round(FIXP_DBL f_inp, INT sf);
899 
900 /*****************************************************************************
901 
902  array for 1/n, n=1..80
903 
904 ****************************************************************************/
905 
906 extern const FIXP_DBL invCount[80];
907 
908 LNK_SECTION_INITCODE
InitInvInt(void)909 inline void InitInvInt(void) {}
910 
911 /**
912  * \brief Calculate the value of 1/i where i is a integer value. It supports
913  *        input values from 1 upto (80-1).
914  * \param intValue Integer input value.
915  * \param FIXP_DBL representation of 1/intValue
916  */
GetInvInt(int intValue)917 inline FIXP_DBL GetInvInt(int intValue) {
918   return invCount[fMin(fMax(intValue, 0), 80 - 1)];
919 }
920 
921 #endif /* FIXPOINT_MATH_H */
922