1/* -----------------------------------------------------------------------------
2Software License for The Fraunhofer FDK AAC Codec Library for Android
3
4© Copyright  1995 - 2018 Fraunhofer-Gesellschaft zur Förderung der angewandten
5Forschung e.V. All rights reserved.
6
7 1.    INTRODUCTION
8The Fraunhofer FDK AAC Codec Library for Android ("FDK AAC Codec") is software
9that implements the MPEG Advanced Audio Coding ("AAC") encoding and decoding
10scheme for digital audio. This FDK AAC Codec software is intended to be used on
11a wide variety of Android devices.
12
13AAC's HE-AAC and HE-AAC v2 versions are regarded as today's most efficient
14general perceptual audio codecs. AAC-ELD is considered the best-performing
15full-bandwidth communications codec by independent studies and is widely
16deployed. AAC has been standardized by ISO and IEC as part of the MPEG
17specifications.
18
19Patent licenses for necessary patent claims for the FDK AAC Codec (including
20those of Fraunhofer) may be obtained through Via Licensing
21(www.vialicensing.com) or through the respective patent owners individually for
22the purpose of encoding or decoding bit streams in products that are compliant
23with the ISO/IEC MPEG audio standards. Please note that most manufacturers of
24Android devices already license these patent claims through Via Licensing or
25directly from the patent owners, and therefore FDK AAC Codec software may
26already be covered under those patent licenses when it is used for those
27licensed purposes only.
28
29Commercially-licensed AAC software libraries, including floating-point versions
30with enhanced sound quality, are also available from Fraunhofer. Users are
31encouraged to check the Fraunhofer website for additional applications
32information and documentation.
33
342.    COPYRIGHT LICENSE
35
36Redistribution and use in source and binary forms, with or without modification,
37are permitted without payment of copyright license fees provided that you
38satisfy the following conditions:
39
40You must retain the complete text of this software license in redistributions of
41the FDK AAC Codec or your modifications thereto in source code form.
42
43You must retain the complete text of this software license in the documentation
44and/or other materials provided with redistributions of the FDK AAC Codec or
45your modifications thereto in binary form. You must make available free of
46charge copies of the complete source code of the FDK AAC Codec and your
47modifications thereto to recipients of copies in binary form.
48
49The name of Fraunhofer may not be used to endorse or promote products derived
50from this library without prior written permission.
51
52You may not charge copyright license fees for anyone to use, copy or distribute
53the FDK AAC Codec software or your modifications thereto.
54
55Your modified versions of the FDK AAC Codec must carry prominent notices stating
56that you changed the software and the date of any change. For modified versions
57of the FDK AAC Codec, the term "Fraunhofer FDK AAC Codec Library for Android"
58must be replaced by the term "Third-Party Modified Version of the Fraunhofer FDK
59AAC Codec Library for Android."
60
613.    NO PATENT LICENSE
62
63NO EXPRESS OR IMPLIED LICENSES TO ANY PATENT CLAIMS, including without
64limitation the patents of Fraunhofer, ARE GRANTED BY THIS SOFTWARE LICENSE.
65Fraunhofer provides no warranty of patent non-infringement with respect to this
66software.
67
68You may use this FDK AAC Codec software or modifications thereto only for
69purposes that are authorized by appropriate patent licenses.
70
714.    DISCLAIMER
72
73This FDK AAC Codec software is provided by Fraunhofer on behalf of the copyright
74holders and contributors "AS IS" and WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES,
75including but not limited to the implied warranties of merchantability and
76fitness for a particular purpose. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
77CONTRIBUTORS BE LIABLE for any direct, indirect, incidental, special, exemplary,
78or consequential damages, including but not limited to procurement of substitute
79goods or services; loss of use, data, or profits, or business interruption,
80however caused and on any theory of liability, whether in contract, strict
81liability, or tort (including negligence), arising in any way out of the use of
82this software, even if advised of the possibility of such damage.
83
845.    CONTACT INFORMATION
85
86Fraunhofer Institute for Integrated Circuits IIS
87Attention: Audio and Multimedia Departments - FDK AAC LL
88Am Wolfsmantel 33
8991058 Erlangen, Germany
90
91www.iis.fraunhofer.de/amm
92amm-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*****************************************************************************/
116LNK_SECTION_CODE_L1
117void 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
140RAM_ALIGN
141LNK_SECTION_CONSTDATA_L1
142static 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
153RAM_ALIGN
154LNK_SECTION_CONSTDATA_L1
155static 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
175positive fractional value. For positive input values the output should be
176treated as a positive integer value. This function does not output negative
177values.
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 */
207LNK_SECTION_DATA_L1
208const 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 */
220LNK_SECTION_DATA_L1
221const 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 */
232LNK_SECTION_DATA_L1
233const 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
246LD_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)
254LNK_SECTION_CONSTDATA_L1
255static 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)
298LNK_SECTION_CONSTDATA_L1
299static 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
371LNK_SECTION_INITCODE
372void InitLdInt() { /* nothing to do! Use preinitialized logarithm table */
373}
374
375#if (LD_INT_TAB_LEN != 0)
376
377LNK_SECTION_CODE_L1
378FIXP_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
402FIXP_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
427FIXP_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
453FIXP_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
481FIXP_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
502FIXP_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}
527FIXP_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
566FIXP_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
593FIXP_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
638FIXP_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
648FIXP_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  /* Calc log2 of base */
654  base_lg2 = fLog2(base_m, base_e, &baselg2_e);
655
656  /* Prepare exp */
657  {
658    INT leadingBits;
659
660    leadingBits = CountLeadingBits(fAbs(exp_m));
661    exp_m = exp_m << leadingBits;
662    exp_e -= leadingBits;
663  }
664
665  /* Calc base pow exp */
666  ans_lg2 = fMult(base_lg2, exp_m);
667  ans_lg2_e = exp_e + baselg2_e;
668
669  /* Calc antilog */
670  result = f2Pow(ans_lg2, ans_lg2_e, result_e);
671
672  return result;
673}
674
675FIXP_DBL fLdPow(FIXP_DBL baseLd_m, INT baseLd_e, FIXP_DBL exp_m, INT exp_e,
676                INT *result_e) {
677  INT ans_lg2_e;
678  FIXP_DBL ans_lg2, result;
679
680  /* Prepare exp */
681  {
682    INT leadingBits;
683
684    leadingBits = CountLeadingBits(fAbs(exp_m));
685    exp_m = exp_m << leadingBits;
686    exp_e -= leadingBits;
687  }
688
689  /* Calc base pow exp */
690  ans_lg2 = fMult(baseLd_m, exp_m);
691  ans_lg2_e = exp_e + baseLd_e;
692
693  /* Calc antilog */
694  result = f2Pow(ans_lg2, ans_lg2_e, result_e);
695
696  return result;
697}
698
699FIXP_DBL fLdPow(FIXP_DBL baseLd_m, INT baseLd_e, FIXP_DBL exp_m, INT exp_e) {
700  FIXP_DBL result_m;
701  int result_e;
702
703  result_m = fLdPow(baseLd_m, baseLd_e, exp_m, exp_e, &result_e);
704
705  return SATURATE_SHIFT(result_m, -result_e, DFRACT_BITS);
706}
707
708FIXP_DBL fPowInt(FIXP_DBL base_m, INT base_e, INT exp, INT *pResult_e) {
709  FIXP_DBL result;
710
711  if (exp != 0) {
712    INT result_e = 0;
713
714    if (base_m != (FIXP_DBL)0) {
715      {
716        INT leadingBits;
717        leadingBits = CountLeadingBits(base_m);
718        base_m <<= leadingBits;
719        base_e -= leadingBits;
720      }
721
722      result = base_m;
723
724      {
725        int i;
726        for (i = 1; i < fAbs(exp); i++) {
727          result = fMult(result, base_m);
728        }
729      }
730
731      if (exp < 0) {
732        /* 1.0 / ans */
733        result = fDivNorm(FL2FXCONST_DBL(0.5f), result, &result_e);
734        result_e++;
735      } else {
736        int ansScale = CountLeadingBits(result);
737        result <<= ansScale;
738        result_e -= ansScale;
739      }
740
741      result_e += exp * base_e;
742
743    } else {
744      result = (FIXP_DBL)0;
745    }
746    *pResult_e = result_e;
747  } else {
748    result = FL2FXCONST_DBL(0.5f);
749    *pResult_e = 1;
750  }
751
752  return result;
753}
754#endif /* FUNCTION_fPow */
755
756#ifndef FUNCTION_fLog2
757FIXP_DBL CalcLog2(FIXP_DBL base_m, INT base_e, INT *result_e) {
758  return fLog2(base_m, base_e, result_e);
759}
760#endif /* FUNCTION_fLog2 */
761
762INT fixp_floorToInt(FIXP_DBL f_inp, INT sf) {
763  FDK_ASSERT(sf >= 0);
764  INT floorInt = (INT)(f_inp >> ((DFRACT_BITS - 1) - sf));
765  return floorInt;
766}
767
768FIXP_DBL fixp_floor(FIXP_DBL f_inp, INT sf) {
769  FDK_ASSERT(sf >= 0);
770  INT floorInt = fixp_floorToInt(f_inp, sf);
771  FIXP_DBL f_floor = (FIXP_DBL)(floorInt << ((DFRACT_BITS - 1) - sf));
772  return f_floor;
773}
774
775INT fixp_ceilToInt(FIXP_DBL f_inp, INT sf)  // sf  mantissaBits left of dot
776{
777  FDK_ASSERT(sf >= 0);
778  INT sx = (DFRACT_BITS - 1) - sf;  // sx  mantissaBits right of dot
779  INT inpINT = (INT)f_inp;
780
781  INT mask = (0x1 << sx) - 1;
782  INT ceilInt = (INT)(f_inp >> sx);
783
784  if (inpINT & mask) {
785    ceilInt++;  // increment only, if there is at least one set mantissaBit
786                // right of dot [in inpINT]
787  }
788
789  return ceilInt;
790}
791
792FIXP_DBL fixp_ceil(FIXP_DBL f_inp, INT sf) {
793  FDK_ASSERT(sf >= 0);
794  INT sx = (DFRACT_BITS - 1) - sf;
795  INT ceilInt = fixp_ceilToInt(f_inp, sf);
796  ULONG mask = (ULONG)0x1 << (DFRACT_BITS - 1);  // 0x80000000
797  ceilInt = ceilInt
798            << sx;  // no fract warn bec. shift into saturation done on int side
799
800  if ((f_inp > FL2FXCONST_DBL(0.0f)) && (ceilInt & mask)) {
801    --ceilInt;
802  }
803  FIXP_DBL f_ceil = (FIXP_DBL)ceilInt;
804
805  return f_ceil;
806}
807
808/*****************************************************************************
809   fixp_truncateToInt()
810     Just remove the fractional part which is located right of decimal point
811     Same as which is done when a float is casted to (INT) :
812     result_INTtype = (INT)b_floatTypeInput;
813
814   returns INT
815*****************************************************************************/
816INT fixp_truncateToInt(FIXP_DBL f_inp, INT sf)  // sf  mantissaBits left  of dot
817                                                // (without sign)  e.g. at width
818                                                // 32 this would be [sign]7.
819                                                // supposed sf equals 8.
820{
821  FDK_ASSERT(sf >= 0);
822  INT sx = (DFRACT_BITS - 1) - sf;  // sx  mantissaBits right of dot
823                                    // at width 32 this would be        .24
824                                    // supposed sf equals 8.
825  INT fbaccu = (INT)f_inp;
826  INT mask = (0x1 << sx);
827
828  if ((fbaccu < 0) && (fbaccu & (mask - 1))) {
829    fbaccu = fbaccu + mask;
830  }
831
832  fbaccu = fbaccu >> sx;
833  return fbaccu;
834}
835
836/*****************************************************************************
837   fixp_truncate()
838     Just remove the fractional part which is located right of decimal point
839
840   returns FIXP_DBL
841*****************************************************************************/
842FIXP_DBL fixp_truncate(FIXP_DBL f_inp, INT sf) {
843  FDK_ASSERT(sf >= 0);
844  INT truncateInt = fixp_truncateToInt(f_inp, sf);
845  FIXP_DBL f_truncate = (FIXP_DBL)(truncateInt << ((DFRACT_BITS - 1) - sf));
846  return f_truncate;
847}
848
849/*****************************************************************************
850  fixp_roundToInt()
851    round [typical rounding]
852
853    See fct roundRef() [which is the reference]
854  returns INT
855*****************************************************************************/
856INT fixp_roundToInt(FIXP_DBL f_inp, INT sf) {
857  FDK_ASSERT(sf >= 0);
858  INT sx = DFRACT_BITS - 1 - sf;
859  INT inp = (INT)f_inp;
860  INT mask1 = (0x1 << (sx - 1));
861  INT mask2 = (0x1 << (sx)) - 1;
862  INT mask3 = 0x7FFFFFFF;
863  INT iam = inp & mask2;
864  INT rnd;
865
866  if ((inp < 0) && !(iam == mask1))
867    rnd = inp + mask1;
868  else if ((inp > 0) && !(inp == mask3))
869    rnd = inp + mask1;
870  else
871    rnd = inp;
872
873  rnd = rnd >> sx;
874
875  if (inp == mask3) rnd++;
876
877  return rnd;
878}
879
880/*****************************************************************************
881  fixp_round()
882    round [typical rounding]
883
884    See fct roundRef() [which is the reference]
885  returns FIXP_DBL
886*****************************************************************************/
887FIXP_DBL fixp_round(FIXP_DBL f_inp, INT sf) {
888  FDK_ASSERT(sf >= 0);
889  INT sx = DFRACT_BITS - 1 - sf;
890  INT r = fixp_roundToInt(f_inp, sf);
891  ULONG mask = (ULONG)0x1 << (DFRACT_BITS - 1);  // 0x80000000
892  r = r << sx;
893
894  if ((f_inp > FL2FXCONST_DBL(0.0f)) && (r & mask)) {
895    --r;
896  }
897
898  FIXP_DBL f_round = (FIXP_DBL)r;
899  return f_round;
900}
901