1 /* -----------------------------------------------------------------------------
2 Software License for The Fraunhofer FDK AAC Codec Library for Android
3 
4 © Copyright  1995 - 2020 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 #include "fixpoint_math.h"
104 
105 /*
106  * Hardware specific implementations
107  */
108 
109 /*
110  * Fallback implementations
111  */
112 
113 /*****************************************************************************
114   functionname: LdDataVector
115 *****************************************************************************/
116 LNK_SECTION_CODE_L1
LdDataVector(FIXP_DBL * srcVector,FIXP_DBL * destVector,INT n)117 void LdDataVector(FIXP_DBL *srcVector, FIXP_DBL *destVector, INT n) {
118   INT i;
119   for (i = 0; i < n; i++) {
120     destVector[i] = fLog2(srcVector[i], 0);
121   }
122 }
123 
124 #define MAX_POW2_PRECISION 8
125 #ifndef SINETABLE_16BIT
126 #define POW2_PRECISION MAX_POW2_PRECISION
127 #else
128 #define POW2_PRECISION 5
129 #endif
130 
131 /*
132   Taylor series coefficients of the function x^2. The first coefficient is
133   ommited (equal to 1.0).
134 
135   pow2Coeff[i-1] = (1/i!) d^i(2^x)/dx^i, i=1..MAX_POW2_PRECISION
136   To evaluate the taylor series around x = 0, the coefficients are: 1/!i *
137   ln(2)^i
138  */
139 #ifndef POW2COEFF_16BIT
140 RAM_ALIGN
141 LNK_SECTION_CONSTDATA_L1
142 static const FIXP_DBL pow2Coeff[MAX_POW2_PRECISION] = {
143     FL2FXCONST_DBL(0.693147180559945309417232121458177),   /* ln(2)^1 /1! */
144     FL2FXCONST_DBL(0.240226506959100712333551263163332),   /* ln(2)^2 /2! */
145     FL2FXCONST_DBL(0.0555041086648215799531422637686218),  /* ln(2)^3 /3! */
146     FL2FXCONST_DBL(0.00961812910762847716197907157365887), /* ln(2)^4 /4! */
147     FL2FXCONST_DBL(0.00133335581464284434234122219879962), /* ln(2)^5 /5! */
148     FL2FXCONST_DBL(1.54035303933816099544370973327423e-4), /* ln(2)^6 /6! */
149     FL2FXCONST_DBL(1.52527338040598402800254390120096e-5), /* ln(2)^7 /7! */
150     FL2FXCONST_DBL(1.32154867901443094884037582282884e-6)  /* ln(2)^8 /8! */
151 };
152 #else
153 RAM_ALIGN
154 LNK_SECTION_CONSTDATA_L1
155 static const FIXP_SGL pow2Coeff[MAX_POW2_PRECISION] = {
156     FL2FXCONST_SGL(0.693147180559945309417232121458177),   /* ln(2)^1 /1! */
157     FL2FXCONST_SGL(0.240226506959100712333551263163332),   /* ln(2)^2 /2! */
158     FL2FXCONST_SGL(0.0555041086648215799531422637686218),  /* ln(2)^3 /3! */
159     FL2FXCONST_SGL(0.00961812910762847716197907157365887), /* ln(2)^4 /4! */
160     FL2FXCONST_SGL(0.00133335581464284434234122219879962), /* ln(2)^5 /5! */
161     FL2FXCONST_SGL(1.54035303933816099544370973327423e-4), /* ln(2)^6 /6! */
162     FL2FXCONST_SGL(1.52527338040598402800254390120096e-5), /* ln(2)^7 /7! */
163     FL2FXCONST_SGL(1.32154867901443094884037582282884e-6)  /* ln(2)^8 /8! */
164 };
165 #endif
166 
167 /*****************************************************************************
168 
169   functionname: CalcInvLdData
170   description:  Delivers the inverse of function CalcLdData().
171                 Delivers 2^(op*LD_DATA_SCALING)
172   input:        Input op is assumed to be fractional -1.0 < op < 1.0
173   output:       For op == 0, the result is MAXVAL_DBL (almost 1.0).
174                 For negative input values the output should be treated as a
175 positive fractional value. For positive input values the output should be
176 treated as a positive integer value. This function does not output negative
177 values.
178 
179 *****************************************************************************/
180 /* Date: 06-JULY-2012 Arthur Tritthart, IIS Fraunhofer Erlangen */
181 /* Version with 3 table lookup and 1 linear interpolations      */
182 /* Algorithm: compute power of 2, argument x is in Q7.25 format */
183 /*  result = 2^(x/64)                                           */
184 /*  We split exponent (x/64) into 5 components:                 */
185 /*  integer part:      represented by b31..b25  (exp)           */
186 /*  fractional part 1: represented by b24..b20  (lookup1)       */
187 /*  fractional part 2: represented by b19..b15  (lookup2)       */
188 /*  fractional part 3: represented by b14..b10  (lookup3)       */
189 /*  fractional part 4: represented by b09..b00  (frac)          */
190 /*  => result = (lookup1*lookup2*(lookup3+C1*frac)<<3)>>exp     */
191 /* Due to the fact, that all lookup values contain a factor 0.5 */
192 /* the result has to be shifted by 3 to the right also.         */
193 /* Table exp2_tab_long contains the log2 for 0 to 1.0 in steps  */
194 /* of 1/32, table exp2w_tab_long the log2 for 0 to 1/32 in steps*/
195 /* of 1/1024, table exp2x_tab_long the log2 for 0 to 1/1024 in  */
196 /* steps of 1/32768. Since the 2-logarithm of very very small   */
197 /* negative value is rather linear, we can use interpolation.   */
198 /* Limitations:                                                 */
199 /* For x <= 0, the result is fractional positive                */
200 /* For x > 0, the result is integer in range 1...7FFF.FFFF      */
201 /* For x < -31/64, we have to clear the result                  */
202 /* For x = 0, the result is ~1.0 (0x7FFF.FFFF)                  */
203 /* For x >= 31/64, the result is 0x7FFF.FFFF                    */
204 
205 /* This table is used for lookup 2^x with   */
206 /* x in range [0...1.0[ in steps of 1/32 */
207 LNK_SECTION_DATA_L1
208 const UINT exp2_tab_long[32] = {
209     0x40000000, 0x4166C34C, 0x42D561B4, 0x444C0740, 0x45CAE0F2, 0x47521CC6,
210     0x48E1E9BA, 0x4A7A77D4, 0x4C1BF829, 0x4DC69CDD, 0x4F7A9930, 0x51382182,
211     0x52FF6B55, 0x54D0AD5A, 0x56AC1F75, 0x5891FAC1, 0x5A82799A, 0x5C7DD7A4,
212     0x5E8451D0, 0x60962665, 0x62B39509, 0x64DCDEC3, 0x6712460B, 0x69540EC9,
213     0x6BA27E65, 0x6DFDDBCC, 0x70666F76, 0x72DC8374, 0x75606374, 0x77F25CCE,
214     0x7A92BE8B, 0x7D41D96E
215     // 0x80000000
216 };
217 
218 /* This table is used for lookup 2^x with   */
219 /* x in range [0...1/32[ in steps of 1/1024 */
220 LNK_SECTION_DATA_L1
221 const UINT exp2w_tab_long[32] = {
222     0x40000000, 0x400B1818, 0x4016321B, 0x40214E0C, 0x402C6BE9, 0x40378BB4,
223     0x4042AD6D, 0x404DD113, 0x4058F6A8, 0x40641E2B, 0x406F479E, 0x407A7300,
224     0x4085A051, 0x4090CF92, 0x409C00C4, 0x40A733E6, 0x40B268FA, 0x40BD9FFF,
225     0x40C8D8F5, 0x40D413DD, 0x40DF50B8, 0x40EA8F86, 0x40F5D046, 0x410112FA,
226     0x410C57A2, 0x41179E3D, 0x4122E6CD, 0x412E3152, 0x41397DCC, 0x4144CC3B,
227     0x41501CA0, 0x415B6EFB,
228     // 0x4166C34C,
229 };
230 /* This table is used for lookup 2^x with   */
231 /* x in range [0...1/1024[ in steps of 1/32768 */
232 LNK_SECTION_DATA_L1
233 const UINT exp2x_tab_long[32] = {
234     0x40000000, 0x400058B9, 0x4000B173, 0x40010A2D, 0x400162E8, 0x4001BBA3,
235     0x4002145F, 0x40026D1B, 0x4002C5D8, 0x40031E95, 0x40037752, 0x4003D011,
236     0x400428CF, 0x4004818E, 0x4004DA4E, 0x4005330E, 0x40058BCE, 0x4005E48F,
237     0x40063D51, 0x40069613, 0x4006EED5, 0x40074798, 0x4007A05B, 0x4007F91F,
238     0x400851E4, 0x4008AAA8, 0x4009036E, 0x40095C33, 0x4009B4FA, 0x400A0DC0,
239     0x400A6688, 0x400ABF4F,
240     // 0x400B1818
241 };
242 
243 /*****************************************************************************
244     functionname: InitLdInt and CalcLdInt
245     description:  Create and access table with integer LdData (0 to
246 LD_INT_TAB_LEN)
247 *****************************************************************************/
248 #ifndef LD_INT_TAB_LEN
249 #define LD_INT_TAB_LEN \
250   193 /* Default tab length. Lower value should be set in fix.h */
251 #endif
252 
253 #if (LD_INT_TAB_LEN <= 120)
254 LNK_SECTION_CONSTDATA_L1
255 static const FIXP_DBL ldIntCoeff[] = {
256     (FIXP_DBL)0x80000001, (FIXP_DBL)0x00000000, (FIXP_DBL)0x02000000,
257     (FIXP_DBL)0x032b8034, (FIXP_DBL)0x04000000, (FIXP_DBL)0x04a4d3c2,
258     (FIXP_DBL)0x052b8034, (FIXP_DBL)0x059d5da0, (FIXP_DBL)0x06000000,
259     (FIXP_DBL)0x06570069, (FIXP_DBL)0x06a4d3c2, (FIXP_DBL)0x06eb3a9f,
260     (FIXP_DBL)0x072b8034, (FIXP_DBL)0x0766a009, (FIXP_DBL)0x079d5da0,
261     (FIXP_DBL)0x07d053f7, (FIXP_DBL)0x08000000, (FIXP_DBL)0x082cc7ee,
262     (FIXP_DBL)0x08570069, (FIXP_DBL)0x087ef05b, (FIXP_DBL)0x08a4d3c2,
263     (FIXP_DBL)0x08c8ddd4, (FIXP_DBL)0x08eb3a9f, (FIXP_DBL)0x090c1050,
264     (FIXP_DBL)0x092b8034, (FIXP_DBL)0x0949a785, (FIXP_DBL)0x0966a009,
265     (FIXP_DBL)0x0982809d, (FIXP_DBL)0x099d5da0, (FIXP_DBL)0x09b74949,
266     (FIXP_DBL)0x09d053f7, (FIXP_DBL)0x09e88c6b, (FIXP_DBL)0x0a000000,
267     (FIXP_DBL)0x0a16bad3, (FIXP_DBL)0x0a2cc7ee, (FIXP_DBL)0x0a423162,
268     (FIXP_DBL)0x0a570069, (FIXP_DBL)0x0a6b3d79, (FIXP_DBL)0x0a7ef05b,
269     (FIXP_DBL)0x0a92203d, (FIXP_DBL)0x0aa4d3c2, (FIXP_DBL)0x0ab7110e,
270     (FIXP_DBL)0x0ac8ddd4, (FIXP_DBL)0x0ada3f60, (FIXP_DBL)0x0aeb3a9f,
271     (FIXP_DBL)0x0afbd42b, (FIXP_DBL)0x0b0c1050, (FIXP_DBL)0x0b1bf312,
272     (FIXP_DBL)0x0b2b8034, (FIXP_DBL)0x0b3abb40, (FIXP_DBL)0x0b49a785,
273     (FIXP_DBL)0x0b584822, (FIXP_DBL)0x0b66a009, (FIXP_DBL)0x0b74b1fd,
274     (FIXP_DBL)0x0b82809d, (FIXP_DBL)0x0b900e61, (FIXP_DBL)0x0b9d5da0,
275     (FIXP_DBL)0x0baa708f, (FIXP_DBL)0x0bb74949, (FIXP_DBL)0x0bc3e9ca,
276     (FIXP_DBL)0x0bd053f7, (FIXP_DBL)0x0bdc899b, (FIXP_DBL)0x0be88c6b,
277     (FIXP_DBL)0x0bf45e09, (FIXP_DBL)0x0c000000, (FIXP_DBL)0x0c0b73cb,
278     (FIXP_DBL)0x0c16bad3, (FIXP_DBL)0x0c21d671, (FIXP_DBL)0x0c2cc7ee,
279     (FIXP_DBL)0x0c379085, (FIXP_DBL)0x0c423162, (FIXP_DBL)0x0c4caba8,
280     (FIXP_DBL)0x0c570069, (FIXP_DBL)0x0c6130af, (FIXP_DBL)0x0c6b3d79,
281     (FIXP_DBL)0x0c7527b9, (FIXP_DBL)0x0c7ef05b, (FIXP_DBL)0x0c88983f,
282     (FIXP_DBL)0x0c92203d, (FIXP_DBL)0x0c9b8926, (FIXP_DBL)0x0ca4d3c2,
283     (FIXP_DBL)0x0cae00d2, (FIXP_DBL)0x0cb7110e, (FIXP_DBL)0x0cc0052b,
284     (FIXP_DBL)0x0cc8ddd4, (FIXP_DBL)0x0cd19bb0, (FIXP_DBL)0x0cda3f60,
285     (FIXP_DBL)0x0ce2c97d, (FIXP_DBL)0x0ceb3a9f, (FIXP_DBL)0x0cf39355,
286     (FIXP_DBL)0x0cfbd42b, (FIXP_DBL)0x0d03fda9, (FIXP_DBL)0x0d0c1050,
287     (FIXP_DBL)0x0d140ca0, (FIXP_DBL)0x0d1bf312, (FIXP_DBL)0x0d23c41d,
288     (FIXP_DBL)0x0d2b8034, (FIXP_DBL)0x0d3327c7, (FIXP_DBL)0x0d3abb40,
289     (FIXP_DBL)0x0d423b08, (FIXP_DBL)0x0d49a785, (FIXP_DBL)0x0d510118,
290     (FIXP_DBL)0x0d584822, (FIXP_DBL)0x0d5f7cff, (FIXP_DBL)0x0d66a009,
291     (FIXP_DBL)0x0d6db197, (FIXP_DBL)0x0d74b1fd, (FIXP_DBL)0x0d7ba190,
292     (FIXP_DBL)0x0d82809d, (FIXP_DBL)0x0d894f75, (FIXP_DBL)0x0d900e61,
293     (FIXP_DBL)0x0d96bdad, (FIXP_DBL)0x0d9d5da0, (FIXP_DBL)0x0da3ee7f,
294     (FIXP_DBL)0x0daa708f, (FIXP_DBL)0x0db0e412, (FIXP_DBL)0x0db74949,
295     (FIXP_DBL)0x0dbda072, (FIXP_DBL)0x0dc3e9ca, (FIXP_DBL)0x0dca258e};
296 
297 #elif (LD_INT_TAB_LEN <= 193)
298 LNK_SECTION_CONSTDATA_L1
299 static const FIXP_DBL ldIntCoeff[] = {
300     (FIXP_DBL)0x80000001, (FIXP_DBL)0x00000000, (FIXP_DBL)0x02000000,
301     (FIXP_DBL)0x032b8034, (FIXP_DBL)0x04000000, (FIXP_DBL)0x04a4d3c2,
302     (FIXP_DBL)0x052b8034, (FIXP_DBL)0x059d5da0, (FIXP_DBL)0x06000000,
303     (FIXP_DBL)0x06570069, (FIXP_DBL)0x06a4d3c2, (FIXP_DBL)0x06eb3a9f,
304     (FIXP_DBL)0x072b8034, (FIXP_DBL)0x0766a009, (FIXP_DBL)0x079d5da0,
305     (FIXP_DBL)0x07d053f7, (FIXP_DBL)0x08000000, (FIXP_DBL)0x082cc7ee,
306     (FIXP_DBL)0x08570069, (FIXP_DBL)0x087ef05b, (FIXP_DBL)0x08a4d3c2,
307     (FIXP_DBL)0x08c8ddd4, (FIXP_DBL)0x08eb3a9f, (FIXP_DBL)0x090c1050,
308     (FIXP_DBL)0x092b8034, (FIXP_DBL)0x0949a785, (FIXP_DBL)0x0966a009,
309     (FIXP_DBL)0x0982809d, (FIXP_DBL)0x099d5da0, (FIXP_DBL)0x09b74949,
310     (FIXP_DBL)0x09d053f7, (FIXP_DBL)0x09e88c6b, (FIXP_DBL)0x0a000000,
311     (FIXP_DBL)0x0a16bad3, (FIXP_DBL)0x0a2cc7ee, (FIXP_DBL)0x0a423162,
312     (FIXP_DBL)0x0a570069, (FIXP_DBL)0x0a6b3d79, (FIXP_DBL)0x0a7ef05b,
313     (FIXP_DBL)0x0a92203d, (FIXP_DBL)0x0aa4d3c2, (FIXP_DBL)0x0ab7110e,
314     (FIXP_DBL)0x0ac8ddd4, (FIXP_DBL)0x0ada3f60, (FIXP_DBL)0x0aeb3a9f,
315     (FIXP_DBL)0x0afbd42b, (FIXP_DBL)0x0b0c1050, (FIXP_DBL)0x0b1bf312,
316     (FIXP_DBL)0x0b2b8034, (FIXP_DBL)0x0b3abb40, (FIXP_DBL)0x0b49a785,
317     (FIXP_DBL)0x0b584822, (FIXP_DBL)0x0b66a009, (FIXP_DBL)0x0b74b1fd,
318     (FIXP_DBL)0x0b82809d, (FIXP_DBL)0x0b900e61, (FIXP_DBL)0x0b9d5da0,
319     (FIXP_DBL)0x0baa708f, (FIXP_DBL)0x0bb74949, (FIXP_DBL)0x0bc3e9ca,
320     (FIXP_DBL)0x0bd053f7, (FIXP_DBL)0x0bdc899b, (FIXP_DBL)0x0be88c6b,
321     (FIXP_DBL)0x0bf45e09, (FIXP_DBL)0x0c000000, (FIXP_DBL)0x0c0b73cb,
322     (FIXP_DBL)0x0c16bad3, (FIXP_DBL)0x0c21d671, (FIXP_DBL)0x0c2cc7ee,
323     (FIXP_DBL)0x0c379085, (FIXP_DBL)0x0c423162, (FIXP_DBL)0x0c4caba8,
324     (FIXP_DBL)0x0c570069, (FIXP_DBL)0x0c6130af, (FIXP_DBL)0x0c6b3d79,
325     (FIXP_DBL)0x0c7527b9, (FIXP_DBL)0x0c7ef05b, (FIXP_DBL)0x0c88983f,
326     (FIXP_DBL)0x0c92203d, (FIXP_DBL)0x0c9b8926, (FIXP_DBL)0x0ca4d3c2,
327     (FIXP_DBL)0x0cae00d2, (FIXP_DBL)0x0cb7110e, (FIXP_DBL)0x0cc0052b,
328     (FIXP_DBL)0x0cc8ddd4, (FIXP_DBL)0x0cd19bb0, (FIXP_DBL)0x0cda3f60,
329     (FIXP_DBL)0x0ce2c97d, (FIXP_DBL)0x0ceb3a9f, (FIXP_DBL)0x0cf39355,
330     (FIXP_DBL)0x0cfbd42b, (FIXP_DBL)0x0d03fda9, (FIXP_DBL)0x0d0c1050,
331     (FIXP_DBL)0x0d140ca0, (FIXP_DBL)0x0d1bf312, (FIXP_DBL)0x0d23c41d,
332     (FIXP_DBL)0x0d2b8034, (FIXP_DBL)0x0d3327c7, (FIXP_DBL)0x0d3abb40,
333     (FIXP_DBL)0x0d423b08, (FIXP_DBL)0x0d49a785, (FIXP_DBL)0x0d510118,
334     (FIXP_DBL)0x0d584822, (FIXP_DBL)0x0d5f7cff, (FIXP_DBL)0x0d66a009,
335     (FIXP_DBL)0x0d6db197, (FIXP_DBL)0x0d74b1fd, (FIXP_DBL)0x0d7ba190,
336     (FIXP_DBL)0x0d82809d, (FIXP_DBL)0x0d894f75, (FIXP_DBL)0x0d900e61,
337     (FIXP_DBL)0x0d96bdad, (FIXP_DBL)0x0d9d5da0, (FIXP_DBL)0x0da3ee7f,
338     (FIXP_DBL)0x0daa708f, (FIXP_DBL)0x0db0e412, (FIXP_DBL)0x0db74949,
339     (FIXP_DBL)0x0dbda072, (FIXP_DBL)0x0dc3e9ca, (FIXP_DBL)0x0dca258e,
340     (FIXP_DBL)0x0dd053f7, (FIXP_DBL)0x0dd6753e, (FIXP_DBL)0x0ddc899b,
341     (FIXP_DBL)0x0de29143, (FIXP_DBL)0x0de88c6b, (FIXP_DBL)0x0dee7b47,
342     (FIXP_DBL)0x0df45e09, (FIXP_DBL)0x0dfa34e1, (FIXP_DBL)0x0e000000,
343     (FIXP_DBL)0x0e05bf94, (FIXP_DBL)0x0e0b73cb, (FIXP_DBL)0x0e111cd2,
344     (FIXP_DBL)0x0e16bad3, (FIXP_DBL)0x0e1c4dfb, (FIXP_DBL)0x0e21d671,
345     (FIXP_DBL)0x0e275460, (FIXP_DBL)0x0e2cc7ee, (FIXP_DBL)0x0e323143,
346     (FIXP_DBL)0x0e379085, (FIXP_DBL)0x0e3ce5d8, (FIXP_DBL)0x0e423162,
347     (FIXP_DBL)0x0e477346, (FIXP_DBL)0x0e4caba8, (FIXP_DBL)0x0e51daa8,
348     (FIXP_DBL)0x0e570069, (FIXP_DBL)0x0e5c1d0b, (FIXP_DBL)0x0e6130af,
349     (FIXP_DBL)0x0e663b74, (FIXP_DBL)0x0e6b3d79, (FIXP_DBL)0x0e7036db,
350     (FIXP_DBL)0x0e7527b9, (FIXP_DBL)0x0e7a1030, (FIXP_DBL)0x0e7ef05b,
351     (FIXP_DBL)0x0e83c857, (FIXP_DBL)0x0e88983f, (FIXP_DBL)0x0e8d602e,
352     (FIXP_DBL)0x0e92203d, (FIXP_DBL)0x0e96d888, (FIXP_DBL)0x0e9b8926,
353     (FIXP_DBL)0x0ea03232, (FIXP_DBL)0x0ea4d3c2, (FIXP_DBL)0x0ea96df0,
354     (FIXP_DBL)0x0eae00d2, (FIXP_DBL)0x0eb28c7f, (FIXP_DBL)0x0eb7110e,
355     (FIXP_DBL)0x0ebb8e96, (FIXP_DBL)0x0ec0052b, (FIXP_DBL)0x0ec474e4,
356     (FIXP_DBL)0x0ec8ddd4, (FIXP_DBL)0x0ecd4012, (FIXP_DBL)0x0ed19bb0,
357     (FIXP_DBL)0x0ed5f0c4, (FIXP_DBL)0x0eda3f60, (FIXP_DBL)0x0ede8797,
358     (FIXP_DBL)0x0ee2c97d, (FIXP_DBL)0x0ee70525, (FIXP_DBL)0x0eeb3a9f,
359     (FIXP_DBL)0x0eef69ff, (FIXP_DBL)0x0ef39355, (FIXP_DBL)0x0ef7b6b4,
360     (FIXP_DBL)0x0efbd42b, (FIXP_DBL)0x0effebcd, (FIXP_DBL)0x0f03fda9,
361     (FIXP_DBL)0x0f0809cf, (FIXP_DBL)0x0f0c1050, (FIXP_DBL)0x0f10113b,
362     (FIXP_DBL)0x0f140ca0, (FIXP_DBL)0x0f18028d, (FIXP_DBL)0x0f1bf312,
363     (FIXP_DBL)0x0f1fde3d, (FIXP_DBL)0x0f23c41d, (FIXP_DBL)0x0f27a4c0,
364     (FIXP_DBL)0x0f2b8034};
365 
366 #else
367 #error "ldInt table size too small"
368 
369 #endif
370 
371 LNK_SECTION_INITCODE
InitLdInt()372 void InitLdInt() { /* nothing to do! Use preinitialized logarithm table */
373 }
374 
375 #if (LD_INT_TAB_LEN != 0)
376 
377 LNK_SECTION_CODE_L1
CalcLdInt(INT i)378 FIXP_DBL CalcLdInt(INT i) {
379   /* calculates ld(op)/LD_DATA_SCALING */
380   /* op is assumed to be an integer value between 1 and LD_INT_TAB_LEN */
381 
382   FDK_ASSERT((LD_INT_TAB_LEN > 0) &&
383              ((FIXP_DBL)ldIntCoeff[0] ==
384               (FIXP_DBL)0x80000001)); /* tab has to be initialized */
385 
386   if ((i > 0) && (i < LD_INT_TAB_LEN))
387     return ldIntCoeff[i];
388   else {
389     return (0);
390   }
391 }
392 #endif /* (LD_INT_TAB_LEN!=0)  */
393 
394 #if !defined(FUNCTION_schur_div)
395 /*****************************************************************************
396 
397     functionname: schur_div
398     description:  delivers op1/op2 with op3-bit accuracy
399 
400 *****************************************************************************/
401 
schur_div(FIXP_DBL num,FIXP_DBL denum,INT count)402 FIXP_DBL schur_div(FIXP_DBL num, FIXP_DBL denum, INT count) {
403   INT L_num = (LONG)num >> 1;
404   INT L_denum = (LONG)denum >> 1;
405   INT div = 0;
406   INT k = count;
407 
408   FDK_ASSERT(num >= (FIXP_DBL)0);
409   FDK_ASSERT(denum > (FIXP_DBL)0);
410   FDK_ASSERT(num <= denum);
411 
412   if (L_num != 0)
413     while (--k) {
414       div <<= 1;
415       L_num <<= 1;
416       if (L_num >= L_denum) {
417         L_num -= L_denum;
418         div++;
419       }
420     }
421   return (FIXP_DBL)(div << (DFRACT_BITS - count));
422 }
423 
424 #endif /* !defined(FUNCTION_schur_div) */
425 
426 #ifndef FUNCTION_fMultNorm
fMultNorm(FIXP_DBL f1,FIXP_DBL f2,INT * result_e)427 FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2, INT *result_e) {
428   INT product = 0;
429   INT norm_f1, norm_f2;
430 
431   if ((f1 == (FIXP_DBL)0) || (f2 == (FIXP_DBL)0)) {
432     *result_e = 0;
433     return (FIXP_DBL)0;
434   }
435   norm_f1 = CountLeadingBits(f1);
436   f1 = f1 << norm_f1;
437   norm_f2 = CountLeadingBits(f2);
438   f2 = f2 << norm_f2;
439 
440   if ((f1 == (FIXP_DBL)MINVAL_DBL) && (f2 == (FIXP_DBL)MINVAL_DBL)) {
441     product = -((FIXP_DBL)MINVAL_DBL >> 1);
442     *result_e = -(norm_f1 + norm_f2 - 1);
443   } else {
444     product = fMult(f1, f2);
445     *result_e = -(norm_f1 + norm_f2);
446   }
447 
448   return (FIXP_DBL)product;
449 }
450 #endif
451 
452 #ifndef FUNCTION_fDivNorm
fDivNorm(FIXP_DBL L_num,FIXP_DBL L_denum,INT * result_e)453 FIXP_DBL fDivNorm(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e) {
454   FIXP_DBL div;
455   INT norm_num, norm_den;
456 
457   FDK_ASSERT(L_num >= (FIXP_DBL)0);
458   FDK_ASSERT(L_denum > (FIXP_DBL)0);
459 
460   if (L_num == (FIXP_DBL)0) {
461     *result_e = 0;
462     return ((FIXP_DBL)0);
463   }
464 
465   norm_num = CountLeadingBits(L_num);
466   L_num = L_num << norm_num;
467   L_num = L_num >> 1;
468   *result_e = -norm_num + 1;
469 
470   norm_den = CountLeadingBits(L_denum);
471   L_denum = L_denum << norm_den;
472   *result_e -= -norm_den;
473 
474   div = schur_div(L_num, L_denum, FRACT_BITS);
475 
476   return div;
477 }
478 #endif /* !FUNCTION_fDivNorm */
479 
480 #ifndef FUNCTION_fDivNorm
fDivNorm(FIXP_DBL num,FIXP_DBL denom)481 FIXP_DBL fDivNorm(FIXP_DBL num, FIXP_DBL denom) {
482   INT e;
483   FIXP_DBL res;
484 
485   FDK_ASSERT(denom >= num);
486 
487   res = fDivNorm(num, denom, &e);
488 
489   /* Avoid overflow since we must output a value with exponent 0
490      there is no other choice than saturating to almost 1.0f */
491   if (res == (FIXP_DBL)(1 << (DFRACT_BITS - 2)) && e == 1) {
492     res = (FIXP_DBL)MAXVAL_DBL;
493   } else {
494     res = scaleValue(res, e);
495   }
496 
497   return res;
498 }
499 #endif /* !FUNCTION_fDivNorm */
500 
501 #ifndef FUNCTION_fDivNormSigned
fDivNormSigned(FIXP_DBL num,FIXP_DBL denom)502 FIXP_DBL fDivNormSigned(FIXP_DBL num, FIXP_DBL denom) {
503   INT e;
504   FIXP_DBL res;
505   int sign;
506 
507   if (denom == (FIXP_DBL)0) {
508     return (FIXP_DBL)MAXVAL_DBL;
509   }
510 
511   sign = ((num >= (FIXP_DBL)0) != (denom >= (FIXP_DBL)0));
512   res = fDivNormSigned(num, denom, &e);
513 
514   /* Saturate since we must output a value with exponent 0 */
515   if ((e > 0) && (fAbs(res) >= FL2FXCONST_DBL(0.5))) {
516     if (sign) {
517       res = (FIXP_DBL)MINVAL_DBL;
518     } else {
519       res = (FIXP_DBL)MAXVAL_DBL;
520     }
521   } else {
522     res = scaleValue(res, e);
523   }
524 
525   return res;
526 }
fDivNormSigned(FIXP_DBL L_num,FIXP_DBL L_denum,INT * result_e)527 FIXP_DBL fDivNormSigned(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e) {
528   FIXP_DBL div;
529   INT norm_num, norm_den;
530   int sign;
531 
532   sign = ((L_num >= (FIXP_DBL)0) != (L_denum >= (FIXP_DBL)0));
533 
534   if (L_num == (FIXP_DBL)0) {
535     *result_e = 0;
536     return ((FIXP_DBL)0);
537   }
538   if (L_denum == (FIXP_DBL)0) {
539     *result_e = 14;
540     return ((FIXP_DBL)MAXVAL_DBL);
541   }
542 
543   norm_num = CountLeadingBits(L_num);
544   L_num = L_num << norm_num;
545   L_num = L_num >> 2;
546   L_num = fAbs(L_num);
547   *result_e = -norm_num + 1;
548 
549   norm_den = CountLeadingBits(L_denum);
550   L_denum = L_denum << norm_den;
551   L_denum = L_denum >> 1;
552   L_denum = fAbs(L_denum);
553   *result_e -= -norm_den;
554 
555   div = schur_div(L_num, L_denum, FRACT_BITS);
556 
557   if (sign) {
558     div = -div;
559   }
560 
561   return div;
562 }
563 #endif /* FUNCTION_fDivNormSigned */
564 
565 #ifndef FUNCTION_fDivNormHighPrec
fDivNormHighPrec(FIXP_DBL num,FIXP_DBL denom,INT * result_e)566 FIXP_DBL fDivNormHighPrec(FIXP_DBL num, FIXP_DBL denom, INT *result_e) {
567   FIXP_DBL div;
568   INT norm_num, norm_den;
569 
570   FDK_ASSERT(num >= (FIXP_DBL)0);
571   FDK_ASSERT(denom > (FIXP_DBL)0);
572 
573   if (num == (FIXP_DBL)0) {
574     *result_e = 0;
575     return ((FIXP_DBL)0);
576   }
577 
578   norm_num = CountLeadingBits(num);
579   num = num << norm_num;
580   num = num >> 1;
581   *result_e = -norm_num + 1;
582 
583   norm_den = CountLeadingBits(denom);
584   denom = denom << norm_den;
585   *result_e -= -norm_den;
586 
587   div = schur_div(num, denom, 31);
588   return div;
589 }
590 #endif /* !FUNCTION_fDivNormHighPrec */
591 
592 #ifndef FUNCTION_fPow
f2Pow(const FIXP_DBL exp_m,const INT exp_e,INT * result_e)593 FIXP_DBL f2Pow(const FIXP_DBL exp_m, const INT exp_e, INT *result_e) {
594   FIXP_DBL frac_part, result_m;
595   INT int_part;
596 
597   if (exp_e > 0) {
598     INT exp_bits = DFRACT_BITS - 1 - exp_e;
599     int_part = exp_m >> exp_bits;
600     frac_part = exp_m - (FIXP_DBL)(int_part << exp_bits);
601     frac_part = frac_part << exp_e;
602   } else {
603     int_part = 0;
604     frac_part = exp_m >> -exp_e;
605   }
606 
607   /* Best accuracy is around 0, so try to get there with the fractional part. */
608   if (frac_part > FL2FXCONST_DBL(0.5f)) {
609     int_part = int_part + 1;
610     frac_part = frac_part + FL2FXCONST_DBL(-1.0f);
611   }
612   if (frac_part < FL2FXCONST_DBL(-0.5f)) {
613     int_part = int_part - 1;
614     frac_part = -(FL2FXCONST_DBL(-1.0f) - frac_part);
615   }
616 
617   /* "+ 1" compensates fMultAddDiv2() of the polynomial evaluation below. */
618   *result_e = int_part + 1;
619 
620   /* Evaluate taylor polynomial which approximates 2^x */
621   {
622     FIXP_DBL p;
623 
624     /* result_m ~= 2^frac_part */
625     p = frac_part;
626     /* First taylor series coefficient a_0 = 1.0, scaled by 0.5 due to
627      * fMultDiv2(). */
628     result_m = FL2FXCONST_DBL(1.0f / 2.0f);
629     for (INT i = 0; i < POW2_PRECISION; i++) {
630       /* next taylor series term: a_i * x^i, x=0 */
631       result_m = fMultAddDiv2(result_m, pow2Coeff[i], p);
632       p = fMult(p, frac_part);
633     }
634   }
635   return result_m;
636 }
637 
f2Pow(const FIXP_DBL exp_m,const INT exp_e)638 FIXP_DBL f2Pow(const FIXP_DBL exp_m, const INT exp_e) {
639   FIXP_DBL result_m;
640   INT result_e;
641 
642   result_m = f2Pow(exp_m, exp_e, &result_e);
643   result_e = fixMin(DFRACT_BITS - 1, fixMax(-(DFRACT_BITS - 1), result_e));
644 
645   return scaleValue(result_m, result_e);
646 }
647 
fPow(FIXP_DBL base_m,INT base_e,FIXP_DBL exp_m,INT exp_e,INT * result_e)648 FIXP_DBL fPow(FIXP_DBL base_m, INT base_e, FIXP_DBL exp_m, INT exp_e,
649               INT *result_e) {
650   INT ans_lg2_e, baselg2_e;
651   FIXP_DBL base_lg2, ans_lg2, result;
652 
653   if (base_m <= (FIXP_DBL)0) {
654     result = (FIXP_DBL)0;
655     *result_e = 0;
656     return result;
657   }
658 
659   /* Calc log2 of base */
660   base_lg2 = fLog2(base_m, base_e, &baselg2_e);
661 
662   /* Prepare exp */
663   {
664     INT leadingBits;
665 
666     leadingBits = CountLeadingBits(fAbs(exp_m));
667     exp_m = exp_m << leadingBits;
668     exp_e -= leadingBits;
669   }
670 
671   /* Calc base pow exp */
672   ans_lg2 = fMult(base_lg2, exp_m);
673   ans_lg2_e = exp_e + baselg2_e;
674 
675   /* Calc antilog */
676   result = f2Pow(ans_lg2, ans_lg2_e, result_e);
677 
678   return result;
679 }
680 
fLdPow(FIXP_DBL baseLd_m,INT baseLd_e,FIXP_DBL exp_m,INT exp_e,INT * result_e)681 FIXP_DBL fLdPow(FIXP_DBL baseLd_m, INT baseLd_e, FIXP_DBL exp_m, INT exp_e,
682                 INT *result_e) {
683   INT ans_lg2_e;
684   FIXP_DBL ans_lg2, result;
685 
686   /* Prepare exp */
687   {
688     INT leadingBits;
689 
690     leadingBits = CountLeadingBits(fAbs(exp_m));
691     exp_m = exp_m << leadingBits;
692     exp_e -= leadingBits;
693   }
694 
695   /* Calc base pow exp */
696   ans_lg2 = fMult(baseLd_m, exp_m);
697   ans_lg2_e = exp_e + baseLd_e;
698 
699   /* Calc antilog */
700   result = f2Pow(ans_lg2, ans_lg2_e, result_e);
701 
702   return result;
703 }
704 
fLdPow(FIXP_DBL baseLd_m,INT baseLd_e,FIXP_DBL exp_m,INT exp_e)705 FIXP_DBL fLdPow(FIXP_DBL baseLd_m, INT baseLd_e, FIXP_DBL exp_m, INT exp_e) {
706   FIXP_DBL result_m;
707   int result_e;
708 
709   result_m = fLdPow(baseLd_m, baseLd_e, exp_m, exp_e, &result_e);
710 
711   return SATURATE_SHIFT(result_m, -result_e, DFRACT_BITS);
712 }
713 
fPowInt(FIXP_DBL base_m,INT base_e,INT exp,INT * pResult_e)714 FIXP_DBL fPowInt(FIXP_DBL base_m, INT base_e, INT exp, INT *pResult_e) {
715   FIXP_DBL result;
716 
717   if (exp != 0) {
718     INT result_e = 0;
719 
720     if (base_m != (FIXP_DBL)0) {
721       {
722         INT leadingBits;
723         leadingBits = CountLeadingBits(base_m);
724         base_m <<= leadingBits;
725         base_e -= leadingBits;
726       }
727 
728       result = base_m;
729 
730       {
731         int i;
732         for (i = 1; i < fAbs(exp); i++) {
733           result = fMult(result, base_m);
734         }
735       }
736 
737       if (exp < 0) {
738         /* 1.0 / ans */
739         result = fDivNorm(FL2FXCONST_DBL(0.5f), result, &result_e);
740         result_e++;
741       } else {
742         int ansScale = CountLeadingBits(result);
743         result <<= ansScale;
744         result_e -= ansScale;
745       }
746 
747       result_e += exp * base_e;
748 
749     } else {
750       result = (FIXP_DBL)0;
751     }
752     *pResult_e = result_e;
753   } else {
754     result = FL2FXCONST_DBL(0.5f);
755     *pResult_e = 1;
756   }
757 
758   return result;
759 }
760 #endif /* FUNCTION_fPow */
761 
762 #ifndef FUNCTION_fLog2
CalcLog2(FIXP_DBL base_m,INT base_e,INT * result_e)763 FIXP_DBL CalcLog2(FIXP_DBL base_m, INT base_e, INT *result_e) {
764   return fLog2(base_m, base_e, result_e);
765 }
766 #endif /* FUNCTION_fLog2 */
767 
fixp_floorToInt(FIXP_DBL f_inp,INT sf)768 INT fixp_floorToInt(FIXP_DBL f_inp, INT sf) {
769   FDK_ASSERT(sf >= 0);
770   INT floorInt = (INT)(f_inp >> ((DFRACT_BITS - 1) - sf));
771   return floorInt;
772 }
773 
fixp_floor(FIXP_DBL f_inp,INT sf)774 FIXP_DBL fixp_floor(FIXP_DBL f_inp, INT sf) {
775   FDK_ASSERT(sf >= 0);
776   INT floorInt = fixp_floorToInt(f_inp, sf);
777   FIXP_DBL f_floor = (FIXP_DBL)(floorInt << ((DFRACT_BITS - 1) - sf));
778   return f_floor;
779 }
780 
fixp_ceilToInt(FIXP_DBL f_inp,INT sf)781 INT fixp_ceilToInt(FIXP_DBL f_inp, INT sf)  // sf  mantissaBits left of dot
782 {
783   FDK_ASSERT(sf >= 0);
784   INT sx = (DFRACT_BITS - 1) - sf;  // sx  mantissaBits right of dot
785   INT inpINT = (INT)f_inp;
786 
787   INT mask = (0x1 << sx) - 1;
788   INT ceilInt = (INT)(f_inp >> sx);
789 
790   if (inpINT & mask) {
791     ceilInt++;  // increment only, if there is at least one set mantissaBit
792                 // right of dot [in inpINT]
793   }
794 
795   return ceilInt;
796 }
797 
fixp_ceil(FIXP_DBL f_inp,INT sf)798 FIXP_DBL fixp_ceil(FIXP_DBL f_inp, INT sf) {
799   FDK_ASSERT(sf >= 0);
800   INT sx = (DFRACT_BITS - 1) - sf;
801   INT ceilInt = fixp_ceilToInt(f_inp, sf);
802   ULONG mask = (ULONG)0x1 << (DFRACT_BITS - 1);  // 0x80000000
803   ceilInt = ceilInt
804             << sx;  // no fract warn bec. shift into saturation done on int side
805 
806   if ((f_inp > FL2FXCONST_DBL(0.0f)) && (ceilInt & mask)) {
807     --ceilInt;
808   }
809   FIXP_DBL f_ceil = (FIXP_DBL)ceilInt;
810 
811   return f_ceil;
812 }
813 
814 /*****************************************************************************
815    fixp_truncateToInt()
816      Just remove the fractional part which is located right of decimal point
817      Same as which is done when a float is casted to (INT) :
818      result_INTtype = (INT)b_floatTypeInput;
819 
820    returns INT
821 *****************************************************************************/
fixp_truncateToInt(FIXP_DBL f_inp,INT sf)822 INT fixp_truncateToInt(FIXP_DBL f_inp, INT sf)  // sf  mantissaBits left  of dot
823                                                 // (without sign)  e.g. at width
824                                                 // 32 this would be [sign]7.
825                                                 // supposed sf equals 8.
826 {
827   FDK_ASSERT(sf >= 0);
828   INT sx = (DFRACT_BITS - 1) - sf;  // sx  mantissaBits right of dot
829                                     // at width 32 this would be        .24
830                                     // supposed sf equals 8.
831   INT fbaccu = (INT)f_inp;
832   INT mask = (0x1 << sx);
833 
834   if ((fbaccu < 0) && (fbaccu & (mask - 1))) {
835     fbaccu = fbaccu + mask;
836   }
837 
838   fbaccu = fbaccu >> sx;
839   return fbaccu;
840 }
841 
842 /*****************************************************************************
843    fixp_truncate()
844      Just remove the fractional part which is located right of decimal point
845 
846    returns FIXP_DBL
847 *****************************************************************************/
fixp_truncate(FIXP_DBL f_inp,INT sf)848 FIXP_DBL fixp_truncate(FIXP_DBL f_inp, INT sf) {
849   FDK_ASSERT(sf >= 0);
850   INT truncateInt = fixp_truncateToInt(f_inp, sf);
851   FIXP_DBL f_truncate = (FIXP_DBL)(truncateInt << ((DFRACT_BITS - 1) - sf));
852   return f_truncate;
853 }
854 
855 /*****************************************************************************
856   fixp_roundToInt()
857     round [typical rounding]
858 
859     See fct roundRef() [which is the reference]
860   returns INT
861 *****************************************************************************/
fixp_roundToInt(FIXP_DBL f_inp,INT sf)862 INT fixp_roundToInt(FIXP_DBL f_inp, INT sf) {
863   FDK_ASSERT(sf >= 0);
864   INT sx = DFRACT_BITS - 1 - sf;
865   INT inp = (INT)f_inp;
866   INT mask1 = (0x1 << (sx - 1));
867   INT mask2 = (0x1 << (sx)) - 1;
868   INT mask3 = 0x7FFFFFFF;
869   INT iam = inp & mask2;
870   INT rnd;
871 
872   if ((inp < 0) && !(iam == mask1))
873     rnd = inp + mask1;
874   else if ((inp > 0) && !(inp == mask3))
875     rnd = inp + mask1;
876   else
877     rnd = inp;
878 
879   rnd = rnd >> sx;
880 
881   if (inp == mask3) rnd++;
882 
883   return rnd;
884 }
885 
886 /*****************************************************************************
887   fixp_round()
888     round [typical rounding]
889 
890     See fct roundRef() [which is the reference]
891   returns FIXP_DBL
892 *****************************************************************************/
fixp_round(FIXP_DBL f_inp,INT sf)893 FIXP_DBL fixp_round(FIXP_DBL f_inp, INT sf) {
894   FDK_ASSERT(sf >= 0);
895   INT sx = DFRACT_BITS - 1 - sf;
896   INT r = fixp_roundToInt(f_inp, sf);
897   ULONG mask = (ULONG)0x1 << (DFRACT_BITS - 1);  // 0x80000000
898   r = r << sx;
899 
900   if ((f_inp > FL2FXCONST_DBL(0.0f)) && (r & mask)) {
901     --r;
902   }
903 
904   FIXP_DBL f_round = (FIXP_DBL)r;
905   return f_round;
906 }
907