Line data Source code
1 : /******************************************************************************************************
2 :
3 : (C) 2022-2025 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB,
4 : Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,
5 : Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,
6 : Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other
7 : contributors to this repository. All Rights Reserved.
8 :
9 : This software is protected by copyright law and by international treaties.
10 : The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB,
11 : Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,
12 : Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,
13 : Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other
14 : contributors to this repository retain full ownership rights in their respective contributions in
15 : the software. This notice grants no license of any kind, including but not limited to patent
16 : license, nor is any license granted by implication, estoppel or otherwise.
17 :
18 : Contributors are required to enter into the IVAS codec Public Collaboration agreement before making
19 : contributions.
20 :
21 : This software is provided "AS IS", without any express or implied warranties. The software is in the
22 : development stage. It is intended exclusively for experts who have experience with such software and
23 : solely for the purpose of inspection. All implied warranties of non-infringement, merchantability
24 : and fitness for a particular purpose are hereby disclaimed and excluded.
25 :
26 : Any dispute, controversy or claim arising under or in relation to providing this software shall be
27 : submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in
28 : accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and
29 : the United Nations Convention on Contracts on the International Sales of Goods.
30 :
31 : *******************************************************************************************************/
32 :
33 : /*====================================================================================
34 : EVS Codec 3GPP TS26.443 Nov 04, 2021. Version 12.14.0 / 13.10.0 / 14.6.0 / 15.4.0 / 16.3.0
35 : ====================================================================================*/
36 :
37 : #include <assert.h>
38 : #include <stdint.h>
39 : #include "options.h"
40 : #include "complex_basop.h"
41 : #include "basop_util.h"
42 : #include "basop32.h"
43 : #include "rom_com.h"
44 : #include "rom_basic_math.h"
45 : #include "basop_settings.h"
46 : #include "cnst.h"
47 :
48 : extern const Word32 SqrtTable[32]; // Q31
49 : extern const Word16 SqrtDiffTable[32]; /* Q15 */
50 :
51 : extern const Word32 ISqrtTable[32];
52 : extern const Word16 ISqrtDiffTable[32];
53 :
54 : extern const Word32 InvTable[32];
55 : extern const Word16 InvDiffTable[32];
56 :
57 13372015 : Word32 BASOP_Util_Log2( Word32 x )
58 : {
59 : Word32 exp;
60 : Word16 exp_e;
61 : Word16 nIn;
62 : Word16 accuSqr;
63 : Word32 accuRes;
64 :
65 :
66 13372015 : assert( x >= 0 );
67 :
68 13372015 : if ( x == 0 )
69 : {
70 :
71 42416 : return ( (Word32) MIN_32 );
72 : }
73 :
74 : /* normalize input, calculate integer part */
75 13329599 : exp_e = norm_l( x );
76 13329599 : x = L_shl( x, exp_e );
77 13329599 : exp = L_deposit_l( exp_e );
78 :
79 : /* calculate (1-normalized_input) */
80 13329599 : nIn = extract_h( L_sub( MAX_32, x ) );
81 :
82 : /* approximate ln() for fractional part (nIn *c0 + nIn^2*c1 + nIn^3*c2 + ... + nIn^8 *c7) */
83 :
84 : /* iteration 1, no need for accumulation */
85 13329599 : accuRes = L_mult( nIn, ldCoeff[0] ); /* nIn^i * coeff[0] */
86 13329599 : accuSqr = mult( nIn, nIn ); /* nIn^2, nIn^3 .... */
87 :
88 : /* iteration 2 */
89 13329599 : accuRes = L_mac( accuRes, accuSqr, ldCoeff[1] ); /* nIn^i * coeff[1] */
90 13329599 : accuSqr = mult( accuSqr, nIn ); /* nIn^2, nIn^3 .... */
91 :
92 : /* iteration 3 */
93 13329599 : accuRes = L_mac( accuRes, accuSqr, ldCoeff[2] ); /* nIn^i * coeff[2] */
94 13329599 : accuSqr = mult( accuSqr, nIn ); /* nIn^2, nIn^3 .... */
95 :
96 : /* iteration 4 */
97 13329599 : accuRes = L_mac( accuRes, accuSqr, ldCoeff[3] ); /* nIn^i * coeff[3] */
98 13329599 : accuSqr = mult( accuSqr, nIn ); /* nIn^2, nIn^3 .... */
99 :
100 : /* iteration 5 */
101 13329599 : accuRes = L_mac( accuRes, accuSqr, ldCoeff[4] ); /* nIn^i * coeff[4] */
102 13329599 : accuSqr = mult( accuSqr, nIn ); /* nIn^2, nIn^3 .... */
103 :
104 : /* iteration 6 */
105 13329599 : accuRes = L_mac( accuRes, accuSqr, ldCoeff[5] ); /* nIn^i * coeff[5] */
106 13329599 : accuSqr = mult( accuSqr, nIn ); /* nIn^2, nIn^3 .... */
107 :
108 : /* iteration 7, no need to calculate accuSqr any more */
109 13329599 : accuRes = L_mac( accuRes, accuSqr, ldCoeff[6] ); /* nIn^i * coeff[6] */
110 :
111 : /* ld(fractional part) = ln(fractional part)/ln(2), 1/ln(2) = (1 + 0.44269504) */
112 13329599 : accuRes = L_mac0( L_shr( accuRes, 1 ), extract_h( accuRes ), 14506 );
113 :
114 13329599 : accuRes = L_shr( accuRes, LD_DATA_SCALE - 1 ); /* fractional part/LD_DATA_SCALE */
115 13329599 : exp = L_shl( exp, ( 31 - LD_DATA_SCALE ) ); /* integer part/LD_DATA_SCALE */
116 13329599 : accuRes = L_sub( accuRes, exp ); /* result = integer part + fractional part */
117 :
118 :
119 13329599 : return ( accuRes );
120 : }
121 :
122 7367089 : Word32 BASOP_Util_InvLog2( Word32 x )
123 : {
124 : Word16 frac;
125 : Word16 exp;
126 : Word32 retVal;
127 : UWord32 index3;
128 : UWord32 index2;
129 : UWord32 index1;
130 : UWord32 lookup3f;
131 : UWord32 lookup12;
132 : UWord32 lookup;
133 :
134 :
135 7367089 : if ( x < -1040187392l /*-31.0/64.0 Q31*/ )
136 : {
137 :
138 0 : return 0;
139 : }
140 7367089 : test();
141 7367089 : if ( ( GE_32( x, 1040187392l /*31.0/64.0 Q31*/ ) ) || ( x == 0 ) )
142 : {
143 :
144 0 : return 0x7FFFFFFF;
145 : }
146 :
147 7367089 : frac = extract_l( L_and( x, 0x3FF ) );
148 :
149 7367089 : index3 = L_and( L_shr( x, 10 ), 0x1F );
150 7367089 : index2 = L_and( L_shr( x, 15 ), 0x1F );
151 7367089 : index1 = L_and( L_shr( x, 20 ), 0x1F );
152 :
153 7367089 : exp = extract_l( L_shr( x, 25 ) );
154 7367089 : if ( x > 0 )
155 : {
156 160 : exp = sub( 31, exp );
157 : }
158 7367089 : if ( x < 0 )
159 : {
160 7366929 : exp = negate( exp );
161 : }
162 :
163 7367089 : lookup3f = L_add( exp2x_tab_long[index3], L_shr( Mpy_32_16_1( 0x0016302F, frac ), 1 ) );
164 7367089 : lookup12 = Mpy_32_32( exp2_tab_long[index1], exp2w_tab_long[index2] );
165 7367089 : lookup = Mpy_32_32( lookup12, lookup3f );
166 :
167 7367089 : retVal = L_shr( lookup, sub( exp, 3 ) );
168 :
169 :
170 7367089 : return retVal;
171 : }
172 :
173 0 : Word32 BASOP_Util_Log10( Word32 x, Word16 e )
174 : {
175 0 : test();
176 0 : IF( e >= 0 && LE_16( e, 31 ) )
177 : {
178 0 : IF( EQ_32( x, L_shl_sat( 1, sub( 31, e ) ) ) )
179 : {
180 0 : return 0;
181 : }
182 : }
183 0 : Word32 res = BASOP_Util_Log2( x );
184 0 : res = L_add( Mpy_32_32( res, 646456993 /* log10(2) in Q31 */ ), Mpy_32_32( L_shl( e, 24 ), 1292913986 /* log10(2) in Q32 */ ) ); // Adjusting for the exponent mismatch: multiplying first so as to avoid saturation
185 : /* log10(2) is used in Q32 to support exponent till 127 in Mpy_32_32( L_shl( e, 24 ), 1292913986 )*/
186 0 : return res;
187 : }
188 :
189 0 : Word32 BASOP_Util_Loge( Word32 x, Word16 e )
190 : {
191 0 : Word32 res = BASOP_Util_Log2( x );
192 0 : res = L_add( Mpy_32_32( res, 1488522235 /* loge(2) in Q31 */ ), Mpy_32_32( L_shl( e, 25 ), 1488522235 /* loge(2) in Q31 */ ) ); // Adjusting for the exponent mismatch: multiplying first so as to avoid saturation
193 0 : return res;
194 : }
195 :
196 33631873 : Word16 BASOP_Util_Add_MantExp /*!< Exponent of result */
197 : ( Word16 a_m, /*!< Mantissa of 1st operand a */
198 : Word16 a_e, /*!< Exponent of 1st operand a */
199 : Word16 b_m, /*!< Mantissa of 2nd operand b */
200 : Word16 b_e, /*!< Exponent of 2nd operand b */
201 : Word16 *ptrSum_m ) /*!< Mantissa of result */
202 : {
203 : Word32 L_lm, L_hm;
204 : Word16 shift;
205 :
206 :
207 : /* Compare exponents: the difference is limited to +/- 15
208 : The Word16 mantissa of the operand with higher exponent is moved into the low
209 : part of a Word32 and shifted left by the exponent difference. Then, the
210 : unshifted mantissa of the operand with the lower exponent is added to the lower
211 : 16 bits. The addition result is normalized and the upper Word16 of the result represents
212 : the mantissa to return. The returned exponent takes into account all shift operations
213 : including the final 16-bit extraction.
214 : Note: The resulting mantissa may be inaccurate in the case, where the mantissa of the operand
215 : with higher exponent is not really left-aligned, while the mantissa of the operand with
216 : lower exponent is so. If in such a case, the difference in exponents is more than 15,
217 : an inaccuracy is introduced.
218 : Example:
219 : A: a_e = 20, a_m = 0x0001
220 : B: b_e = 0, b_m = 0x4000
221 : correct: A+B=1*2^20+1*2^14=0x0010.0000+0x0000.4000=0x0010.4000=0x4100*2^6
222 : previously: A+B=1*2^20+1*2^14=0x0001+0x0000=0x0001*2^20
223 : this version: A+B=1*2^20+1*2^14=0x0000.8000+0x0000.4000=0x6000*2^6
224 : */
225 :
226 33631873 : shift = sub( a_e, b_e );
227 33631873 : if ( shift >= 0 )
228 18687733 : shift = s_min( 15, shift );
229 :
230 33631873 : if ( shift < 0 )
231 14944140 : shift = s_max( -15, shift );
232 33631873 : a_e = s_max( a_e, b_e );
233 33631873 : L_hm = L_deposit_l( a_m ); /* mantissa belonging to higher exponent */
234 33631873 : L_lm = L_deposit_l( a_m ); /* mantissa belonging to lower exponent */
235 33631873 : if ( shift >= 0 )
236 18687733 : L_lm = L_deposit_l( b_m );
237 33631873 : if ( shift < 0 )
238 14944140 : L_hm = L_deposit_l( b_m );
239 :
240 33631873 : if ( shift > 0 )
241 16044420 : shift = negate( shift );
242 :
243 33631873 : L_hm = L_shr( L_hm, shift ); /* shift left due to negative shift parameter */
244 33631873 : a_e = add( a_e, shift );
245 33631873 : L_hm = L_add( L_hm, L_lm );
246 33631873 : shift = norm_l( L_hm );
247 33631873 : L_hm = L_shl( L_hm, shift );
248 33631873 : *ptrSum_m = extract_h( L_hm );
249 33631873 : move16();
250 :
251 33631873 : a_e = sub( a_e, shift );
252 33631873 : if ( L_hm )
253 33393065 : a_e = add( a_e, 16 );
254 33631873 : return ( a_e );
255 : }
256 :
257 :
258 547163 : void BASOP_Util_Divide_MantExp( Word16 a_m, /*!< Mantissa of dividend a */
259 : Word16 a_e, /*!< Exponent of dividend a */
260 : Word16 b_m, /*!< Mantissa of divisor b */
261 : Word16 b_e, /*!< Exponent of divisor b */
262 : Word16 *ptrResult_m, /*!< Mantissa of quotient a/b */
263 : Word16 *ptrResult_e /*!< Exponent of quotient a/b */
264 : )
265 : {
266 : Word16 index, frac;
267 : Word16 preShift, postShift;
268 : Word16 m;
269 : Word32 m32;
270 : #ifdef BASOP_NOGLOB_DECLARE_LOCAL
271 547163 : Flag Overflow = 0;
272 : #endif
273 :
274 :
275 547163 : assert( b_m != 0 );
276 :
277 : /* normalize b */
278 547163 : preShift = norm_s( b_m );
279 547163 : m = shl( b_m, preShift );
280 :
281 : /* make b positive */
282 : BASOP_SATURATE_WARNING_OFF_EVS;
283 547163 : if ( m < 0 )
284 0 : m = negate( m );
285 : BASOP_SATURATE_WARNING_ON_EVS;
286 :
287 : /* get table index (upper 6 bits minus 16) */
288 : /* index = (m >> 9) - 32; */
289 547163 : index = mac_r( -32768 - ( 32 << 16 ), m, 1 << 6 );
290 :
291 : /* get fractional part for interpolation (lower 9 bits) */
292 547163 : frac = shl( s_and( m, 0x1FF ), 1 ); /* Q10 */
293 :
294 : /* interpolate 1/b */
295 547163 : m = msu_r( InvTable[index], InvDiffTable[index], frac );
296 :
297 : /* restore sign */
298 547163 : if ( b_m < 0 )
299 0 : m = negate( m );
300 :
301 : /* multiply with a */
302 547163 : m32 = L_mult( a_m, m );
303 :
304 : /* normalize result */
305 547163 : postShift = norm_l( m32 );
306 547163 : m = round_fx_o( L_shl( m32, postShift ), &Overflow );
307 :
308 : /* exponent */
309 547163 : *ptrResult_e = sub( add( add( a_e, sub( 1, b_e ) ), preShift ), postShift );
310 547163 : move16();
311 :
312 547163 : *ptrResult_m = m;
313 547163 : move16();
314 547163 : }
315 :
316 :
317 : /* local function for Sqrt16 */
318 76543394 : static Word16 Sqrt16_common( Word16 m,
319 : Word16 e )
320 : {
321 : Word16 index, frac;
322 : Flag Overflow;
323 :
324 :
325 76543394 : assert( ( m >= 0x4000 ) || ( m == 0 ) );
326 :
327 : /* get table index (upper 6 bits minus 32) */
328 : /* index = (m >> 9) - 32; */
329 76543394 : index = mac_r( -32768 - ( 32 << 16 ), m, 1 << 6 );
330 :
331 : /* get fractional part for interpolation (lower 9 bits) */
332 76543394 : frac = s_and( m, 0x1FF ); /* Q9 */
333 :
334 : /* interpolate */
335 76543394 : if ( m != 0 )
336 : {
337 : BASOP_SATURATE_WARNING_OFF_EVS;
338 75429485 : m = mac_ro( SqrtTable[index], SqrtDiffTable[index], frac, &Overflow );
339 : BASOP_SATURATE_WARNING_ON_EVS;
340 : }
341 :
342 : /* handle odd exponents */
343 76543394 : if ( s_and( e, 1 ) != 0 )
344 34095336 : m = mult_r( m, 0x5a82 );
345 :
346 76543394 : return m;
347 : }
348 :
349 :
350 : /* local function for Sqrt32 and Sqrt32norm */
351 691075412 : static Word32 Sqrt32_common( Word32 m,
352 : Word16 e )
353 : {
354 : Word16 m16, index, frac;
355 : #ifdef BASOP_NOGLOB_DECLARE_LOCAL
356 691075412 : Flag Overflow = 0;
357 : #endif
358 :
359 691075412 : assert( ( m >= 0x40000000 ) || ( m == 0 ) );
360 :
361 691075412 : m16 = round_fx_o( m, &Overflow );
362 :
363 : /* get table index (upper 6 bits minus 32) */
364 : /* index = (m16 >> 9) - 32; */
365 691075412 : index = mac_r( -32768 - ( 32 << 16 ), m16, 1 << 6 );
366 :
367 : /* get fractional part for interpolation (lower 9 bits) */
368 691075412 : frac = s_and( m16, 0x1FF ); /* Q9 */
369 :
370 : /* interpolate */
371 691075412 : if ( m != 0 )
372 : {
373 : BASOP_SATURATE_WARNING_OFF_EVS;
374 644238551 : m = L_mac_sat( SqrtTable[index], SqrtDiffTable[index], frac );
375 : BASOP_SATURATE_WARNING_ON_EVS;
376 : }
377 :
378 : /* handle odd exponents */
379 691075412 : if ( s_and( e, 1 ) != 0 )
380 328877020 : m = Mpy_32_16_1( m, 0x5a82 );
381 :
382 691075412 : return m;
383 : }
384 :
385 : /* local function for ISqrt16 and ISqrt16norm */
386 39199670 : static Word16 ISqrt16_common( Word16 m,
387 : Word16 e )
388 : {
389 : Word16 index, frac;
390 :
391 39199670 : assert( m >= 0x4000 );
392 :
393 : /* get table index (upper 6 bits minus 32) */
394 : /* index = (m >> 9) - 32; */
395 39199670 : index = mac_r( -32768 - ( 32 << 16 ), m, 1 << 6 );
396 :
397 : /* get fractional part for interpolation (lower 9 bits) */
398 39199670 : frac = s_and( m, 0x1FF ); /* Q9 */
399 :
400 : /* interpolate */
401 39199670 : m = msu_r( ISqrtTable[index], ISqrtDiffTable[index], frac );
402 :
403 : /* handle even exponents */
404 39199670 : if ( s_and( e, 1 ) == 0 )
405 : {
406 17253465 : m = mult_r( m, 0x5a82 );
407 : }
408 :
409 39199670 : return m;
410 : }
411 :
412 : /* local function for ISqrt32 and ISqrt32norm */
413 66972872 : static Word32 ISqrt32_common( Word32 m,
414 : Word16 e )
415 : {
416 : Word16 m16, index, frac;
417 : #ifdef BASOP_NOGLOB_DECLARE_LOCAL
418 66972872 : Flag Overflow = 0;
419 : #endif
420 :
421 66972872 : assert( m >= 0x40000000 );
422 : #ifdef BASOP_NOGLOB_DECLARE_LOCAL
423 66972872 : m16 = round_fx_o( m, &Overflow );
424 : #else
425 : m16 = round_fx( m );
426 : #endif
427 :
428 : /* get table index (upper 6 bits minus 32) */
429 : /* index = (m16 >> 25) - 32; */
430 66972872 : index = mac_r( -32768 - ( 32 << 16 ), m16, 1 << 6 );
431 :
432 : /* get fractional part for interpolation (lower 9 bits) */
433 66972872 : frac = s_and( m16, 0x1FF ); /* Q9 */
434 :
435 : /* interpolate */
436 66972872 : m = L_msu( ISqrtTable[index], ISqrtDiffTable[index], frac );
437 :
438 : /* handle even exponents */
439 66972872 : if ( s_and( e, 1 ) == 0 )
440 : {
441 38503015 : m = Mpy_32_16_1( m, 0x5a82 );
442 : }
443 :
444 66972872 : return m;
445 : }
446 :
447 76003967 : Word16 Sqrt16( /*!< output mantissa */
448 : Word16 mantissa, /*!< input mantissa */
449 : Word16 *exponent /*!< pointer to exponent */
450 : )
451 : {
452 : Word16 preShift, e;
453 :
454 76003967 : assert( mantissa >= 0 );
455 :
456 : /* normalize */
457 76003967 : preShift = norm_s( mantissa );
458 :
459 76003967 : e = sub( *exponent, preShift );
460 76003967 : mantissa = shl( mantissa, preShift );
461 :
462 : /* calc mantissa */
463 76003967 : mantissa = Sqrt16_common( mantissa, e );
464 :
465 : /* e = (e + 1) >> 1 */
466 76003967 : *exponent = mult_r( e, 1 << 14 );
467 76003967 : move16();
468 :
469 76003967 : return mantissa;
470 : }
471 :
472 :
473 539427 : Word16 Sqrt16norm( /*!< output mantissa */
474 : Word16 mantissa, /*!< normalized input mantissa */
475 : Word16 *exponent /*!< pointer to exponent */
476 : )
477 : {
478 :
479 539427 : assert( ( mantissa >= 0x4000 ) || ( mantissa == 0 ) );
480 :
481 : /* calc mantissa */
482 539427 : mantissa = Sqrt16_common( mantissa, *exponent );
483 :
484 : /* e = (e + 1) >> 1 */
485 539427 : *exponent = mult_r( *exponent, 1 << 14 );
486 539427 : move16();
487 :
488 539427 : return mantissa;
489 : }
490 :
491 :
492 39199670 : Word16 ISqrt16( /*!< output mantissa */
493 : Word16 mantissa, /*!< input mantissa */
494 : Word16 *exponent /*!< pointer to exponent */
495 : )
496 : {
497 : Word16 preShift, e;
498 :
499 39199670 : assert( mantissa > 0 );
500 :
501 : /* normalize */
502 39199670 : preShift = norm_s( mantissa );
503 :
504 39199670 : e = sub( *exponent, preShift );
505 39199670 : mantissa = shl( mantissa, preShift );
506 :
507 : /* calc mantissa */
508 39199670 : mantissa = ISqrt16_common( mantissa, e );
509 :
510 : /* e = (2 - e) >> 1 */
511 39199670 : *exponent = msu_r( 1L << 15, e, 1 << 14 );
512 39199670 : move16();
513 :
514 39199670 : return mantissa;
515 : }
516 :
517 691075412 : Word32 Sqrt32( /*!< output mantissa */
518 : Word32 mantissa, /*!< input mantissa */
519 : Word16 *exponent /*!< pointer to exponent */
520 : )
521 : {
522 : Word16 preShift, e;
523 :
524 691075412 : assert( mantissa >= 0 );
525 :
526 : /* normalize */
527 691075412 : preShift = norm_l( mantissa );
528 :
529 691075412 : e = sub( *exponent, preShift );
530 691075412 : mantissa = L_shl( mantissa, preShift );
531 :
532 : /* calc mantissa */
533 691075412 : mantissa = Sqrt32_common( mantissa, e );
534 :
535 : /* e = (e + 1) >> 1 */
536 691075412 : *exponent = mult_r( e, 1 << 14 );
537 691075412 : move16();
538 :
539 691075412 : return mantissa;
540 : }
541 :
542 :
543 66958830 : Word32 ISqrt32( /*!< output mantissa */
544 : Word32 mantissa, /*!< input mantissa */
545 : Word16 *exponent /*!< pointer to exponent */
546 : )
547 : {
548 : Word16 preShift, e;
549 :
550 66958830 : assert( mantissa > 0 );
551 :
552 : /* normalize */
553 66958830 : preShift = norm_l( mantissa );
554 :
555 66958830 : e = sub( *exponent, preShift );
556 66958830 : mantissa = L_shl( mantissa, preShift );
557 :
558 : /* calc mantissa */
559 66958830 : mantissa = ISqrt32_common( mantissa, e );
560 :
561 : /* e = (2 - e) >> 1 */
562 66958830 : *exponent = msu_r( 1L << 15, e, 1 << 14 );
563 66958830 : move16();
564 :
565 66958830 : return mantissa;
566 : }
567 :
568 14042 : Word32 ISqrt32norm( /*!< output mantissa */
569 : Word32 mantissa, /*!< normalized input mantissa */
570 : Word16 *exponent /*!< pointer to exponent */
571 : )
572 : {
573 :
574 14042 : assert( mantissa >= 0x40000000 );
575 :
576 : /* calc mantissa */
577 14042 : mantissa = ISqrt32_common( mantissa, *exponent );
578 :
579 : /* e = (2 - e) >> 1 */
580 14042 : *exponent = msu_r( 1L << 15, *exponent, 1 << 14 );
581 14042 : move16();
582 :
583 14042 : return mantissa;
584 : }
585 :
586 5605044 : Word16 Inv16( /*!< output mantissa */
587 : Word16 mantissa, /*!< input mantissa */
588 : Word16 *exponent /*!< pointer to exponent */
589 : )
590 : {
591 : Word16 index, frac;
592 : Word16 preShift;
593 : Word16 m, e;
594 :
595 :
596 5605044 : assert( mantissa != 0 );
597 :
598 : /* absolute */
599 : BASOP_SATURATE_WARNING_OFF_EVS;
600 5605044 : m = abs_s( mantissa );
601 : BASOP_SATURATE_WARNING_ON_EVS;
602 :
603 : /* normalize */
604 5605044 : preShift = norm_s( m );
605 :
606 5605044 : e = sub( *exponent, preShift );
607 5605044 : m = shl( m, preShift );
608 :
609 : /* get table index (upper 6 bits minus 32) */
610 : /* index = (m >> 9) - 32; */
611 5605044 : index = mac_r( -32768 - ( 32 << 16 ), m, 1 << 6 );
612 :
613 : /* get fractional part for interpolation (lower 9 bits) */
614 5605044 : frac = shl( s_and( m, 0x1FF ), 1 ); /* Q10 */
615 :
616 : /* interpolate */
617 5605044 : m = msu_r( InvTable[index], InvDiffTable[index], frac );
618 :
619 : /* restore sign */
620 5605044 : if ( mantissa < 0 )
621 2342510 : m = negate( m );
622 :
623 : /* e = 1 - e */
624 5605044 : *exponent = sub( 1, e );
625 5605044 : move16();
626 :
627 5605044 : return m;
628 : }
629 :
630 :
631 15662656 : void BASOP_Util_Sqrt_InvSqrt_MantExp( Word16 mantissa, /*!< mantissa */
632 : Word16 exponent, /*!< expoinent */
633 : Word16 *sqrt_mant, /*!< Pointer to sqrt mantissa */
634 : Word16 *sqrt_exp, /*!< Pointer to sqrt exponent */
635 : Word16 *isqrt_mant, /*!< Pointer to 1/sqrt mantissa */
636 : Word16 *isqrt_exp /*!< Pointer to 1/sqrt exponent */
637 : )
638 : {
639 : Word16 index, frac;
640 : Word16 preShift;
641 : Word16 m, mi, e_odd;
642 :
643 :
644 15662656 : assert( mantissa > 0 );
645 :
646 : /* normalize */
647 15662656 : preShift = norm_s( mantissa );
648 :
649 15662656 : exponent = sub( exponent, preShift );
650 15662656 : mantissa = shl( mantissa, preShift );
651 :
652 : /* get table index (upper 6 bits minus 32) */
653 : /* index = (m >> 9) - 32; */
654 15662656 : index = mac_r( -32768 - ( 32 << 16 ), mantissa, 1 << 6 );
655 :
656 : /* get fractional part for interpolation (lower 9 bits) */
657 15662656 : frac = s_and( mantissa, 0x1FF ); /* Q9 */
658 :
659 : /* interpolate */
660 : BASOP_SATURATE_WARNING_OFF_EVS;
661 15662656 : m = mac_r_sat( SqrtTable[index], SqrtDiffTable[index], frac );
662 15662656 : mi = msu_r_sat( ISqrtTable[index], ISqrtDiffTable[index], frac );
663 : BASOP_SATURATE_WARNING_ON_EVS;
664 :
665 : /* handle even/odd exponents */
666 15662656 : e_odd = s_and( exponent, 1 );
667 15662656 : if ( e_odd != 0 )
668 7923627 : m = mult_r( m, 0x5a82 );
669 15662656 : if ( e_odd == 0 )
670 7739029 : mi = mult_r( mi, 0x5a82 );
671 :
672 : /* e = (e + 1) >> 1 */
673 15662656 : *sqrt_exp = mult_r( exponent, 1 << 14 );
674 15662656 : move16();
675 :
676 : /* e = (2 - e) >> 1 */
677 15662656 : *isqrt_exp = msu_r( 1L << 15, exponent, 1 << 14 );
678 15662656 : move16();
679 :
680 : /* Write result */
681 15662656 : *sqrt_mant = m;
682 15662656 : move16();
683 15662656 : *isqrt_mant = mi;
684 15662656 : move16();
685 15662656 : }
686 :
687 :
688 : /********************************************************************/
689 : /*!
690 : \brief Calculates the scalefactor needed to normalize input array
691 :
692 : The scalefactor needed to normalize the Word16 input array is returned <br>
693 : If the input array contains only '0', a scalefactor 0 is returned <br>
694 : Scaling factor is determined wrt a normalized target x: 16384 <= x <= 32767 for positive x <br>
695 : and -32768 <= x <= -16384 for negative x
696 : */
697 :
698 4422327 : Word16 getScaleFactor16( /* o: measured headroom in range [0..15], 0 if all x[i] == 0 */
699 : const Word16 *x, /* i: array containing 16-bit data */
700 : const Word16 len_x ) /* i: length of the array to scan */
701 : {
702 : Word16 i, i_min, i_max;
703 : Word16 x_min, x_max;
704 :
705 :
706 4422327 : x_max = 0;
707 4422327 : move16();
708 4422327 : x_min = 0;
709 4422327 : move16();
710 5592398681 : FOR( i = 0; i < len_x; i++ )
711 : {
712 5587976354 : if ( x[i] >= 0 )
713 3477872200 : x_max = s_max( x_max, x[i] );
714 5587976354 : if ( x[i] < 0 )
715 2110104154 : x_min = s_min( x_min, x[i] );
716 : }
717 :
718 4422327 : i_max = 0x10;
719 4422327 : move16();
720 4422327 : i_min = 0x10;
721 4422327 : move16();
722 :
723 4422327 : if ( x_max != 0 )
724 3918276 : i_max = norm_s( x_max );
725 :
726 4422327 : if ( x_min != 0 )
727 3868337 : i_min = norm_s( x_min );
728 :
729 4422327 : i = s_and( s_min( i_max, i_min ), 0xF );
730 :
731 :
732 4422327 : return i;
733 : }
734 :
735 :
736 : /********************************************************************/
737 : /*!
738 : \brief Calculates the scalefactor needed to normalize input array
739 :
740 : The scalefactor needed to normalize the Word32 input array is returned <br>
741 : If the input array contains only '0', a scalefactor 0 is returned <br>
742 : Scaling factor is determined wrt a normalized target x: 1073741824 <= x <= 2147483647 for positive x <br>
743 : and -2147483648 <= x <= -1073741824 for negative x
744 : */
745 :
746 91187833 : Word16 getScaleFactor32( /* o: measured headroom in range [0..31], 0 if all x[i] == 0 */
747 : const Word32 *x, /* i: array containing 32-bit data */
748 : const Word16 len_x ) /* i: length of the array to scan */
749 : {
750 : Word16 i, i_min, i_max;
751 : Word32 x_min, x_max;
752 :
753 :
754 91187833 : x_max = 0;
755 91187833 : move32();
756 91187833 : x_min = 0;
757 91187833 : move32();
758 9511545955 : FOR( i = 0; i < len_x; i++ )
759 : {
760 9420358122 : if ( x[i] >= 0 )
761 6294291380 : x_max = L_max( x_max, x[i] );
762 9420358122 : if ( x[i] < 0 )
763 3126066742 : x_min = L_min( x_min, x[i] );
764 : }
765 :
766 91187833 : i_max = 0x20;
767 91187833 : move16();
768 91187833 : i_min = 0x20;
769 91187833 : move16();
770 :
771 91187833 : if ( x_max != 0 )
772 84486026 : i_max = norm_l( x_max );
773 :
774 91187833 : if ( x_min != 0 )
775 78308332 : i_min = norm_l( x_min );
776 :
777 91187833 : i = s_and( s_min( i_max, i_min ), 0x1F );
778 :
779 :
780 91187833 : return i;
781 : }
782 :
783 19405 : Word16 getScaleFactor32_copy( /* o: measured headroom in range [0..31], 0 if all x[i] == 0 */
784 : const Word32 *x, /* i: array containing 32-bit data */
785 : const Word32 len_x ) /* i: length of the array to scan */
786 : {
787 : Word32 i;
788 : Word16 i_min, i_max;
789 : Word32 x_min, x_max;
790 :
791 :
792 19405 : x_max = 0;
793 19405 : move32();
794 19405 : x_min = 0;
795 19405 : move32();
796 125907085 : FOR( i = 0; i < len_x; i++ )
797 : {
798 125887680 : if ( x[i] >= 0 )
799 105057578 : x_max = L_max( x_max, x[i] );
800 125887680 : if ( x[i] < 0 )
801 20830102 : x_min = L_min( x_min, x[i] );
802 : }
803 :
804 19405 : i_max = 0x20;
805 19405 : move16();
806 19405 : i_min = 0x20;
807 19405 : move16();
808 :
809 19405 : if ( x_max != 0 )
810 19399 : i_max = norm_l( x_max );
811 :
812 19405 : if ( x_min != 0 )
813 19405 : i_min = norm_l( x_min );
814 :
815 19405 : i_max = s_and( s_min( i_max, i_min ), 0x1F );
816 :
817 :
818 19405 : return i_max;
819 : }
820 :
821 0 : Word16 normalize16( Word16 mantissa, Word16 *pexponent )
822 : {
823 : Word16 tmp;
824 :
825 0 : tmp = norm_s( mantissa );
826 0 : mantissa = shl( mantissa, tmp );
827 0 : move16();
828 0 : *pexponent = sub( *pexponent, tmp );
829 :
830 :
831 0 : return mantissa;
832 : }
833 159419 : Word16 divide3216( Word32 x, Word16 y )
834 : {
835 : Word16 z;
836 :
837 :
838 159419 : z = 0;
839 159419 : move16();
840 159419 : if ( 0 == y )
841 : {
842 0 : return 0x7fff;
843 : }
844 :
845 159419 : IF( x != 0 )
846 : {
847 : Word16 den, sign;
848 : Word32 num;
849 159419 : num = L_abs( x );
850 159419 : den = abs_s( y );
851 :
852 159419 : sign = s_and( s_xor( extract_h( x ), y ), -32768 /* 0x8000 */ );
853 :
854 159419 : z = div_l( num, den );
855 159419 : if ( 0 != sign )
856 : {
857 0 : z = negate( z );
858 : }
859 : }
860 :
861 :
862 159419 : return z;
863 : }
864 :
865 78240 : Word16 divide1616( Word16 x, Word16 y )
866 : {
867 : Word16 z, num, den, sign;
868 :
869 :
870 78240 : num = abs_s( x );
871 78240 : den = abs_s( y );
872 :
873 78240 : sign = s_and( s_xor( x, y ), -32768 /* 0x8000 */ );
874 :
875 78240 : move16();
876 78240 : z = 0x7fff;
877 78240 : if ( LT_16( num, den ) )
878 78240 : z = div_s( num, den );
879 :
880 78240 : if ( 0 != sign )
881 : {
882 10944 : z = negate( z );
883 : }
884 :
885 :
886 78240 : return z;
887 : }
888 :
889 4131335 : Word16 divide3232( Word32 L_num, Word32 L_denom )
890 : {
891 : Word16 z;
892 : Word32 sign;
893 : #ifdef BASOP_NOGLOB_DECLARE_LOCAL
894 4131335 : Flag Overflow = 0;
895 : #endif
896 :
897 :
898 4131335 : sign = L_and( L_xor( L_num, L_denom ), (Word32) 0x80000000 );
899 :
900 4131335 : L_num = L_abs( L_num );
901 4131335 : L_denom = L_abs( L_denom );
902 :
903 : /* limit the range of denominator to Word16 */
904 4131335 : z = s_min( norm_l( L_num ), norm_l( L_denom ) );
905 4131335 : L_num = L_shl( L_num, z );
906 4131335 : L_denom = L_shl( L_denom, z );
907 :
908 : /* round_fx instead of extract_h improves spectral distortion in E_UTIL_lev_dur (schur version). */
909 4131335 : z = div_l( L_num, round_fx_o( L_denom, &Overflow ) );
910 4131335 : if ( 0 != sign )
911 : {
912 929674 : z = negate( z );
913 : }
914 :
915 :
916 4131335 : return z;
917 : }
918 :
919 42883461 : Word16 BASOP_Util_Divide3232_uu_1616_Scale( Word32 x, Word32 y, Word16 *s )
920 : {
921 : Word16 z;
922 : Word16 sx;
923 : Word16 sy;
924 : Word16 x16;
925 : Word16 y16;
926 :
927 :
928 42883461 : assert( x >= 0 );
929 42883461 : assert( y > 0 );
930 :
931 42883461 : if ( x == 0 )
932 : {
933 3947721 : *s = 0;
934 3947721 : move16();
935 :
936 :
937 3947721 : return ( 0 );
938 : }
939 :
940 38935740 : sx = norm_l( x );
941 38935740 : sy = norm_l( y );
942 :
943 38935740 : x16 = extract_h( L_shl( x, sx ) );
944 38935740 : y16 = extract_h( L_shl( y, sy ) );
945 :
946 38935740 : if ( GT_16( x16, y16 ) )
947 : {
948 13110212 : sx = sub( sx, 1 );
949 : }
950 :
951 38935740 : if ( LT_16( y16, x16 ) )
952 : {
953 13110212 : x16 = mult_r( x16, 0x4000 );
954 : }
955 :
956 :
957 38935740 : z = div_s( x16, y16 );
958 38935740 : move16();
959 38935740 : *s = sub( sy, sx );
960 :
961 :
962 38935740 : return ( z );
963 : }
964 :
965 0 : Word32 div_w( Word32 L_num, Word32 L_den )
966 : {
967 0 : Word32 L_var_out = 0;
968 : Word16 iteration;
969 0 : move32();
970 :
971 :
972 0 : IF( L_den == 0 )
973 : {
974 : /* printf("Division by 0 in div_l, Fatal error in "); printStack(); */
975 0 : return ( 0 );
976 : }
977 :
978 0 : test();
979 0 : IF( ( L_num < 0 ) || ( L_den < 0 ) )
980 : {
981 : /* printf("Division Error in div_l, Fatal error in "); printStack(); */
982 0 : return ( 0 );
983 : }
984 : Word64 W_num, W_den;
985 0 : W_num = W_deposit32_h( L_num );
986 0 : W_den = W_deposit32_h( L_den );
987 :
988 0 : IF( W_sub( W_num, W_den ) >= 0 )
989 : {
990 0 : return MAX_32;
991 : }
992 : ELSE
993 : {
994 0 : W_num = W_shr( W_num, 1 );
995 0 : W_den = W_shr( W_den, 1 );
996 :
997 0 : FOR( iteration = 0; iteration < 31; iteration++ )
998 : {
999 0 : L_var_out = L_shl( L_var_out, 1 );
1000 0 : W_num = W_shl( W_num, 1 );
1001 :
1002 0 : IF( W_sub( W_num, W_den ) >= 0 )
1003 : {
1004 0 : W_num = W_sub( W_num, W_den );
1005 0 : L_var_out = L_add( L_var_out, 1 );
1006 : }
1007 : }
1008 :
1009 0 : return L_var_out;
1010 : }
1011 : }
1012 :
1013 0 : Word32 BASOP_Util_Divide3232_Scale_cadence( Word32 x, Word32 y, Word16 *s )
1014 : {
1015 : Word32 z;
1016 : Word16 sx;
1017 : Word16 sy;
1018 : Word32 sign;
1019 :
1020 : /* assert (x >= (Word32)0); */
1021 0 : assert( y != (Word32) 0 );
1022 :
1023 0 : sign = 0;
1024 0 : move16();
1025 :
1026 0 : IF( x < 0 )
1027 : {
1028 0 : x = L_negate( x );
1029 0 : sign = L_xor( sign, 1 );
1030 : }
1031 :
1032 0 : IF( y < 0 )
1033 : {
1034 0 : y = L_negate( y );
1035 0 : sign = L_xor( sign, 1 );
1036 : }
1037 :
1038 0 : IF( x == (Word32) 0 )
1039 : {
1040 0 : *s = 0;
1041 0 : return ( (Word32) 0 );
1042 : }
1043 :
1044 0 : sx = norm_l( x );
1045 0 : x = L_shl( x, sx );
1046 0 : x = L_shr( x, 1 );
1047 0 : move16();
1048 0 : *s = sub( 1, sx );
1049 :
1050 0 : sy = norm_l( y );
1051 0 : y = L_shl( y, sy );
1052 0 : move16();
1053 0 : *s = add( *s, sy );
1054 :
1055 0 : z = div_w( x, y );
1056 :
1057 0 : if ( sign != 0 )
1058 : {
1059 0 : z = L_negate( z );
1060 : }
1061 :
1062 0 : return z;
1063 : }
1064 :
1065 : Word32 div_w_newton( Word32 num, Word32 den );
1066 : /*
1067 : Table of 256 precalculated estimates to be used by the "div_w_newton"
1068 : function using the Newton/Raphson method.
1069 : Note: The first table value (for denominator near 0x40000000) is not fully
1070 : accurate and should not be used.
1071 : */
1072 : Word32 div_w_newton_lookup[256] = {
1073 : /* Precalculated rounded results for 0x40000000 / b with b in [0x40000000 ... 0x7FFFFFFF] */
1074 : 0x7FFFFFFF, // 1.000000000000000 i=0 0.5 / 0.5+0/512 (b=0x40000000)
1075 : 0x7F807F80, // 0.996108949416342 i=1 0.5 / 0.5+1/512 (b=0x40400000)
1076 : 0x7F01FC07, // 0.992248062015504 i=2 0.5 / 0.5+2/512 (b=0x40800000)
1077 : 0x7E8472A8, // 0.988416988416988 i=3 0.5 / 0.5+3/512 (b=0x40C00000)
1078 : 0x7E07E07E, // 0.984615384615385 i=4 0.5 / 0.5+4/512 (b=0x41000000)
1079 : 0x7D8C42B2, // 0.980842911877395 i=5 0.5 / 0.5+5/512 (b=0x41400000)
1080 : 0x7D119679, // 0.977099236641221 i=6 0.5 / 0.5+6/512 (b=0x41800000)
1081 : 0x7C97D910, // 0.973384030418251 i=7 0.5 / 0.5+7/512 (b=0x41C00000)
1082 : 0x7C1F07C1, // 0.969696969696970 i=8 0.5 / 0.5+8/512 (b=0x42000000)
1083 : 0x7BA71FE1, // 0.966037735849057 i=9 0.5 / 0.5+9/512 (b=0x42400000)
1084 : 0x7B301ECC, // 0.962406015037594 i=10 0.5 / 0.5+10/512 (b=0x42800000)
1085 : 0x7ABA01EA, // 0.958801498127341 i=11 0.5 / 0.5+11/512 (b=0x42C00000)
1086 : 0x7A44C6AF, // 0.955223880597015 i=12 0.5 / 0.5+12/512 (b=0x43000000)
1087 : 0x79D06A96, // 0.951672862453532 i=13 0.5 / 0.5+13/512 (b=0x43400000)
1088 : 0x795CEB24, // 0.948148148148148 i=14 0.5 / 0.5+14/512 (b=0x43800000)
1089 : 0x78EA45E7, // 0.944649446494465 i=15 0.5 / 0.5+15/512 (b=0x43C00000)
1090 : 0x78787878, // 0.941176470588235 i=16 0.5 / 0.5+16/512 (b=0x44000000)
1091 : 0x78078078, // 0.937728937728938 i=17 0.5 / 0.5+17/512 (b=0x44400000)
1092 : 0x77975B8F, // 0.934306569343066 i=18 0.5 / 0.5+18/512 (b=0x44800000)
1093 : 0x77280772, // 0.930909090909091 i=19 0.5 / 0.5+19/512 (b=0x44C00000)
1094 : 0x76B981DA, // 0.927536231884058 i=20 0.5 / 0.5+20/512 (b=0x45000000)
1095 : 0x764BC88C, // 0.924187725631769 i=21 0.5 / 0.5+21/512 (b=0x45400000)
1096 : 0x75DED952, // 0.920863309352518 i=22 0.5 / 0.5+22/512 (b=0x45800000)
1097 : 0x7572B201, // 0.917562724014337 i=23 0.5 / 0.5+23/512 (b=0x45C00000)
1098 : 0x75075075, // 0.914285714285714 i=24 0.5 / 0.5+24/512 (b=0x46000000)
1099 : 0x749CB28F, // 0.911032028469751 i=25 0.5 / 0.5+25/512 (b=0x46400000)
1100 : 0x7432D63D, // 0.907801418439716 i=26 0.5 / 0.5+26/512 (b=0x46800000)
1101 : 0x73C9B971, // 0.904593639575972 i=27 0.5 / 0.5+27/512 (b=0x46C00000)
1102 : 0x73615A24, // 0.901408450704225 i=28 0.5 / 0.5+28/512 (b=0x47000000)
1103 : 0x72F9B658, // 0.898245614035088 i=29 0.5 / 0.5+29/512 (b=0x47400000)
1104 : 0x7292CC15, // 0.895104895104895 i=30 0.5 / 0.5+30/512 (b=0x47800000)
1105 : 0x722C996B, // 0.891986062717770 i=31 0.5 / 0.5+31/512 (b=0x47C00000)
1106 : 0x71C71C71, // 0.888888888888889 i=32 0.5 / 0.5+32/512 (b=0x48000000)
1107 : 0x71625344, // 0.885813148788927 i=33 0.5 / 0.5+33/512 (b=0x48400000)
1108 : 0x70FE3C07, // 0.882758620689655 i=34 0.5 / 0.5+34/512 (b=0x48800000)
1109 : 0x709AD4E4, // 0.879725085910653 i=35 0.5 / 0.5+35/512 (b=0x48C00000)
1110 : 0x70381C0E, // 0.876712328767123 i=36 0.5 / 0.5+36/512 (b=0x49000000)
1111 : 0x6FD60FBA, // 0.873720136518771 i=37 0.5 / 0.5+37/512 (b=0x49400000)
1112 : 0x6F74AE26, // 0.870748299319728 i=38 0.5 / 0.5+38/512 (b=0x49800000)
1113 : 0x6F13F596, // 0.867796610169492 i=39 0.5 / 0.5+39/512 (b=0x49C00000)
1114 : 0x6EB3E453, // 0.864864864864865 i=40 0.5 / 0.5+40/512 (b=0x4A000000)
1115 : 0x6E5478AC, // 0.861952861952862 i=41 0.5 / 0.5+41/512 (b=0x4A400000)
1116 : 0x6DF5B0F7, // 0.859060402684564 i=42 0.5 / 0.5+42/512 (b=0x4A800000)
1117 : 0x6D978B8E, // 0.856187290969900 i=43 0.5 / 0.5+43/512 (b=0x4AC00000)
1118 : 0x6D3A06D3, // 0.853333333333333 i=44 0.5 / 0.5+44/512 (b=0x4B000000)
1119 : 0x6CDD212B, // 0.850498338870432 i=45 0.5 / 0.5+45/512 (b=0x4B400000)
1120 : 0x6C80D901, // 0.847682119205298 i=46 0.5 / 0.5+46/512 (b=0x4B800000)
1121 : 0x6C252CC7, // 0.844884488448845 i=47 0.5 / 0.5+47/512 (b=0x4BC00000)
1122 : 0x6BCA1AF2, // 0.842105263157895 i=48 0.5 / 0.5+48/512 (b=0x4C000000)
1123 : 0x6B6FA1FE, // 0.839344262295082 i=49 0.5 / 0.5+49/512 (b=0x4C400000)
1124 : 0x6B15C06B, // 0.836601307189543 i=50 0.5 / 0.5+50/512 (b=0x4C800000)
1125 : 0x6ABC74BE, // 0.833876221498371 i=51 0.5 / 0.5+51/512 (b=0x4CC00000)
1126 : 0x6A63BD81, // 0.831168831168831 i=52 0.5 / 0.5+52/512 (b=0x4D000000)
1127 : 0x6A0B9944, // 0.828478964401295 i=53 0.5 / 0.5+53/512 (b=0x4D400000)
1128 : 0x69B4069B, // 0.825806451612903 i=54 0.5 / 0.5+54/512 (b=0x4D800000)
1129 : 0x695D041D, // 0.823151125401929 i=55 0.5 / 0.5+55/512 (b=0x4DC00000)
1130 : 0x69069069, // 0.820512820512820 i=56 0.5 / 0.5+56/512 (b=0x4E000000)
1131 : 0x68B0AA1F, // 0.817891373801917 i=57 0.5 / 0.5+57/512 (b=0x4E400000)
1132 : 0x685B4FE5, // 0.815286624203822 i=58 0.5 / 0.5+58/512 (b=0x4E800000)
1133 : 0x68068068, // 0.812698412698413 i=59 0.5 / 0.5+59/512 (b=0x4EC00000)
1134 : 0x67B23A54, // 0.810126582278481 i=60 0.5 / 0.5+60/512 (b=0x4F000000)
1135 : 0x675E7C5D, // 0.807570977917981 i=61 0.5 / 0.5+61/512 (b=0x4F400000)
1136 : 0x670B453B, // 0.805031446540881 i=62 0.5 / 0.5+62/512 (b=0x4F800000)
1137 : 0x66B893A9, // 0.802507836990596 i=63 0.5 / 0.5+63/512 (b=0x4FC00000)
1138 : 0x66666666, // 0.800000000000000 i=64 0.5 / 0.5+64/512 (b=0x50000000)
1139 : 0x6614BC36, // 0.797507788161994 i=65 0.5 / 0.5+65/512 (b=0x50400000)
1140 : 0x65C393E0, // 0.795031055900621 i=66 0.5 / 0.5+66/512 (b=0x50800000)
1141 : 0x6572EC2F, // 0.792569659442725 i=67 0.5 / 0.5+67/512 (b=0x50C00000)
1142 : 0x6522C3F3, // 0.790123456790123 i=68 0.5 / 0.5+68/512 (b=0x51000000)
1143 : 0x64D319FE, // 0.787692307692308 i=69 0.5 / 0.5+69/512 (b=0x51400000)
1144 : 0x6483ED27, // 0.785276073619632 i=70 0.5 / 0.5+70/512 (b=0x51800000)
1145 : 0x64353C48, // 0.782874617737003 i=71 0.5 / 0.5+71/512 (b=0x51C00000)
1146 : 0x63E7063E, // 0.780487804878049 i=72 0.5 / 0.5+72/512 (b=0x52000000)
1147 : 0x639949EB, // 0.778115501519757 i=73 0.5 / 0.5+73/512 (b=0x52400000)
1148 : 0x634C0634, // 0.775757575757576 i=74 0.5 / 0.5+74/512 (b=0x52800000)
1149 : 0x62FF3A01, // 0.773413897280967 i=75 0.5 / 0.5+75/512 (b=0x52C00000)
1150 : 0x62B2E43D, // 0.771084337349398 i=76 0.5 / 0.5+76/512 (b=0x53000000)
1151 : 0x626703D8, // 0.768768768768769 i=77 0.5 / 0.5+77/512 (b=0x53400000)
1152 : 0x621B97C2, // 0.766467065868264 i=78 0.5 / 0.5+78/512 (b=0x53800000)
1153 : 0x61D09EF3, // 0.764179104477612 i=79 0.5 / 0.5+79/512 (b=0x53C00000)
1154 : 0x61861861, // 0.761904761904762 i=80 0.5 / 0.5+80/512 (b=0x54000000)
1155 : 0x613C0309, // 0.759643916913947 i=81 0.5 / 0.5+81/512 (b=0x54400000)
1156 : 0x60F25DEA, // 0.757396449704142 i=82 0.5 / 0.5+82/512 (b=0x54800000)
1157 : 0x60A92806, // 0.755162241887906 i=83 0.5 / 0.5+83/512 (b=0x54C00000)
1158 : 0x60606060, // 0.752941176470588 i=84 0.5 / 0.5+84/512 (b=0x55000000)
1159 : 0x60180601, // 0.750733137829912 i=85 0.5 / 0.5+85/512 (b=0x55400000)
1160 : 0x5FD017F4, // 0.748538011695906 i=86 0.5 / 0.5+86/512 (b=0x55800000)
1161 : 0x5F889545, // 0.746355685131195 i=87 0.5 / 0.5+87/512 (b=0x55C00000)
1162 : 0x5F417D05, // 0.744186046511628 i=88 0.5 / 0.5+88/512 (b=0x56000000)
1163 : 0x5EFACE48, // 0.742028985507246 i=89 0.5 / 0.5+89/512 (b=0x56400000)
1164 : 0x5EB48823, // 0.739884393063584 i=90 0.5 / 0.5+90/512 (b=0x56800000)
1165 : 0x5E6EA9AE, // 0.737752161383285 i=91 0.5 / 0.5+91/512 (b=0x56C00000)
1166 : 0x5E293205, // 0.735632183908046 i=92 0.5 / 0.5+92/512 (b=0x57000000)
1167 : 0x5DE42046, // 0.733524355300860 i=93 0.5 / 0.5+93/512 (b=0x57400000)
1168 : 0x5D9F7390, // 0.731428571428571 i=94 0.5 / 0.5+94/512 (b=0x57800000)
1169 : 0x5D5B2B08, // 0.729344729344729 i=95 0.5 / 0.5+95/512 (b=0x57C00000)
1170 : 0x5D1745D1, // 0.727272727272727 i=96 0.5 / 0.5+96/512 (b=0x58000000)
1171 : 0x5CD3C315, // 0.725212464589235 i=97 0.5 / 0.5+97/512 (b=0x58400000)
1172 : 0x5C90A1FD, // 0.723163841807910 i=98 0.5 / 0.5+98/512 (b=0x58800000)
1173 : 0x5C4DE1B6, // 0.721126760563380 i=99 0.5 / 0.5+99/512 (b=0x58C00000)
1174 : 0x5C0B8170, // 0.719101123595506 i=100 0.5 / 0.5+100/512 (b=0x59000000)
1175 : 0x5BC9805B, // 0.717086834733894 i=101 0.5 / 0.5+101/512 (b=0x59400000)
1176 : 0x5B87DDAD, // 0.715083798882682 i=102 0.5 / 0.5+102/512 (b=0x59800000)
1177 : 0x5B46989A, // 0.713091922005571 i=103 0.5 / 0.5+103/512 (b=0x59C00000)
1178 : 0x5B05B05B, // 0.711111111111111 i=104 0.5 / 0.5+104/512 (b=0x5A000000)
1179 : 0x5AC5242A, // 0.709141274238227 i=105 0.5 / 0.5+105/512 (b=0x5A400000)
1180 : 0x5A84F345, // 0.707182320441989 i=106 0.5 / 0.5+106/512 (b=0x5A800000)
1181 : 0x5A451CEA, // 0.705234159779614 i=107 0.5 / 0.5+107/512 (b=0x5AC00000)
1182 : 0x5A05A05A, // 0.703296703296703 i=108 0.5 / 0.5+108/512 (b=0x5B000000)
1183 : 0x59C67CD8, // 0.701369863013699 i=109 0.5 / 0.5+109/512 (b=0x5B400000)
1184 : 0x5987B1A9, // 0.699453551912568 i=110 0.5 / 0.5+110/512 (b=0x5B800000)
1185 : 0x59493E14, // 0.697547683923706 i=111 0.5 / 0.5+111/512 (b=0x5BC00000)
1186 : 0x590B2164, // 0.695652173913043 i=112 0.5 / 0.5+112/512 (b=0x5C000000)
1187 : 0x58CD5AE2, // 0.693766937669377 i=113 0.5 / 0.5+113/512 (b=0x5C400000)
1188 : 0x588FE9DC, // 0.691891891891892 i=114 0.5 / 0.5+114/512 (b=0x5C800000)
1189 : 0x5852CDA0, // 0.690026954177898 i=115 0.5 / 0.5+115/512 (b=0x5CC00000)
1190 : 0x58160581, // 0.688172043010753 i=116 0.5 / 0.5+116/512 (b=0x5D000000)
1191 : 0x57D990D0, // 0.686327077747989 i=117 0.5 / 0.5+117/512 (b=0x5D400000)
1192 : 0x579D6EE3, // 0.684491978609626 i=118 0.5 / 0.5+118/512 (b=0x5D800000)
1193 : 0x57619F0F, // 0.682666666666667 i=119 0.5 / 0.5+119/512 (b=0x5DC00000)
1194 : 0x572620AE, // 0.680851063829787 i=120 0.5 / 0.5+120/512 (b=0x5E000000)
1195 : 0x56EAF319, // 0.679045092838196 i=121 0.5 / 0.5+121/512 (b=0x5E400000)
1196 : 0x56B015AC, // 0.677248677248677 i=122 0.5 / 0.5+122/512 (b=0x5E800000)
1197 : 0x567587C4, // 0.675461741424802 i=123 0.5 / 0.5+123/512 (b=0x5EC00000)
1198 : 0x563B48C2, // 0.673684210526316 i=124 0.5 / 0.5+124/512 (b=0x5F000000)
1199 : 0x56015805, // 0.671916010498688 i=125 0.5 / 0.5+125/512 (b=0x5F400000)
1200 : 0x55C7B4F1, // 0.670157068062827 i=126 0.5 / 0.5+126/512 (b=0x5F800000)
1201 : 0x558E5EE9, // 0.668407310704961 i=127 0.5 / 0.5+127/512 (b=0x5FC00000)
1202 : 0x55555555, // 0.666666666666667 i=128 0.5 / 0.5+128/512 (b=0x60000000)
1203 : 0x551C979A, // 0.664935064935065 i=129 0.5 / 0.5+129/512 (b=0x60400000)
1204 : 0x54E42523, // 0.663212435233161 i=130 0.5 / 0.5+130/512 (b=0x60800000)
1205 : 0x54ABFD5A, // 0.661498708010336 i=131 0.5 / 0.5+131/512 (b=0x60C00000)
1206 : 0x54741FAB, // 0.659793814432990 i=132 0.5 / 0.5+132/512 (b=0x61000000)
1207 : 0x543C8B84, // 0.658097686375321 i=133 0.5 / 0.5+133/512 (b=0x61400000)
1208 : 0x54054054, // 0.656410256410256 i=134 0.5 / 0.5+134/512 (b=0x61800000)
1209 : 0x53CE3D8B, // 0.654731457800512 i=135 0.5 / 0.5+135/512 (b=0x61C00000)
1210 : 0x5397829C, // 0.653061224489796 i=136 0.5 / 0.5+136/512 (b=0x62000000)
1211 : 0x53610EFB, // 0.651399491094148 i=137 0.5 / 0.5+137/512 (b=0x62400000)
1212 : 0x532AE21C, // 0.649746192893401 i=138 0.5 / 0.5+138/512 (b=0x62800000)
1213 : 0x52F4FB76, // 0.648101265822785 i=139 0.5 / 0.5+139/512 (b=0x62C00000)
1214 : 0x52BF5A81, // 0.646464646464647 i=140 0.5 / 0.5+140/512 (b=0x63000000)
1215 : 0x5289FEB5, // 0.644836272040302 i=141 0.5 / 0.5+141/512 (b=0x63400000)
1216 : 0x5254E78E, // 0.643216080402010 i=142 0.5 / 0.5+142/512 (b=0x63800000)
1217 : 0x52201488, // 0.641604010025063 i=143 0.5 / 0.5+143/512 (b=0x63C00000)
1218 : 0x51EB851E, // 0.640000000000000 i=144 0.5 / 0.5+144/512 (b=0x64000000)
1219 : 0x51B738D1, // 0.638403990024938 i=145 0.5 / 0.5+145/512 (b=0x64400000)
1220 : 0x51832F1F, // 0.636815920398010 i=146 0.5 / 0.5+146/512 (b=0x64800000)
1221 : 0x514F678B, // 0.635235732009926 i=147 0.5 / 0.5+147/512 (b=0x64C00000)
1222 : 0x511BE195, // 0.633663366336634 i=148 0.5 / 0.5+148/512 (b=0x65000000)
1223 : 0x50E89CC2, // 0.632098765432099 i=149 0.5 / 0.5+149/512 (b=0x65400000)
1224 : 0x50B59897, // 0.630541871921182 i=150 0.5 / 0.5+150/512 (b=0x65800000)
1225 : 0x5082D499, // 0.628992628992629 i=151 0.5 / 0.5+151/512 (b=0x65C00000)
1226 : 0x50505050, // 0.627450980392157 i=152 0.5 / 0.5+152/512 (b=0x66000000)
1227 : 0x501E0B44, // 0.625916870415648 i=153 0.5 / 0.5+153/512 (b=0x66400000)
1228 : 0x4FEC04FE, // 0.624390243902439 i=154 0.5 / 0.5+154/512 (b=0x66800000)
1229 : 0x4FBA3D0A, // 0.622871046228710 i=155 0.5 / 0.5+155/512 (b=0x66C00000)
1230 : 0x4F88B2F3, // 0.621359223300971 i=156 0.5 / 0.5+156/512 (b=0x67000000)
1231 : 0x4F576646, // 0.619854721549637 i=157 0.5 / 0.5+157/512 (b=0x67400000)
1232 : 0x4F265691, // 0.618357487922705 i=158 0.5 / 0.5+158/512 (b=0x67800000)
1233 : 0x4EF58364, // 0.616867469879518 i=159 0.5 / 0.5+159/512 (b=0x67C00000)
1234 : 0x4EC4EC4E, // 0.615384615384615 i=160 0.5 / 0.5+160/512 (b=0x68000000)
1235 : 0x4E9490E1, // 0.613908872901679 i=161 0.5 / 0.5+161/512 (b=0x68400000)
1236 : 0x4E6470B0, // 0.612440191387560 i=162 0.5 / 0.5+162/512 (b=0x68800000)
1237 : 0x4E348B4D, // 0.610978520286396 i=163 0.5 / 0.5+163/512 (b=0x68C00000)
1238 : 0x4E04E04E, // 0.609523809523810 i=164 0.5 / 0.5+164/512 (b=0x69000000)
1239 : 0x4DD56F47, // 0.608076009501188 i=165 0.5 / 0.5+165/512 (b=0x69400000)
1240 : 0x4DA637CF, // 0.606635071090047 i=166 0.5 / 0.5+166/512 (b=0x69800000)
1241 : 0x4D77397E, // 0.605200945626478 i=167 0.5 / 0.5+167/512 (b=0x69C00000)
1242 : 0x4D4873EC, // 0.603773584905660 i=168 0.5 / 0.5+168/512 (b=0x6A000000)
1243 : 0x4D19E6B3, // 0.602352941176471 i=169 0.5 / 0.5+169/512 (b=0x6A400000)
1244 : 0x4CEB916D, // 0.600938967136150 i=170 0.5 / 0.5+170/512 (b=0x6A800000)
1245 : 0x4CBD73B5, // 0.599531615925059 i=171 0.5 / 0.5+171/512 (b=0x6AC00000)
1246 : 0x4C8F8D28, // 0.598130841121495 i=172 0.5 / 0.5+172/512 (b=0x6B000000)
1247 : 0x4C61DD63, // 0.596736596736597 i=173 0.5 / 0.5+173/512 (b=0x6B400000)
1248 : 0x4C346404, // 0.595348837209302 i=174 0.5 / 0.5+174/512 (b=0x6B800000)
1249 : 0x4C0720AB, // 0.593967517401392 i=175 0.5 / 0.5+175/512 (b=0x6BC00000)
1250 : 0x4BDA12F6, // 0.592592592592593 i=176 0.5 / 0.5+176/512 (b=0x6C000000)
1251 : 0x4BAD3A87, // 0.591224018475751 i=177 0.5 / 0.5+177/512 (b=0x6C400000)
1252 : 0x4B809701, // 0.589861751152074 i=178 0.5 / 0.5+178/512 (b=0x6C800000)
1253 : 0x4B542804, // 0.588505747126437 i=179 0.5 / 0.5+179/512 (b=0x6CC00000)
1254 : 0x4B27ED36, // 0.587155963302752 i=180 0.5 / 0.5+180/512 (b=0x6D000000)
1255 : 0x4AFBE639, // 0.585812356979405 i=181 0.5 / 0.5+181/512 (b=0x6D400000)
1256 : 0x4AD012B4, // 0.584474885844749 i=182 0.5 / 0.5+182/512 (b=0x6D800000)
1257 : 0x4AA4724B, // 0.583143507972665 i=183 0.5 / 0.5+183/512 (b=0x6DC00000)
1258 : 0x4A7904A7, // 0.581818181818182 i=184 0.5 / 0.5+184/512 (b=0x6E000000)
1259 : 0x4A4DC96E, // 0.580498866213152 i=185 0.5 / 0.5+185/512 (b=0x6E400000)
1260 : 0x4A22C04A, // 0.579185520361991 i=186 0.5 / 0.5+186/512 (b=0x6E800000)
1261 : 0x49F7E8E2, // 0.577878103837472 i=187 0.5 / 0.5+187/512 (b=0x6EC00000)
1262 : 0x49CD42E2, // 0.576576576576577 i=188 0.5 / 0.5+188/512 (b=0x6F000000)
1263 : 0x49A2CDF3, // 0.575280898876405 i=189 0.5 / 0.5+189/512 (b=0x6F400000)
1264 : 0x497889C2, // 0.573991031390135 i=190 0.5 / 0.5+190/512 (b=0x6F800000)
1265 : 0x494E75FA, // 0.572706935123042 i=191 0.5 / 0.5+191/512 (b=0x6FC00000)
1266 : 0x49249249, // 0.571428571428571 i=192 0.5 / 0.5+192/512 (b=0x70000000)
1267 : 0x48FADE5C, // 0.570155902004454 i=193 0.5 / 0.5+193/512 (b=0x70400000)
1268 : 0x48D159E2, // 0.568888888888889 i=194 0.5 / 0.5+194/512 (b=0x70800000)
1269 : 0x48A8048A, // 0.567627494456763 i=195 0.5 / 0.5+195/512 (b=0x70C00000)
1270 : 0x487EDE04, // 0.566371681415929 i=196 0.5 / 0.5+196/512 (b=0x71000000)
1271 : 0x4855E601, // 0.565121412803532 i=197 0.5 / 0.5+197/512 (b=0x71400000)
1272 : 0x482D1C31, // 0.563876651982379 i=198 0.5 / 0.5+198/512 (b=0x71800000)
1273 : 0x48048048, // 0.562637362637363 i=199 0.5 / 0.5+199/512 (b=0x71C00000)
1274 : 0x47DC11F7, // 0.561403508771930 i=200 0.5 / 0.5+200/512 (b=0x72000000)
1275 : 0x47B3D0F1, // 0.560175054704595 i=201 0.5 / 0.5+201/512 (b=0x72400000)
1276 : 0x478BBCEC, // 0.558951965065502 i=202 0.5 / 0.5+202/512 (b=0x72800000)
1277 : 0x4763D59C, // 0.557734204793028 i=203 0.5 / 0.5+203/512 (b=0x72C00000)
1278 : 0x473C1AB6, // 0.556521739130435 i=204 0.5 / 0.5+204/512 (b=0x73000000)
1279 : 0x47148BF0, // 0.555314533622560 i=205 0.5 / 0.5+205/512 (b=0x73400000)
1280 : 0x46ED2901, // 0.554112554112554 i=206 0.5 / 0.5+206/512 (b=0x73800000)
1281 : 0x46C5F19F, // 0.552915766738661 i=207 0.5 / 0.5+207/512 (b=0x73C00000)
1282 : 0x469EE584, // 0.551724137931034 i=208 0.5 / 0.5+208/512 (b=0x74000000)
1283 : 0x46780467, // 0.550537634408602 i=209 0.5 / 0.5+209/512 (b=0x74400000)
1284 : 0x46514E02, // 0.549356223175966 i=210 0.5 / 0.5+210/512 (b=0x74800000)
1285 : 0x462AC20E, // 0.548179871520343 i=211 0.5 / 0.5+211/512 (b=0x74C00000)
1286 : 0x46046046, // 0.547008547008547 i=212 0.5 / 0.5+212/512 (b=0x75000000)
1287 : 0x45DE2864, // 0.545842217484009 i=213 0.5 / 0.5+213/512 (b=0x75400000)
1288 : 0x45B81A25, // 0.544680851063830 i=214 0.5 / 0.5+214/512 (b=0x75800000)
1289 : 0x45923543, // 0.543524416135881 i=215 0.5 / 0.5+215/512 (b=0x75C00000)
1290 : 0x456C797D, // 0.542372881355932 i=216 0.5 / 0.5+216/512 (b=0x76000000)
1291 : 0x4546E68F, // 0.541226215644820 i=217 0.5 / 0.5+217/512 (b=0x76400000)
1292 : 0x45217C38, // 0.540084388185654 i=218 0.5 / 0.5+218/512 (b=0x76800000)
1293 : 0x44FC3A34, // 0.538947368421053 i=219 0.5 / 0.5+219/512 (b=0x76C00000)
1294 : 0x44D72044, // 0.537815126050420 i=220 0.5 / 0.5+220/512 (b=0x77000000)
1295 : 0x44B22E27, // 0.536687631027254 i=221 0.5 / 0.5+221/512 (b=0x77400000)
1296 : 0x448D639D, // 0.535564853556485 i=222 0.5 / 0.5+222/512 (b=0x77800000)
1297 : 0x4468C066, // 0.534446764091858 i=223 0.5 / 0.5+223/512 (b=0x77C00000)
1298 : 0x44444444, // 0.533333333333333 i=224 0.5 / 0.5+224/512 (b=0x78000000)
1299 : 0x441FEEF8, // 0.532224532224532 i=225 0.5 / 0.5+225/512 (b=0x78400000)
1300 : 0x43FBC043, // 0.531120331950207 i=226 0.5 / 0.5+226/512 (b=0x78800000)
1301 : 0x43D7B7EA, // 0.530020703933747 i=227 0.5 / 0.5+227/512 (b=0x78C00000)
1302 : 0x43B3D5AF, // 0.528925619834711 i=228 0.5 / 0.5+228/512 (b=0x79000000)
1303 : 0x43901956, // 0.527835051546392 i=229 0.5 / 0.5+229/512 (b=0x79400000)
1304 : 0x436C82A2, // 0.526748971193416 i=230 0.5 / 0.5+230/512 (b=0x79800000)
1305 : 0x43491158, // 0.525667351129363 i=231 0.5 / 0.5+231/512 (b=0x79C00000)
1306 : 0x4325C53E, // 0.524590163934426 i=232 0.5 / 0.5+232/512 (b=0x7A000000)
1307 : 0x43029E1A, // 0.523517382413088 i=233 0.5 / 0.5+233/512 (b=0x7A400000)
1308 : 0x42DF9BB0, // 0.522448979591837 i=234 0.5 / 0.5+234/512 (b=0x7A800000)
1309 : 0x42BCBDC8, // 0.521384928716904 i=235 0.5 / 0.5+235/512 (b=0x7AC00000)
1310 : 0x429A0429, // 0.520325203252033 i=236 0.5 / 0.5+236/512 (b=0x7B000000)
1311 : 0x42776E9A, // 0.519269776876268 i=237 0.5 / 0.5+237/512 (b=0x7B400000)
1312 : 0x4254FCE4, // 0.518218623481781 i=238 0.5 / 0.5+238/512 (b=0x7B800000)
1313 : 0x4232AECD, // 0.517171717171717 i=239 0.5 / 0.5+239/512 (b=0x7BC00000)
1314 : 0x42108421, // 0.516129032258065 i=240 0.5 / 0.5+240/512 (b=0x7C000000)
1315 : 0x41EE7CA6, // 0.515090543259557 i=241 0.5 / 0.5+241/512 (b=0x7C400000)
1316 : 0x41CC9829, // 0.514056224899598 i=242 0.5 / 0.5+242/512 (b=0x7C800000)
1317 : 0x41AAD671, // 0.513026052104208 i=243 0.5 / 0.5+243/512 (b=0x7CC00000)
1318 : 0x4189374B, // 0.512000000000000 i=244 0.5 / 0.5+244/512 (b=0x7D000000)
1319 : 0x4167BA81, // 0.510978043912176 i=245 0.5 / 0.5+245/512 (b=0x7D400000)
1320 : 0x41465FDF, // 0.509960159362550 i=246 0.5 / 0.5+246/512 (b=0x7D800000)
1321 : 0x41252730, // 0.508946322067594 i=247 0.5 / 0.5+247/512 (b=0x7DC00000)
1322 : 0x41041041, // 0.507936507936508 i=248 0.5 / 0.5+248/512 (b=0x7E000000)
1323 : 0x40E31ADE, // 0.506930693069307 i=249 0.5 / 0.5+249/512 (b=0x7E400000)
1324 : 0x40C246D4, // 0.505928853754941 i=250 0.5 / 0.5+250/512 (b=0x7E800000)
1325 : 0x40A193F1, // 0.504930966469428 i=251 0.5 / 0.5+251/512 (b=0x7EC00000)
1326 : 0x40810204, // 0.503937007874016 i=252 0.5 / 0.5+252/512 (b=0x7F000000)
1327 : 0x406090D9, // 0.502946954813359 i=253 0.5 / 0.5+253/512 (b=0x7F400000)
1328 : 0x40404040, // 0.501960784313725 i=254 0.5 / 0.5+254/512 (b=0x7F800000)
1329 : 0x40201008 // 0.500978473581213 i=255 0.5 / 0.5+255/512 (b=0x7FC00000)
1330 : };
1331 :
1332 :
1333 : /*
1334 : * Fractional multiplication of signed a and b, both in Q31. The result is doubled.
1335 : * Note: in this test, saturation is not needed.
1336 : * BASOP weights: 3
1337 : */
1338 :
1339 1025744049 : static Word32 L_dmult( Word32 L_var1, Word32 L_var2 )
1340 : {
1341 1025744049 : Word64 L64_var1 = W_mult0_32_32( L_var1, L_var2 );
1342 1025744049 : L64_var1 = W_shr( L64_var1, 30 );
1343 1025744049 : return W_extract_l( L64_var1 );
1344 : }
1345 :
1346 : /*
1347 : * 32 by 32 bit division, following the Newton / Raphson method.
1348 : * Usage of this low-level procedure by the caller:
1349 : * 1. Numerator <num> can use the full range of signed 32-bit datatypes: 0x80000000...0x7FFFFFFF
1350 : * Note: Since <num> is not normalized here, but it is multplied to the reciprocal of the
1351 : * denominator, the caller should use a normalized <num> and handle any exponent outside.
1352 : * 2. Denominator <den> must be normalized into range 0x40000001 to 0x7FFFFFFF (all positive)
1353 : * Note: In case of den=0x40000000, the caller is not allowed to call the division routine,
1354 : since the result is known.
1355 : * Note: num / 0x40000000 equals to num with exp += 1.
1356 : * 3. The result is in range 0x40000000 to 0x7FFFFFFF, finally multiplied by <num>.
1357 : * BASOP weights: 24 (incl. L_dmult)
1358 : */
1359 :
1360 341914683 : Word32 div_w_newton( Word32 num, Word32 den )
1361 : {
1362 : Word32 x0, x1, x2, x3, diff, result;
1363 :
1364 341914683 : x0 = div_w_newton_lookup[sub( extract_l( L_shr( den, 22 ) ), 256 )];
1365 341914683 : move32();
1366 :
1367 341914683 : diff = L_sub( 0x40000000, Mpy_32_32( den, x0 ) );
1368 :
1369 341914683 : x1 = L_add( x0, L_dmult( x0, diff ) );
1370 341914683 : diff = L_sub( 0x40000000, Mpy_32_32( den, x1 ) );
1371 :
1372 341914683 : x2 = L_add( x1, L_dmult( x1, diff ) );
1373 341914683 : diff = L_sub( 0x40000000, Mpy_32_32( den, x2 ) );
1374 :
1375 341914683 : x3 = L_add( x2, L_dmult( x2, diff ) );
1376 :
1377 341914683 : result = Mpy_32_32( num, x3 );
1378 :
1379 341914683 : return result;
1380 : }
1381 :
1382 : /*
1383 : * 32 / 32 division
1384 : * Usage of this global procedure by the caller:
1385 : * 1. Numerator <num> can use the full range of signed 32-bit datatypes: 0x80000000...0x7FFFFFFF
1386 : * Note: Since <num> is not normalized here, but it is multplied to the reciprocal of the
1387 : * denominator, the caller should use a normalized <num> and handle any exponent outside.
1388 : * 2. Denominator can use the full range of signed 32-bit datatypes: 0x80000000...0x7FFFFFFF
1389 : * except den=0x00000000. In case of 0x80000000, it becomes internally negated to 0x7FFFFFFF.
1390 : * 3. The result is 0x00000000 (*s=0)for <num> equals 0x00000000.
1391 : * The result is rather left aligned, with up to 1 bit headroom.
1392 : * 4. The result exponent is stored in s[0]
1393 : * BASOP weights: 41 (incl. div_w_newton)
1394 : */
1395 :
1396 359547276 : Word32 BASOP_Util_Divide3232_Scale_newton( Word32 x, Word32 y, Word16 *s )
1397 : {
1398 : Word32 z;
1399 : Word16 sx;
1400 : Word16 sy;
1401 : Word32 sign;
1402 :
1403 359547276 : assert( y != (Word32) 0 );
1404 :
1405 : /* Early exit, if numerator is zero */
1406 359547276 : IF( x == (Word32) 0 )
1407 : {
1408 2346924 : *s = 0;
1409 2346924 : return ( (Word32) 0 );
1410 : }
1411 : #ifdef FIX_USAN_BASOP_UTIL_DIVIDE3232
1412 357200352 : IF( EQ_32( y, (Word32) 0x80000000 ) )
1413 : #else
1414 : IF( EQ_32( y, 0x80000000 ) )
1415 : #endif
1416 : {
1417 : /* Division by -1.0: same as negation of numerator */
1418 : /* Return normalized negated numerator */
1419 24270 : sx = norm_l( x );
1420 24270 : x = L_shl( x, sx );
1421 24270 : *s = negate( sx );
1422 24270 : return L_negate( x );
1423 : }
1424 357176082 : sign = y;
1425 357176082 : move32();
1426 357176082 : if ( y < 0 )
1427 : {
1428 11320819 : y = L_negate( y );
1429 : }
1430 :
1431 : /* Normalize numerator */
1432 357176082 : sx = norm_l( x );
1433 357176082 : x = L_shl( x, sx );
1434 :
1435 : /* Normalize denominator */
1436 357176082 : sy = norm_l( y );
1437 357176082 : y = L_shl( y, sy );
1438 :
1439 : /* Store exponent: + 1 for div_w_newton computing 0.5*num/den */
1440 357176082 : *s = sub( add( sy, 1 ), sx );
1441 357176082 : move16();
1442 :
1443 : /* Special treatment for den=0x40000000 */
1444 : /* Result is known: z=2*num */
1445 357176082 : IF( EQ_32( y, 0x40000000 ) )
1446 : {
1447 15261399 : if ( sign < 0 )
1448 : {
1449 3 : x = L_negate( x );
1450 : }
1451 15261399 : return x;
1452 : }
1453 :
1454 : /* Invoke division applying Newton/Raphson-Algorithm */
1455 341914683 : z = div_w_newton( x, y );
1456 :
1457 341914683 : if ( sign < 0 )
1458 : {
1459 11320816 : z = L_negate( z );
1460 : }
1461 :
1462 341914683 : return z;
1463 : }
1464 :
1465 254949633 : Word16 BASOP_Util_Divide3232_Scale( Word32 x, Word32 y, Word16 *s )
1466 : {
1467 : Word16 z;
1468 : Word16 sy;
1469 :
1470 :
1471 254949633 : sy = norm_l( y );
1472 254949633 : if ( sy > 0 )
1473 : {
1474 227872625 : sy = sub( sy, 1 );
1475 : }
1476 254949633 : y = L_shl( y, sy );
1477 :
1478 254949633 : z = BASOP_Util_Divide3216_Scale( x, extract_h( y ), s );
1479 254949633 : move16();
1480 254949633 : *s = add( *s, sy );
1481 :
1482 :
1483 254949633 : return ( z );
1484 : }
1485 :
1486 :
1487 14741673 : Word16 BASOP_Util_Divide1616_Scale( Word16 x, Word16 y, Word16 *s )
1488 : {
1489 : Word16 z;
1490 : Word16 sx;
1491 : Word16 sy;
1492 : Word16 sign;
1493 :
1494 :
1495 : /* assert (x >= (Word16)0); */
1496 14741673 : assert( y != (Word16) 0 );
1497 :
1498 14741673 : sign = 0;
1499 14741673 : move16();
1500 :
1501 14741673 : IF( x < 0 )
1502 : {
1503 653531 : x = negate( x );
1504 653531 : sign = s_xor( sign, 1 );
1505 : }
1506 :
1507 14741673 : IF( y < 0 )
1508 : {
1509 0 : y = negate( y );
1510 0 : sign = s_xor( sign, 1 );
1511 : }
1512 :
1513 14741673 : IF( x == (Word16) 0 )
1514 : {
1515 758256 : move16();
1516 758256 : *s = 0;
1517 :
1518 :
1519 758256 : return ( (Word16) 0 );
1520 : }
1521 :
1522 13983417 : sx = norm_s( x );
1523 13983417 : x = shl( x, sx );
1524 13983417 : x = shr( x, 1 );
1525 13983417 : move16();
1526 13983417 : *s = sub( 1, sx );
1527 :
1528 13983417 : sy = norm_s( y );
1529 13983417 : y = shl( y, sy );
1530 13983417 : move16();
1531 13983417 : *s = add( *s, sy );
1532 :
1533 13983417 : z = div_s( x, y );
1534 :
1535 13983417 : if ( sign != 0 )
1536 : {
1537 653531 : z = negate( z );
1538 : }
1539 :
1540 :
1541 13983417 : return z;
1542 : }
1543 :
1544 :
1545 0 : void set_val_Word16(
1546 : Word16 X[],
1547 : const Word16 val,
1548 : Word16 n )
1549 : {
1550 : Word16 i;
1551 :
1552 :
1553 0 : FOR( i = 0; i < n; i++ )
1554 : {
1555 0 : X[i] = val;
1556 0 : move16();
1557 : }
1558 :
1559 :
1560 0 : return;
1561 : }
1562 :
1563 0 : void set_val_Word32(
1564 : Word32 X[],
1565 : const Word32 val,
1566 : Word16 n )
1567 : {
1568 : Word16 i;
1569 :
1570 :
1571 0 : FOR( i = 0; i < n; i++ )
1572 : {
1573 0 : X[i] = val;
1574 0 : move32();
1575 : }
1576 :
1577 :
1578 0 : return;
1579 : }
1580 :
1581 580567 : Word16 mult0(
1582 : Word16 x,
1583 : Word16 y )
1584 : {
1585 580567 : return extract_l( L_mult0( x, y ) );
1586 : }
1587 :
1588 0 : void copyWord8( const Word8 *src, Word8 *dst, const Word32 n )
1589 : {
1590 : Word32 i;
1591 :
1592 :
1593 0 : FOR( i = 0; i < n; i++ )
1594 : {
1595 0 : dst[i] = src[i];
1596 0 : move16();
1597 : }
1598 0 : }
1599 :
1600 :
1601 0 : void set_zero_Word8( Word8 X[], Word32 n )
1602 : {
1603 : Word32 i;
1604 :
1605 :
1606 0 : FOR( i = 0; i < n; i++ )
1607 : {
1608 0 : X[i] = 0;
1609 0 : move16();
1610 : }
1611 0 : }
1612 :
1613 :
1614 0 : Word32 L_mult0_3216( Word32 x, Word16 y )
1615 : {
1616 : UWord16 mpy_low16;
1617 : Word32 mpy_high32;
1618 :
1619 :
1620 0 : Mpy_32_16_ss( x, y, &mpy_high32, &mpy_low16 );
1621 :
1622 0 : mpy_high32 = L_add( L_shl( mpy_high32, 15 ), L_lshr( L_deposit_h( mpy_low16 ), 17 ) );
1623 :
1624 :
1625 0 : return mpy_high32;
1626 : }
1627 :
1628 0 : Word16 BASOP_util_norm_l_dim2_cplx( const Word32 *const *re, /*!< Real part of 32 Bit input */
1629 : const Word32 *const *im, /*!< Imag part if 32 Bit input */
1630 : Word16 startBand, /*!< start band of cplx data */
1631 : Word16 stopBand, /*!< stop band of cplx data */
1632 : Word16 startSlot, /*!< start slot of cplx data */
1633 : Word16 stopSlot /*!< stop slot of cplx data */
1634 : )
1635 : {
1636 : Word16 col;
1637 : Word16 band;
1638 : Word16 maxShift;
1639 : Word32 maxVal;
1640 :
1641 :
1642 0 : maxVal = L_deposit_l( 1 );
1643 :
1644 0 : FOR( col = startSlot; col < stopSlot; col++ )
1645 : {
1646 0 : FOR( band = startBand; band < stopBand; band++ )
1647 : {
1648 0 : maxVal = L_max( maxVal, L_abs( re[col][band] ) );
1649 0 : maxVal = L_max( maxVal, L_abs( im[col][band] ) );
1650 : }
1651 : }
1652 0 : maxShift = norm_l( maxVal );
1653 :
1654 :
1655 0 : return ( maxShift );
1656 : }
1657 :
1658 0 : Word16 BASOP_util_norm_s_bands2shift( Word16 x )
1659 : {
1660 : Word16 shift;
1661 :
1662 0 : shift = sub( WORD16_BITS - 1, norm_s( negate( x ) ) );
1663 :
1664 0 : return ( shift );
1665 : }
1666 :
1667 : #define SINETAB SineTable512_fx
1668 : #define LD 9
1669 :
1670 : /*
1671 : * Calculates coarse lookup values for sine/cosine and residual angle.
1672 : * \param x angle in radians with exponent = 2 or as radix 2 with exponent = 0.
1673 : * \param scale shall always be 2
1674 : * \param sine pointer to where the sine lookup value is stored into
1675 : * \param cosine pointer to where the cosine lookup value is stored into
1676 : * \param flag_radix2 flag indicating radix 2 angle if non-zero.
1677 : */
1678 67107959 : static Word16 fixp_sin_cos_residual_16(
1679 : Word16 x,
1680 : const Word16 scale,
1681 : Word16 *sine,
1682 : Word16 *cosine,
1683 : Word8 flag_radix2 )
1684 : {
1685 : Word16 residual;
1686 : Word16 s;
1687 : Word16 ssign;
1688 : Word16 csign;
1689 67107959 : Word16 tmp, cl = 0, sl = 0;
1690 67107959 : const Word16 shift = 15 - LD - 1 - scale;
1691 :
1692 67107959 : if ( flag_radix2 == 0 )
1693 : {
1694 3283825 : x = mult_r( x, FL2WORD16( 1.0 / EVS_PI ) );
1695 : }
1696 67107959 : s = shr( x, shift );
1697 :
1698 67107959 : residual = s_and( x, ( 1 << shift ) - 1 );
1699 : /* We assume "2+scale" is a constant */
1700 67107959 : residual = shl( residual, 2 + scale );
1701 67107959 : residual = mult_r( residual, FL2WORD16( EVS_PI / 4.0 ) );
1702 :
1703 : /* Sine sign symmetry */
1704 67107959 : ssign = s_and( s, ( 1 << LD ) << 1 );
1705 :
1706 : /* Cosine sign symmetry */
1707 67107959 : csign = s_and( add( s, ( 1 << LD ) ), ( 1 << LD ) << 1 );
1708 :
1709 : /* Modulo EVS_PI */
1710 67107959 : s = s_and( s, ( 2 << LD ) - 1 );
1711 :
1712 : /* EVS_PI/2 symmetry */
1713 67107959 : s = s_min( s, sub( 2 << LD, s ) );
1714 :
1715 : {
1716 67107959 : tmp = s_min( sub( 1 << LD, s ), s );
1717 67107959 : s = sub( tmp, s );
1718 :
1719 67107959 : if ( !s )
1720 : {
1721 42697424 : move16();
1722 42697424 : sl = SINETAB[tmp].v.im;
1723 : }
1724 67107959 : if ( !s )
1725 : {
1726 42697424 : move16();
1727 42697424 : cl = SINETAB[tmp].v.re;
1728 : }
1729 67107959 : if ( s )
1730 : {
1731 24410535 : move16();
1732 24410535 : sl = SINETAB[tmp].v.re;
1733 : }
1734 67107959 : if ( s )
1735 : {
1736 24410535 : move16();
1737 24410535 : cl = SINETAB[tmp].v.im;
1738 : }
1739 :
1740 67107959 : if ( ssign )
1741 : {
1742 20441178 : sl = negate( sl );
1743 : }
1744 67107959 : if ( csign )
1745 : {
1746 18405099 : cl = negate( cl );
1747 : }
1748 :
1749 67107959 : move16();
1750 67107959 : move16();
1751 67107959 : *sine = sl;
1752 67107959 : *cosine = cl;
1753 : }
1754 :
1755 67107959 : return residual;
1756 : }
1757 :
1758 3283825 : Word16 getCosWord16( Word16 theta )
1759 : {
1760 : Word16 result, residual, sine, cosine;
1761 :
1762 3283825 : residual = fixp_sin_cos_residual_16( theta, 2, &sine, &cosine, 0 );
1763 : /* This negation prevents the subsequent addition from overflow */
1764 : /* The negation cannot overflow, sine is in range [0x0..0x7FFF] */
1765 3283825 : sine = negate( sine );
1766 3283825 : result = mac_r( L_mult0( sine, residual ), cosine, 16384 );
1767 :
1768 :
1769 3283825 : return result;
1770 : }
1771 :
1772 609818 : Word16 getSinWord16( Word16 theta )
1773 : {
1774 : Word16 sine;
1775 609818 : Word32 theta_new = L_sub( EVS_PI_BY_2_FX, theta );
1776 : Word16 l_theta;
1777 609818 : IF( GT_32( theta_new, EVS_PI_FX ) )
1778 : {
1779 56878 : l_theta = extract_l( L_sub( L_sub( theta_new, EVS_PI_FX ), EVS_PI_FX ) );
1780 : }
1781 552940 : ELSE IF( LT_32( theta_new, -EVS_PI_FX ) )
1782 : {
1783 0 : l_theta = extract_l( L_add( L_add( theta_new, EVS_PI_FX ), EVS_PI_FX ) );
1784 : }
1785 : ELSE
1786 : {
1787 552940 : l_theta = extract_l( theta_new );
1788 : }
1789 609818 : sine = getCosWord16( l_theta );
1790 609818 : IF( EQ_16( sine, ONE_IN_Q14 ) )
1791 : {
1792 7032 : sine = MAX_16;
1793 : }
1794 : ELSE
1795 : {
1796 602786 : sine = shl( sine, 1 );
1797 : }
1798 609818 : return sine;
1799 : }
1800 :
1801 63824134 : Word16 getCosWord16R2(
1802 : Word16 theta )
1803 : {
1804 : Word16 result, residual, sine, cosine;
1805 :
1806 63824134 : residual = fixp_sin_cos_residual_16( theta, 1, &sine, &cosine, 1 );
1807 : /* This negation prevents the subsequent addition from overflow */
1808 : /* The negation cannot overflow, sine is in range [0x0..0x7FFF] */
1809 : BASOP_SATURATE_WARNING_OFF
1810 63824134 : sine = negate( sine );
1811 : /* Saturation has been included based on the recommendation from the PC group */
1812 63824134 : result = msu_r_sat( L_mult( sine, residual ), cosine, -32768 );
1813 : BASOP_SATURATE_WARNING_ON
1814 :
1815 63824134 : return result;
1816 : }
1817 :
1818 :
1819 26909491 : Word16 getSineWord16R2( Word16 theta )
1820 : {
1821 26909491 : IF( EQ_16( theta, (Word16) 0 ) )
1822 : {
1823 7256739 : return 0;
1824 : }
1825 19652752 : ELSE IF( LT_16( theta, (Word16) -24576 ) )
1826 : {
1827 109433 : theta = add( add( theta, (Word16) 32767 ), (Word16) 1 );
1828 : }
1829 19652752 : return getCosWord16R2( sub( 8192, theta ) );
1830 : }
1831 :
1832 : /*
1833 : * Calculate Integer Square Root of 'val'. This is the equivalent of (int)sqrt(val);
1834 : * The return value will be truncated to the lowest integer (throwing away the fractionnal part.
1835 : *
1836 : * There are many ways to do this. The approach here is to use a simple function to get a
1837 : * 1st estimate of (int)sqrt(val) and then correct this estimate if it is too low or too high.
1838 : *
1839 : * Using Word16, the range of 'val' is limited to roughly 2^30.
1840 : *
1841 : * Complexity: Worst=31Clks, Best=27Clks
1842 : */
1843 1694626 : Word16 getSqrtWord32( Word32 val )
1844 : {
1845 : Word32 L_temp, L_temp2;
1846 : Word16 temp, temp2;
1847 : Word16 exp, exp2;
1848 :
1849 : /* Calc Approximation */
1850 1694626 : exp2 = norm_l( val );
1851 1694626 : L_temp2 = L_shl( val, exp2 );
1852 1694626 : exp = sub( 31 - 32, exp2 );
1853 1694626 : L_temp = Isqrt_lc( L_temp2, &exp ); /* 12 clks */
1854 :
1855 1694626 : temp = round_fx_sat( L_temp );
1856 1694626 : L_temp = Mpy_32_16_1( L_temp2, temp ); /* 2 clks */
1857 :
1858 1694626 : L_temp = L_shl( L_temp, sub( exp, exp2 ) );
1859 :
1860 : /* The Approximation Error Range is -1..+7, so Too Low by 1 or Up to Too High by 7 */
1861 1694626 : temp = round_fx( L_temp );
1862 :
1863 : /* Too High? */
1864 1694626 : if ( L_msu0( val, temp, temp ) < 0 )
1865 : {
1866 : /* Reduce by 2 */
1867 1386865 : temp = sub( temp, 2 );
1868 : }
1869 : /* Too High? */
1870 1694626 : if ( L_msu0( val, temp, temp ) < 0 )
1871 : {
1872 : /* Reduce by 2 */
1873 156155 : temp = sub( temp, 2 );
1874 : }
1875 : /* Too High? */
1876 1694626 : if ( L_msu0( val, temp, temp ) < 0 )
1877 : {
1878 : /* Reduce by 2 */
1879 327 : temp = sub( temp, 2 );
1880 : }
1881 : /* Too High? */
1882 1694626 : if ( L_msu0( val, temp, temp ) < 0 )
1883 : {
1884 : /* Reduce by 1 */
1885 3 : temp = sub( temp, 1 );
1886 : }
1887 :
1888 : /* Try +1 */
1889 1694626 : temp2 = add( temp, 1 );
1890 : /* It fits? */
1891 1694626 : if ( L_msu0( val, temp2, temp2 ) >= 0 )
1892 : {
1893 : /* Yes */
1894 1167074 : temp = temp2;
1895 1167074 : move16();
1896 : }
1897 1694626 : return temp;
1898 : }
1899 0 : Word16 findIndexOfMinWord32( Word32 *x, const Word16 len )
1900 : {
1901 : Word16 i, indx;
1902 :
1903 :
1904 0 : indx = 0;
1905 0 : move16();
1906 0 : FOR( i = 1; i < len; i++ )
1907 : {
1908 0 : if ( LT_32( x[i], x[indx] ) )
1909 : {
1910 0 : indx = i;
1911 0 : move16();
1912 : }
1913 : }
1914 :
1915 :
1916 0 : return indx;
1917 : }
1918 :
1919 0 : Word16 findIndexOfMinWord64( Word64 *x, const Word16 len )
1920 : {
1921 : Word16 i, indx;
1922 :
1923 :
1924 0 : indx = 0;
1925 0 : move16();
1926 0 : FOR( i = 1; i < len; i++ )
1927 : {
1928 0 : if ( LT_64( x[i], x[indx] ) )
1929 : {
1930 0 : indx = i;
1931 0 : move16();
1932 : }
1933 : }
1934 :
1935 :
1936 0 : return indx;
1937 : }
1938 :
1939 :
1940 1795930796 : Word16 imult1616( Word16 x, Word16 y )
1941 : {
1942 1795930796 : assert( (int) x * (int) y < 32768 && (int) x * (int) y >= -32768 );
1943 1795930796 : return i_mult( x, y );
1944 : }
1945 :
1946 4012839359 : Word32 imult3216( Word32 x, Word16 y )
1947 : {
1948 : Word32 mh;
1949 : UWord16 ml;
1950 :
1951 4012839359 : Mpy_32_16_ss( x, y, &mh, &ml );
1952 :
1953 4012839359 : mh = L_shl( mh, 15 );
1954 4012839359 : ml = (UWord16) lshr( (Word16) ml, 1 );
1955 :
1956 4012839359 : return L_or( mh, L_deposit_l( ml ) );
1957 : }
1958 :
1959 :
1960 0 : Word16 idiv1616U_IVAS(
1961 : Word16 x,
1962 : Word16 y )
1963 : {
1964 : Word16 s;
1965 : Word16 tmp;
1966 :
1967 : /* make y > x */
1968 0 : s = add( sub( norm_s( y ), norm_s( x ) ), 1 );
1969 0 : s = s_max( s, 0 );
1970 :
1971 : BASOP_SATURATE_WARNING_OFF
1972 0 : y = shl( y, s );
1973 : BASOP_SATURATE_WARNING_ON
1974 :
1975 : /* divide and shift */
1976 0 : tmp = div_s( x, y );
1977 0 : y = shr( tmp, sub( 15, s ) );
1978 :
1979 0 : return y;
1980 : }
1981 1725893 : Word16 idiv1616U( Word16 x, Word16 y )
1982 : {
1983 : Word16 sx, sy;
1984 :
1985 : /* make y > x to meet the requirements for div_s parameters */
1986 1725893 : sx = norm_s( x );
1987 1725893 : sy = norm_s( y );
1988 1725893 : x = shl( x, sx );
1989 1725893 : y = shl( y, sy );
1990 :
1991 1725893 : if ( x >= y )
1992 : {
1993 936448 : x = shr( x, 1 );
1994 936448 : sx = sub( sx, 1 );
1995 : }
1996 :
1997 : /* divide and shift */
1998 1725893 : y = shr( div_s( x, y ), sub( 15, sub( sy, sx ) ) );
1999 :
2000 1725893 : return y;
2001 : }
2002 :
2003 :
2004 48158382 : Word16 idiv1616( Word16 x, Word16 y )
2005 : {
2006 : Word16 s, num, den, sign;
2007 :
2008 :
2009 48158382 : num = abs_s( x );
2010 48158382 : den = abs_s( y );
2011 :
2012 48158382 : sign = s_and( s_xor( x, y ), -32768 /* 0x8000 */ );
2013 :
2014 : /* make num > den */
2015 48158382 : s = add( sub( norm_s( den ), norm_s( num ) ), 1 );
2016 48158382 : s = s_max( s, 0 );
2017 :
2018 48158382 : den = shl( den, s );
2019 :
2020 : /* divide and shift */
2021 48158382 : y = shr( div_s( num, den ), sub( 15, s ) );
2022 :
2023 48158382 : if ( 0 != sign )
2024 : {
2025 139143 : y = negate( y );
2026 : }
2027 :
2028 :
2029 48158382 : return y;
2030 : }
2031 :
2032 0 : Word16 idiv1616_1( Word16 x, Word16 y )
2033 : {
2034 0 : IF( L_mult0( x, y ) < 0 )
2035 : {
2036 0 : return negate( idiv1616( abs_s( x ), abs_s( y ) ) );
2037 : }
2038 0 : ELSE IF( L_mult0( x, y ) > 0 )
2039 : {
2040 0 : return idiv1616( x, y );
2041 : }
2042 : ELSE
2043 : {
2044 0 : return 0;
2045 : }
2046 : }
2047 :
2048 11156181 : Word32 norm_llQ31( /* o : normalized result Q31 */
2049 : Word32 L_c, /* i : upper bits of accu Q-1 */
2050 : Word32 L_sum, /* i : lower bits of accu, unsigned Q31 */
2051 : Word16 *exp /* o : exponent of result in [-32,31] Q0 */
2052 : )
2053 : {
2054 : Word16 i;
2055 : Word32 L_tmp;
2056 : #ifdef BASOP_NOGLOB_DECLARE_LOCAL
2057 11156181 : Flag Overflow = 0;
2058 11156181 : Flag Carry = 0;
2059 : #endif /* BASOP_NOGLOB */
2060 :
2061 : /* Move MSBit of L_sum into L_c */
2062 11156181 : Carry = 0;
2063 11156181 : L_tmp = L_add_co( L_sum, L_sum, &Carry, &Overflow ); /* L_tmp = L_sum << 1 */
2064 11156181 : L_c = L_add_co( L_c, L_c, &Carry, &Overflow );
2065 11156181 : L_add( 0, 0 );
2066 11156181 : test();
2067 11156181 : IF( ( L_c != (Word32) 0L ) && ( L_c != (Word32) 0xFFFFFFFFL ) )
2068 : {
2069 3253333 : i = norm_l( L_c );
2070 3253333 : L_c = L_shl( L_c, i );
2071 3253333 : i = sub( 31, i ); /* positive exponent */
2072 3253333 : L_sum = L_lshr( L_tmp, 1 ); /* L_sum with MSBit=0 */
2073 3253333 : L_sum = L_lshr( L_sum, i );
2074 3253333 : L_sum = L_add( L_c, L_sum );
2075 : }
2076 : ELSE
2077 : {
2078 7902848 : i = -32;
2079 7902848 : move16(); /* default exponent, if total sum=0 */
2080 7902848 : IF( L_sum )
2081 : {
2082 2581784 : i = norm_l( L_sum );
2083 2581784 : L_sum = L_shl( L_sum, i );
2084 2581784 : i = negate( i ); /* negative or zero exponent */
2085 : }
2086 : }
2087 11156181 : *exp = i;
2088 11156181 : move16();
2089 11156181 : return L_sum;
2090 : }
2091 :
2092 : Word32 w_norm_llQ31( Word64 L_sum, Word16 *exp );
2093 422217 : Word32 w_norm_llQ31( /* o : normalized result Q31 */
2094 : Word64 L_sum, /* i : upper and lower bits of accu, unsigned Q31 */
2095 : Word16 *exp /* o : exponent of result in [-32,31] Q0 */
2096 : )
2097 : {
2098 : Word32 L_tmp;
2099 : Word16 exp_val;
2100 422217 : Word64 L64_inp64 = L_sum;
2101 422217 : move64();
2102 :
2103 422217 : L64_inp64 = W_shl( L64_inp64, 1 );
2104 422217 : exp_val = W_norm( L64_inp64 );
2105 422217 : L64_inp64 = W_shl( L64_inp64, exp_val );
2106 422217 : exp_val = sub( 31, exp_val );
2107 422217 : if ( EQ_64( L_sum, 0 ) )
2108 : {
2109 49 : exp_val = -32;
2110 49 : move16();
2111 : }
2112 422217 : *exp = exp_val;
2113 422217 : move16();
2114 :
2115 422217 : L_tmp = W_extract_h( L64_inp64 );
2116 422217 : return L_tmp;
2117 : }
2118 :
2119 29237 : Word32 Dot_product16HQ( /* o : normalized result Q31 */
2120 : const Word32 L_off, /* i : initial sum value Qn */
2121 : const Word16 x[], /* i : x vector Qn */
2122 : const Word16 y[], /* i : y vector Qn */
2123 : const Word16 lg, /* i : vector length, range [0..7FFF] Q0 */
2124 : Word16 *exp /* o : exponent of result in [-32,31] Q0 */
2125 : )
2126 : {
2127 : Word16 i;
2128 : Word32 L_sum;
2129 : Word64 L_sum64;
2130 :
2131 29237 : L_sum64 = W_deposit32_l( L_off );
2132 :
2133 1234701 : FOR( i = 0; i < lg; i++ )
2134 : {
2135 1205464 : L_sum64 = W_mac_16_16( L_sum64, x[i], y[i] );
2136 : }
2137 :
2138 29237 : L_sum = w_norm_llQ31( L_sum64, exp );
2139 29237 : return L_sum;
2140 : }
2141 :
2142 1602 : Word32 Norm32Norm( const Word32 *x, const Word16 headroom, const Word16 length, Word16 *result_e )
2143 : {
2144 : Word32 L_tmp, L_tmp2;
2145 : Word16 i, shift, tmp;
2146 :
2147 1602 : move16();
2148 1602 : shift = headroom;
2149 :
2150 1602 : L_tmp = L_deposit_l( 0 );
2151 :
2152 193372 : FOR( i = 0; i < length; i++ )
2153 : {
2154 191770 : L_tmp2 = L_sub( L_tmp, 0x40000000 );
2155 191770 : if ( L_tmp2 >= 0 )
2156 3223 : shift = sub( shift, 1 );
2157 191770 : if ( L_tmp2 >= 0 )
2158 3223 : L_tmp = L_shr( L_tmp, 2 );
2159 :
2160 191770 : tmp = round_fx_sat( L_shl_sat( x[i], shift ) );
2161 191770 : L_tmp = L_mac0_sat( L_tmp, tmp, tmp ); /* exponent = (1-shift*2) , Q(30+shift*2) */
2162 : }
2163 :
2164 1602 : move16();
2165 1602 : *result_e = sub( 1, shl( shift, 1 ) );
2166 :
2167 1602 : return L_tmp;
2168 : }
2169 :
2170 392980 : Word32 Dot_productSq16HQ( /* o : normalized result Q31 */
2171 : const Word32 L_off, /* i : initial sum value Qn */
2172 : const Word16 x[], /* i : x vector Qn */
2173 : const Word16 lg, /* i : vector length, range [0..7FFF] Q0 */
2174 : Word16 *exp /* o : exponent of result in [-32,31] Q0 */
2175 : )
2176 : {
2177 : Word16 i;
2178 : Word32 L_sum;
2179 : Word64 L_sum64;
2180 :
2181 392980 : L_sum64 = W_deposit32_l( L_off );
2182 :
2183 43957661 : FOR( i = 0; i < lg; i++ )
2184 : {
2185 43564681 : L_sum64 = W_mac_16_16( L_sum64, x[i], x[i] );
2186 : }
2187 392980 : L_sum = w_norm_llQ31( L_sum64, exp );
2188 :
2189 392980 : return L_sum;
2190 : }
2191 :
2192 118831 : Word32 dotp_s_fx( const Word16 *x, const Word16 *y, const Word16 n, Word16 s )
2193 : {
2194 : Word16 i;
2195 : Word16 n2;
2196 : Word32 L_tmp;
2197 : Word32 L_sum;
2198 :
2199 :
2200 118831 : L_sum = 0;
2201 118831 : move32();
2202 :
2203 118831 : n2 = shr( n, 1 );
2204 :
2205 118831 : s = sub( s, 1 );
2206 :
2207 3738823 : FOR( i = 0; i < n2; i++ )
2208 : {
2209 3619992 : L_tmp = L_mult0( x[2 * i], y[2 * i] );
2210 3619992 : L_tmp = L_mac0( L_tmp, x[2 * i + 1], y[2 * i + 1] );
2211 3619992 : L_sum = L_add( L_sum, L_shr( L_tmp, s ) );
2212 : }
2213 :
2214 118831 : IF( s_and( n, 1 ) )
2215 : {
2216 63769 : L_tmp = L_mult0( x[n - 1], y[n - 1] );
2217 63769 : L_sum = L_add( L_sum, L_shr( L_tmp, s ) );
2218 : }
2219 :
2220 :
2221 118831 : return L_sum;
2222 : }
2223 :
2224 :
2225 52943412 : Word32 BASOP_util_Pow2(
2226 : const Word32 exp_m,
2227 : const Word16 exp_e,
2228 : Word16 *result_e )
2229 : {
2230 : static const Word16 pow2Coeff[8] = {
2231 : 22713 /*0.693147180559945309417232121458177 Q15*/, /* ln(2)^1 /1! */
2232 : 7872 /*0.240226506959100712333551263163332 Q15*/, /* ln(2)^2 /2! */
2233 : 1819 /*0.0555041086648215799531422637686218 Q15*/, /* ln(2)^3 /3! */
2234 : 315 /*0.00961812910762847716197907157365887 Q15*/, /* ln(2)^4 /4! */
2235 : 44 /*0.00133335581464284434234122219879962 Q15*/, /* ln(2)^5 /5! */
2236 : 5 /*1.54035303933816099544370973327423e-4 Q15*/, /* ln(2)^6 /6! */
2237 : 0 /*1.52527338040598402800254390120096e-5 Q15*/, /* ln(2)^7 /7! */
2238 : 0 /*1.32154867901443094884037582282884e-6 Q15*/ /* ln(2)^8 /8! */
2239 : };
2240 :
2241 : Word32 frac_part, tmp_frac, result_m;
2242 : Word16 int_part;
2243 :
2244 52943412 : int_part = 0; /* to avoid compilation warnings */
2245 52943412 : frac_part = 0; /* to avoid compilation warnings */
2246 :
2247 52943412 : IF( exp_e > 0 )
2248 : {
2249 : /* "+ 1" compensates L_shr(,1) of the polynomial evaluation at the loop end. */
2250 :
2251 52897575 : int_part = add( 1, extract_l( L_shr( exp_m, sub( 31, exp_e ) ) ) );
2252 52897575 : frac_part = L_lshl( exp_m, exp_e );
2253 52897575 : frac_part = L_and( 0x7FFFFFFF, frac_part );
2254 : }
2255 52943412 : if ( exp_e <= 0 )
2256 45837 : frac_part = L_shl( exp_m, exp_e );
2257 52943412 : if ( exp_e <= 0 )
2258 : {
2259 45837 : int_part = 1;
2260 45837 : move16();
2261 : }
2262 :
2263 : /* Best accuracy is around 0, so try to get there with the fractional part. */
2264 52943412 : IF( ( tmp_frac = L_sub( frac_part, 1073741824l /*0.5 Q31*/ ) ) >= 0 )
2265 : {
2266 26130934 : int_part = add( int_part, 1 );
2267 26130934 : frac_part = L_sub( tmp_frac, 1073741824l /*0.5 Q31*/ );
2268 : }
2269 26812478 : ELSE IF( ( tmp_frac = L_add( frac_part, 1073741824l /*0.5 Q31*/ ) ) < 0 )
2270 : {
2271 0 : int_part = sub( int_part, 1 );
2272 0 : frac_part = L_add( tmp_frac, 1073741824l /*0.5 Q31*/ );
2273 : }
2274 :
2275 : /* Evaluate taylor polynomial which approximates 2^x */
2276 : {
2277 : Word32 p;
2278 : Word16 i;
2279 :
2280 :
2281 : /* First taylor series coefficient a_0 = 1.0, scaled by 0.5 due to L_shr(,1). */
2282 52943412 : result_m = L_add( 1073741824l /*1.0/2.0 Q31*/, L_shr( Mpy_32_16_1( frac_part, pow2Coeff[0] ), 1 ) );
2283 52943412 : p = Mpy_32_32( frac_part, frac_part );
2284 370603884 : FOR( i = 1; i < 7; i++ )
2285 : {
2286 : /* next taylor series term: a_i * x^i, x=0 */
2287 317660472 : result_m = L_add( result_m, L_shr( Mpy_32_16_1( p, pow2Coeff[i] ), 1 ) );
2288 317660472 : p = Mpy_32_32( p, frac_part );
2289 : }
2290 52943412 : result_m = L_add( result_m, L_shr( Mpy_32_16_1( p, pow2Coeff[i] ), 1 ) );
2291 : }
2292 52943412 : *result_e = int_part;
2293 52943412 : move16();
2294 52943412 : return result_m;
2295 : }
2296 :
2297 0 : Word16 findIndexOfMaxWord32( Word32 *x, const Word16 len )
2298 : {
2299 : Word16 i, indx;
2300 :
2301 :
2302 0 : indx = 0;
2303 0 : move16();
2304 0 : FOR( i = 1; i < len; i++ )
2305 : {
2306 0 : if ( GT_32( x[i], x[indx] ) )
2307 : {
2308 0 : indx = i;
2309 0 : move16();
2310 : }
2311 : }
2312 :
2313 :
2314 0 : return indx;
2315 : }
2316 :
2317 50830728 : Word16 getNormReciprocalWord16( Word16 x )
2318 : {
2319 :
2320 50830728 : assert( x < (Word16) ( sizeof( BASOP_util_normReciprocal ) / sizeof( BASOP_util_normReciprocal[0] ) ) );
2321 :
2322 50830728 : return extract_h( BASOP_util_normReciprocal[x] );
2323 : }
2324 68 : Word16 getNormReciprocalWord16Scale( Word16 x, Word16 s )
2325 : {
2326 :
2327 68 : assert( x < (Word16) ( sizeof( BASOP_util_normReciprocal ) / sizeof( BASOP_util_normReciprocal[0] ) ) );
2328 :
2329 68 : return round_fx( L_shl( BASOP_util_normReciprocal[x], s ) );
2330 : }
2331 :
2332 :
2333 : /*! r: result of division x/y, not normalized */
2334 315048296 : Word16 BASOP_Util_Divide3216_Scale(
2335 : Word32 x, /* i : numerator, signed */
2336 : Word16 y, /* i : denominator, signed */
2337 : Word16 *s ) /* o : scaling, 0, if x==0 */
2338 : {
2339 : Word16 z;
2340 : Word16 sx;
2341 : Word16 sy;
2342 : Word16 sign;
2343 :
2344 : /*assert (x > (Word32)0);
2345 : assert (y >= (Word16)0);*/
2346 :
2347 : /* check, if numerator equals zero, return zero then */
2348 315048296 : IF( x == (Word32) 0 )
2349 : {
2350 27281260 : *s = 0;
2351 27281260 : move16();
2352 :
2353 27281260 : return ( (Word16) 0 );
2354 : }
2355 :
2356 287767036 : sign = s_xor( extract_h( x ), y ); /* just to exor the sign bits */
2357 : BASOP_SATURATE_WARNING_OFF
2358 287767036 : x = L_abs( x );
2359 287767036 : y = abs_s( y );
2360 : BASOP_SATURATE_WARNING_ON
2361 287767036 : sx = sub( norm_l( x ), 1 );
2362 287767036 : x = L_shl( x, sx );
2363 287767036 : sy = norm_s( y );
2364 287767036 : y = shl( y, sy );
2365 287767036 : *s = sub( sy, sx );
2366 287767036 : move16();
2367 :
2368 287767036 : z = div_s( round_fx( x ), y );
2369 :
2370 287767036 : if ( sign < 0 ) /* if sign bits differ, negate the result */
2371 : {
2372 3369498 : z = negate( z );
2373 : }
2374 :
2375 287767036 : return z;
2376 : }
2377 :
2378 :
2379 : /*************************************************************************
2380 : *
2381 : * FUNCTION: Log2_norm()
2382 : *
2383 : * PURPOSE: Computes log2(L_x, exp), where L_x is positive and
2384 : * normalized, and exp is the normalisation exponent
2385 : * If L_x is negative or zero, the result is 0.
2386 : *
2387 : * DESCRIPTION:
2388 : * The function Log2(L_x) is approximated by a table and linear
2389 : * interpolation. The following steps are used to compute Log2(L_x)
2390 : *
2391 : * 1- exponent = 30-norm_exponent
2392 : * 2- i = bit25-b31 of L_x; 32<=i<=63 (because of normalization).
2393 : * 3- a = bit10-b24
2394 : * 4- i -=32
2395 : * 5- fraction = table[i]<<16 - (table[i] - table[i+1]) * a * 2
2396 : *
2397 : *************************************************************************/
2398 :
2399 27006 : Word32 BASOP_Util_fPow(
2400 : Word32 base_m,
2401 : Word16 base_e,
2402 : Word32 exp_m,
2403 : Word16 exp_e,
2404 : Word16 *result_e )
2405 : {
2406 :
2407 : Word16 ans_lg2_e, base_lg2_e;
2408 : Word32 base_lg2_m, ans_lg2_m, result_m;
2409 : Word16 shift;
2410 :
2411 27006 : test();
2412 27006 : IF( ( base_m == 0 ) && ( exp_m != 0 ) )
2413 : {
2414 0 : *result_e = 0;
2415 0 : move16();
2416 0 : return 0;
2417 : }
2418 : /* Calc log2 of base */
2419 27006 : shift = norm_l( base_m );
2420 27006 : base_m = L_shl( base_m, shift );
2421 27006 : base_e = sub( base_e, shift );
2422 27006 : base_lg2_m = BASOP_Util_Log2( base_m );
2423 :
2424 : /* shift: max left shift such that neither base_e or base_lg2_m saturate. */
2425 27006 : shift = sub( s_min( norm_s( base_e ), WORD16_BITS - 1 - LD_DATA_SCALE ), 1 );
2426 : /* Compensate shift into exponent of result. */
2427 27006 : base_lg2_e = sub( WORD16_BITS - 1, shift );
2428 27006 : base_lg2_m = L_add( L_shr( base_lg2_m, sub( WORD16_BITS - 1 - LD_DATA_SCALE, shift ) ), L_deposit_h( shl( base_e, shift ) ) );
2429 :
2430 : /* Prepare exp */
2431 27006 : shift = norm_l( exp_m );
2432 27006 : exp_m = L_shl( exp_m, shift );
2433 27006 : exp_e = sub( exp_e, shift );
2434 :
2435 : /* Calc base pow exp */
2436 27006 : ans_lg2_m = Mpy_32_32( base_lg2_m, exp_m );
2437 27006 : ans_lg2_e = add( exp_e, base_lg2_e );
2438 :
2439 : /* Calc antilog */
2440 27006 : result_m = BASOP_util_Pow2( ans_lg2_m, ans_lg2_e, result_e );
2441 :
2442 27006 : return result_m;
2443 : }
2444 :
2445 : /*___________________________________________________________________________
2446 : | |
2447 : | Function Name : Dot_product12_offs() |
2448 : | |
2449 : | Compute scalar product of <x[],y[]> using accumulator. |
2450 : | The parameter 'L_off' is added to the accumulation result. |
2451 : | The result is normalized (in Q31) with exponent (0..30). |
2452 : | Notes: |
2453 : | o data in x[],y[] must provide enough headroom for accumulation |
2454 : | o L_off must correspond in format with product of x,y |
2455 : | Example: 0.01f for Q9 x Q9: 0x0000147B in Q19 |
2456 : | means: L_off = FL2WORD32_SCALE(0.01,31-19) |
2457 : |---------------------------------------------------------------------------|
2458 : | Algorithm: |
2459 : | |
2460 : | dot_product = L_off + sum(x[i]*y[i]) i=0..N-1 |
2461 : |___________________________________________________________________________|
2462 : */
2463 18564 : Word32 Dot_product12_offs( /* (o) Q31: normalized result (1 < val <= -1) */
2464 : const Word16 x[], /* (i) 12bits: x vector */
2465 : const Word16 y[], /* (i) 12bits: y vector */
2466 : const Word16 lg, /* (i) : vector length in range [1..256] */
2467 : Word16 *exp, /* (o) : exponent of result (0..+30) */
2468 : Word32 L_off /* (i) initial summation offset / 2 */
2469 : )
2470 : {
2471 : Word16 i, sft;
2472 : Word32 L_sum;
2473 :
2474 : Word64 L_sum64;
2475 :
2476 18564 : L_sum64 = W_deposit32_l( L_off );
2477 1698212 : FOR( i = 0; i < lg; i++ )
2478 : {
2479 1679648 : L_sum64 = W_mac0_16_16( L_sum64, x[i], y[i] );
2480 : }
2481 18564 : L_sum = W_sat_l( L_sum64 );
2482 :
2483 : /* Normalize acc in Q31 */
2484 :
2485 18564 : sft = norm_l( L_sum );
2486 18564 : if ( exp != NULL )
2487 : {
2488 18564 : L_sum = L_shl( L_sum, sft );
2489 : }
2490 :
2491 : /* exponent = 0..30, when L_sum != 0 */
2492 18564 : if ( L_sum != 0 )
2493 : {
2494 18564 : sft = sub( 31, sft );
2495 : }
2496 :
2497 18564 : if ( exp != NULL )
2498 : {
2499 18564 : *exp = sft;
2500 18564 : move16();
2501 : }
2502 :
2503 18564 : return L_sum;
2504 : }
2505 :
2506 0 : Word32 Dot_product15_offs( /* (o) Q31: normalized result (1 < val <= -1) */
2507 : const Word16 x[], /* (i) 15bits: x vector */
2508 : const Word16 y[], /* (i) 15bits: y vector */
2509 : const Word16 lg, /* (i) : vector length in range [1..256] */
2510 : Word16 *exp, /* (o) : exponent of result (0..+30) */
2511 : Word32 L_off /* (i) initial summation offset */
2512 : )
2513 : {
2514 : Word16 i, sft, fac, ld;
2515 : Word32 L_sum;
2516 :
2517 0 : ld = sub( 14, norm_s( lg ) );
2518 0 : fac = shr( -32768, ld );
2519 0 : L_sum = L_shr( L_off, ld );
2520 :
2521 0 : FOR( i = 0; i < lg; i++ )
2522 : {
2523 0 : L_sum = L_add( L_sum, Mpy_32_16_1( L_msu( 0, y[i], fac ), x[i] ) );
2524 : }
2525 :
2526 : /* Avoid returning 0 */
2527 0 : if ( L_sum == 0 )
2528 : {
2529 0 : L_sum = L_add( L_sum, 1 );
2530 : }
2531 :
2532 : /* Normalize acc in Q31 */
2533 0 : sft = norm_l( L_sum );
2534 0 : L_sum = L_shl( L_sum, sft );
2535 :
2536 : /* exponent = 0..30, when L_sum != 0 */
2537 0 : if ( L_sum != 0 )
2538 : {
2539 0 : sft = add( ld, sub( 30, sft ) );
2540 : }
2541 :
2542 0 : *exp = sft;
2543 0 : move16();
2544 :
2545 0 : return L_sum;
2546 : }
2547 :
2548 169878888 : Word16 BASOP_Util_Cmp_Mant32Exp /*!< o: flag: result of comparison */
2549 : /* 0, if a == b */
2550 : /* 1, if a > b */
2551 : /* -1, if a < b */
2552 : ( Word32 a_m, /*!< i: Mantissa of 1st operand a */
2553 : Word16 a_e, /*!< i: Exponent of 1st operand a */
2554 : Word32 b_m, /*!< i: Mantissa of 2nd operand b */
2555 : Word16 b_e ) /*!< i: Exponent of 2nd operand b */
2556 :
2557 : {
2558 : Word32 diff_m;
2559 : Word16 diff_e, shift, result;
2560 :
2561 :
2562 : /*
2563 : This function compares two input parameters, both represented by a 32-bit mantissa and a 16-bit exponent.
2564 : If both values are identical, 0 is returned.
2565 : If a is greater b, 1 is returned.
2566 : If a is less than b, -1 is returned.
2567 : */
2568 :
2569 : /* Check, if both mantissa and exponents are identical, when normalized: return 0 */
2570 169878888 : shift = norm_l( a_m );
2571 169878888 : if ( shift )
2572 128455306 : a_m = L_shl( a_m, shift );
2573 169878888 : if ( shift )
2574 128455306 : a_e = sub( a_e, shift );
2575 :
2576 169878888 : shift = norm_l( b_m );
2577 169878888 : if ( shift )
2578 144099600 : b_m = L_shl( b_m, shift );
2579 169878888 : if ( shift )
2580 144099600 : b_e = sub( b_e, shift );
2581 :
2582 : /* align exponent, if any mantissa is zero */
2583 169878888 : if ( !a_m )
2584 : {
2585 13480915 : a_e = b_e;
2586 13480915 : move16();
2587 : }
2588 169878888 : if ( !b_m )
2589 : {
2590 12510787 : b_e = a_e;
2591 12510787 : move16();
2592 : }
2593 :
2594 : BASOP_SATURATE_WARNING_OFF_EVS
2595 169878888 : diff_m = L_sub_sat( a_m, b_m );
2596 : BASOP_SATURATE_WARNING_ON_EVS
2597 169878888 : diff_e = sub( a_e, b_e );
2598 :
2599 169878888 : test();
2600 169878888 : IF( diff_m == 0 && diff_e == 0 )
2601 : {
2602 8640313 : return 0;
2603 : }
2604 :
2605 : /* Check sign, exponent and mantissa to identify, whether a is greater b or not */
2606 161238575 : result = sub( 0, 1 );
2607 :
2608 161238575 : IF( a_m >= 0 )
2609 : {
2610 : /* a is positive */
2611 160972115 : if ( b_m < 0 )
2612 : {
2613 0 : result = 1;
2614 0 : move16();
2615 : }
2616 :
2617 160972115 : test();
2618 160972115 : test();
2619 160972115 : test();
2620 160972115 : if ( ( b_m >= 0 ) && ( ( diff_e > 0 ) || ( diff_e == 0 && diff_m > 0 ) ) )
2621 : {
2622 59816337 : result = 1;
2623 59816337 : move16();
2624 : }
2625 : }
2626 : ELSE
2627 : {
2628 : /* a is negative */
2629 266460 : test();
2630 266460 : test();
2631 266460 : test();
2632 266460 : if ( ( b_m < 0 ) && ( ( diff_e < 0 ) || ( diff_e == 0 && diff_m > 0 ) ) )
2633 : {
2634 264460 : result = 1;
2635 264460 : move16();
2636 : }
2637 : }
2638 161238575 : return result;
2639 : }
2640 :
2641 : /*
2642 :
2643 : headroom is introduced into acc
2644 : */
2645 :
2646 :
2647 6398491260 : Word32 BASOP_Util_Add_Mant32Exp /* o : normalized result mantissa */
2648 : ( Word32 a_m, /* i : Mantissa of 1st operand a */
2649 : Word16 a_e, /* i : Exponent of 1st operand a */
2650 : Word32 b_m, /* i : Mantissa of 2nd operand b */
2651 : Word16 b_e, /* i : Exponent of 2nd operand b */
2652 : Word16 *ptr_e ) /* o : exponent of result */
2653 : {
2654 : Word32 L_tmp;
2655 : Word16 shift;
2656 :
2657 : /* Compare exponents: the difference is limited to +/- 30
2658 : The Word32 mantissa of the operand with lower exponent is shifted right by the exponent difference.
2659 : Then, the unshifted mantissa of the operand with the higher exponent is added. The addition result
2660 : is normalized and the result represents the mantissa to return. The returned exponent takes into
2661 : account all shift operations.
2662 : */
2663 :
2664 6398491260 : if ( !a_m )
2665 : {
2666 1321490120 : a_e = b_e;
2667 1321490120 : move16();
2668 : }
2669 :
2670 6398491260 : if ( !b_m )
2671 : {
2672 428987546 : b_e = a_e;
2673 428987546 : move16();
2674 : }
2675 :
2676 6398491260 : shift = sub( a_e, b_e );
2677 6398491260 : shift = s_max( -31, shift );
2678 6398491260 : shift = s_min( 31, shift );
2679 6398491260 : if ( shift < 0 )
2680 : {
2681 : /* exponent of b is greater than exponent of a, shr a_m */
2682 3402275556 : a_m = L_shl( a_m, shift );
2683 : }
2684 6398491260 : if ( shift > 0 )
2685 : {
2686 : /* exponent of a is greater than exponent of b */
2687 1404964329 : b_m = L_shr( b_m, shift );
2688 : }
2689 6398491260 : a_e = add( s_max( a_e, b_e ), 1 );
2690 6398491260 : L_tmp = L_add( L_shr( a_m, 1 ), L_shr( b_m, 1 ) );
2691 6398491260 : shift = norm_l( L_tmp );
2692 6398491260 : if ( shift )
2693 5602641504 : L_tmp = L_shl( L_tmp, shift );
2694 6398491260 : if ( L_tmp == 0 )
2695 : {
2696 612632762 : a_e = 0;
2697 612632762 : move16();
2698 : }
2699 6398491260 : if ( L_tmp != 0 )
2700 5785858498 : a_e = sub( a_e, shift );
2701 6398491260 : *ptr_e = a_e;
2702 :
2703 6398491260 : return ( L_tmp );
2704 : }
2705 :
2706 :
2707 : static const Word16 shift_lc[] = { 9, 10 };
2708 :
2709 0 : Word32 Isqrt_lc1(
2710 : Word32 frac, /* i : Q31: normalized value (1.0 < frac <= 0.5) */
2711 : Word16 *exp /* i/o: exponent (value = frac x 2^exponent) */
2712 : )
2713 : {
2714 : Word16 i, a;
2715 : Word32 L_tmp;
2716 :
2717 0 : IF( frac <= (Word32) 0 )
2718 : {
2719 0 : *exp = 0;
2720 0 : move16();
2721 0 : return 0x7fffffff; /*0x7fffffff*/
2722 : }
2723 :
2724 : /* If exponant odd -> shift right by 10 (otherwise 9) */
2725 0 : L_tmp = L_shr( frac, shift_lc[s_and( *exp, 1 )] );
2726 :
2727 : /* 1) -16384 to shift left and change sign */
2728 : /* 2) 32768 to Add 1 to Exponent like it was divided by 2 */
2729 : /* 3) We let the mac_r add another 0.5 because it imitates */
2730 : /* the behavior of shr on negative number that should */
2731 : /* not be rounded towards negative infinity. */
2732 : /* It replaces: */
2733 : /* *exp = negate(shr(sub(*exp, 1), 1)); move16(); */
2734 0 : *exp = mac_r( 32768, *exp, -16384 );
2735 0 : move16();
2736 :
2737 0 : a = extract_l( L_tmp ); /* Extract b10-b24 */
2738 0 : a = lshr( a, 1 );
2739 :
2740 0 : i = mac_r( L_tmp, -16 * 2 - 1, 16384 ); /* Extract b25-b31 minus 16 */
2741 :
2742 0 : L_tmp = L_msu( L_table_isqrt[i], table_isqrt_diff[i], a ); /* table[i] << 16 - diff*a*2 */
2743 :
2744 0 : return L_tmp;
2745 : }
2746 :
2747 :
2748 10742 : void bufferCopyFx(
2749 : Word16 *src, /*<! Qx pointer to input buffer */
2750 : Word16 *dest, /*<! Qx pointer to output buffer */
2751 : Word16 length, /*<! Q0 length of buffer to copy */
2752 : Word16 Qf_src, /*<! Q0 Q format (frac-bits) of source buffer */
2753 : Word16 Qf_dest, /*<! Q0 Q format (frac-bits) of dest buffer */
2754 : Word16 Q_src, /*<! Q0 exponent of source buffer */
2755 : Word16 Q_dest /*<! Q0 exponent of destination buffer */
2756 : )
2757 : {
2758 : Word16 tmp_16, i;
2759 :
2760 : /*Copy( st->old_exc, exc_buf, st->old_exc_len);*/
2761 10742 : tmp_16 = sub( sub( Qf_src, Qf_dest ), sub( Q_src, Q_dest ) );
2762 10742 : IF( tmp_16 > 0 ) /*if value will be shifted right, do a multiplication with rounding ->preserves more accuracy*/
2763 : {
2764 2968 : tmp_16 = shl( 1, sub( 15, tmp_16 ) );
2765 965432 : FOR( i = 0; i < length; i++ )
2766 : {
2767 962464 : *( dest + i ) = mult_r( *( src + i ), tmp_16 );
2768 962464 : move16();
2769 : }
2770 : }
2771 7774 : ELSE IF( tmp_16 < 0 ) /*leftshift - no accuracy preservation needed*/
2772 : {
2773 162845 : FOR( i = 0; i < length; i++ )
2774 : {
2775 162175 : *( dest + i ) = shr_sat( *( src + i ), tmp_16 );
2776 162175 : move16();
2777 : }
2778 : }
2779 : ELSE /*no shift, simply copy*/
2780 : {
2781 2584459 : FOR( i = 0; i < length; i++ )
2782 : {
2783 2577355 : *( dest + i ) = *( src + i );
2784 2577355 : move16();
2785 : }
2786 : }
2787 10742 : }
2788 :
2789 38198 : Word32 dotWord32_16_Mant32Exp( const Word32 *bufX32, /* i: 32-bit buffer with unknown headroom */
2790 : Word16 bufX32_exp, /* i: exponent of buffer bufX32 */
2791 : const Word16 *bufY16, /* i: 16-bit buffer quite right-aligned */
2792 : Word16 bufY16_exp, /* i: exponent of buffer bufY16 */
2793 : Word16 len, /* i: buffer len to process */
2794 : Word16 *exp ) /* o: result exponent */
2795 : {
2796 : Word32 L_sum;
2797 : Word16 shift, shift1, i;
2798 :
2799 :
2800 38198 : shift = getScaleFactor32( bufX32, len ); /* current available headroom */
2801 38198 : shift = sub( shift, sub( 14, norm_s( len ) ) ); /* reduced required headroom */
2802 38198 : L_sum = 0; /* Clear accu */
2803 38198 : move32();
2804 2387887 : FOR( i = 0; i < len; i++ )
2805 : {
2806 2349689 : L_sum = L_mac0( L_sum, round_fx( L_shl( bufX32[i], shift ) ), bufY16[i] );
2807 : }
2808 38198 : shift1 = norm_l( L_sum );
2809 38198 : L_sum = L_shl( L_sum, shift1 ); /* return value */
2810 :
2811 38198 : shift = sub( add( bufX32_exp, bufY16_exp ), add( shift, shift1 ) );
2812 38198 : shift = add( shift, 1 ); /* compensate for factor of 2 introduced by L_mac0 */
2813 : /* In case of NULL result, we want to have a 0 exponent too */
2814 38198 : if ( L_sum == 0 )
2815 562 : shift = 0;
2816 38198 : *exp = shift;
2817 38198 : move16();
2818 :
2819 :
2820 38198 : return L_sum;
2821 : }
2822 :
2823 566 : Word16 BASOP_Util_lin2dB( Word32 x, Word16 x_e, Word16 fEnergy )
2824 : {
2825 566 : assert( x >= 0 );
2826 :
2827 : /* log2 */
2828 566 : x = L_shr( BASOP_Util_Log2( x ), 25 - 16 ); /* Q16 */
2829 :
2830 : /* add exponent */
2831 566 : x = L_msu( x, x_e, -32768 /* 0x8000 */ );
2832 :
2833 : /* convert log2 to 20*log10 */
2834 566 : x = Mpy_32_16_1( x, 24660 /*6.0206f Q12*/ ); /* Q13 */
2835 :
2836 : /* if energy divide by 2 (->10*log10) */
2837 566 : if ( fEnergy != 0 )
2838 566 : x = L_shr( x, 1 );
2839 :
2840 : /* return dB as 7Q8 */
2841 566 : return round_fx( L_shl( x, 8 - 13 + 16 ) ); /* Q8 */
2842 : }
2843 :
2844 : /* --- fixp_atan() ---- */
2845 : #define Q_ATANINP ( 25 ) /* Input in q25, Output in q14 */
2846 : #define Q_ATANOUT ( 14 )
2847 : #define ATI_SF ( ( 32 - 1 ) - Q_ATANINP ) /* 6 */
2848 : #define ATO_SF ( ( 16 - 1 ) - Q_ATANOUT ) /* 1 ] -pi/2 .. pi/2 [ */
2849 : /* --- fixp_atan2() --- */
2850 : #define Q_ATAN2OUT ( 13 )
2851 : #define AT2O_SF ( ( 16 - 1 ) - Q_ATAN2OUT ) /* 2 ] -pi .. pi ] */
2852 :
2853 :
2854 34022323 : Word16 BASOP_util_atan2( /* o: atan2(y,x) [-pi,pi] Q13 */
2855 : Word32 y, /* i: */
2856 : Word32 x, /* i: */
2857 : Word16 e /* i: exponent difference (exp_y - exp_x) */
2858 : )
2859 : {
2860 : Word16 q;
2861 : Word32 at;
2862 34022323 : Word16 ret = -32768 /*-1.0f Q15*/;
2863 : Word16 sf, sfo, stf;
2864 : Word32 L_sign;
2865 :
2866 34022323 : if ( L_or( y, x ) == 0 )
2867 : {
2868 138838 : return 0;
2869 : }
2870 :
2871 33883485 : IF( x == 0 )
2872 : {
2873 499572 : ret = 12868 /*+EVS_PI/2 Q13*/;
2874 499572 : move16();
2875 499572 : if ( y < 0 )
2876 : {
2877 284735 : ret = negate( ret );
2878 : }
2879 :
2880 499572 : return ret;
2881 : }
2882 :
2883 : /* --- division */
2884 33383913 : L_sign = L_and( L_xor( x, y ), (Word32) 0x80000000 );
2885 :
2886 33383913 : q = 32767 /*1.0f Q15*/; /* y/x = neg/zero = -Inf */
2887 33383913 : sf = 0;
2888 : BASOP_SATURATE_WARNING_OFF_EVS
2889 33383913 : q = BASOP_Util_Divide3232_uu_1616_Scale( L_abs( y ), L_abs( x ), &sf );
2890 : BASOP_SATURATE_WARNING_ON_EVS
2891 :
2892 : BASOP_SATURATE_WARNING_OFF_EVS
2893 33383913 : if ( L_sign < 0 )
2894 16389368 : q = negate( q );
2895 : BASOP_SATURATE_WARNING_ON_EVS
2896 :
2897 33383913 : sfo = add( sf, e );
2898 :
2899 : /* --- atan() */
2900 33383913 : IF( GT_16( sfo, ATI_SF ) )
2901 : {
2902 : /* --- could not calc fixp_atan() here bec of input data out of range
2903 : ==> therefore give back boundary values */
2904 :
2905 1657789 : sfo = s_min( sfo, MAXSFTAB );
2906 :
2907 : /*q = FL2WORD16( 0.0f ); move16();*/
2908 :
2909 1657789 : if ( q > 0 )
2910 : {
2911 133774 : move16();
2912 133774 : q = +f_atan_expand_range[sfo - ATI_SF - 1];
2913 : }
2914 1657789 : if ( q < 0 )
2915 : {
2916 228127 : move16();
2917 228127 : q = -f_atan_expand_range[sfo - ATI_SF - 1];
2918 : }
2919 : }
2920 : ELSE
2921 : {
2922 : /* --- calc of fixp_atan() is possible; input data within range
2923 : ==> set q on fixed scale level as desired from fixp_atan() */
2924 31726124 : stf = sub( sfo, ATI_SF );
2925 :
2926 31726124 : at = L_deposit_h( q );
2927 31726124 : if ( stf < 0 )
2928 31497659 : at = L_shl( at, stf );
2929 :
2930 31726124 : q = BASOP_util_atan( at ); /* ATO_SF*/
2931 : }
2932 :
2933 :
2934 : /* --- atan2() */
2935 :
2936 33383913 : ret = shr( q, ( AT2O_SF - ATO_SF ) ); /* now AT2O_SF for atan2 */
2937 33383913 : IF( x < 0 )
2938 : {
2939 6193827 : if ( y >= 0 )
2940 : {
2941 3213608 : ret = add( ret, 25736 /*EVS_PI Q13*/ );
2942 : }
2943 6193827 : if ( y < 0 )
2944 : {
2945 2980219 : ret = sub( ret, 25736 /* EVS_PI Q13*/ );
2946 : }
2947 : }
2948 :
2949 33383913 : return ret;
2950 : }
2951 :
2952 : /* SNR of fixp_atan() = 56 dB*/
2953 : #define ONEBY3P56 0x26800000 /* 1.0/3.56 in q31*/
2954 : #define P281 0x00026000 /* 0.281 in q19*/
2955 : #define ONEP571 0x6487 /* 1.571 in q14*/
2956 :
2957 31726124 : Word16 BASOP_util_atan( /* o: atan(x) [-pi/2;pi/2] 1Q14 */
2958 : Word32 x /* i: input data (-64;64) 6Q25 */
2959 : )
2960 : {
2961 : Word16 sign, result, exp;
2962 : Word16 res_e;
2963 : Word16 tmp, xx;
2964 :
2965 :
2966 31726124 : sign = 0;
2967 31726124 : move16();
2968 31726124 : if ( x < 0 )
2969 : {
2970 16103537 : sign = 1;
2971 16103537 : move16();
2972 : }
2973 31726124 : x = L_abs( x );
2974 :
2975 : /* calc of arctan */
2976 31726124 : IF( LT_32( x, 1509950l /*0.045f/64.0f Q31*/ ) )
2977 : {
2978 4363373 : result = round_fx( L_shl( x, 5 ) ); /*Q14*/
2979 : /*BASOP_util_atan_16(0.0444059968): max error 0.0000567511, mean 0.000017, abs mean 0.000017*/
2980 : }
2981 : ELSE
2982 27362751 : IF( LT_32( x, ( L_shl( 1, Q_ATANINP ) - 8482560l /*0.00395 Q31*/ ) ) )
2983 : {
2984 13905920 : xx = round_fx( L_shl( x, 6 ) );
2985 13905920 : tmp = mult_r( xx, xx ); /* q15 * q15 - (16-1) = q15*/
2986 13905920 : tmp = mult_r( tmp, 0x1340 ); /* 15 * (ONEBY3P56) q14 - (16-1) = q14*/
2987 13905920 : tmp = add( tmp, 0x4000 ); /*L_shl(1,14) = 524288*/ /* q14 + q14 = q14 */
2988 13905920 : res_e = Q_ATANOUT - 15 + 14 - 16 + 1;
2989 13905920 : move16();
2990 13905920 : if ( GT_16( xx, tmp ) )
2991 : {
2992 2980210 : res_e = add( res_e, 1 );
2993 : }
2994 13905920 : if ( GT_16( xx, tmp ) )
2995 : {
2996 2980210 : xx = shr( xx, 1 );
2997 : }
2998 13905920 : result = div_s( xx, tmp );
2999 13905920 : result = msu_r( 0, result, shl( -32768, res_e ) );
3000 : /*BASOP_util_atan_16(0.7471138239): max error 0.0020029545, mean 0.000715, abs mean 0.000715*/
3001 : }
3002 13456831 : ELSE IF( LT_32( x, 42949673l /*1.28/64.0 Q31*/ ) )
3003 : {
3004 : Word16 delta_fix;
3005 5174445 : Word32 PI_BY_4 = 1686629684l /*3.1415926/4.0 Q31*/ / 2; /* pi/4 in q30 */
3006 :
3007 5174445 : delta_fix = round_fx( L_shl( L_sub( x, 33554432l /*1.0/64.0 Q31*/ ), 5 ) ); /* q30 */
3008 5174445 : result = round_fx( L_sub( L_add( PI_BY_4, L_msu( 0, delta_fix, -16384 ) ), ( L_mult0( delta_fix, delta_fix ) ) ) );
3009 : /* BASOP_Util_fPow(0.7472000122): max error 0.0020237688, mean 0.000026, abs mean 0.000520 */
3010 : }
3011 : ELSE
3012 : {
3013 8282386 : exp = sub( norm_l( x ), 1 );
3014 8282386 : xx = round_fx( L_shl( x, exp ) );
3015 : /* q25+exp * q25+exp - (16-1) = q19+2*exp*/
3016 8282386 : tmp = mac_r( L_shl( P281, shl( exp, 1 ) ), xx, xx ); /* q19+2*exp + q19+2*exp = q19+2*exp*/
3017 8282386 : res_e = norm_s( tmp );
3018 8282386 : result = div_s( xx, shl( tmp, res_e ) );
3019 8282386 : result = shl( result, add( add( Q_ATANOUT - Q_ATANINP /*-exp*/ + 19 /*+2*exp*/ - 16 + 1, res_e ), exp ) );
3020 8282386 : result = sub( ONEP571, result ); /* q14 + q14 = q14*/
3021 : /*BASOP_Util_fPow(1.2799999714): max error 0.0020168927, mean 0.000066, abs mean 0.000072*/
3022 : }
3023 :
3024 31726124 : if ( sign )
3025 : {
3026 16103537 : result = negate( result );
3027 : }
3028 :
3029 31726124 : return ( result );
3030 : }
3031 :
3032 : /* compare two positive normalized 16 bit mantissa/exponent values */
3033 : /* return value: positive if first value greater, negative if second value greater, zero if equal */
3034 537846 : Word16 compMantExp16Unorm( Word16 m1, Word16 e1, Word16 m2, Word16 e2 )
3035 : {
3036 : Word16 tmp;
3037 :
3038 537846 : assert( ( m1 >= 0x4000 ) && ( m2 >= 0x4000 ) ); /* comparisons below work only for normalized mantissas */
3039 :
3040 537846 : tmp = sub( e1, e2 );
3041 537846 : if ( tmp == 0 )
3042 308716 : tmp = sub( m1, m2 );
3043 :
3044 537846 : return tmp;
3045 : }
3046 :
3047 2029753370 : cmplx CL_scale_t( cmplx x, Word16 y )
3048 : {
3049 : cmplx result;
3050 2029753370 : result.re = Mpy_32_16_1( x.re, y );
3051 2029753370 : result.im = Mpy_32_16_1( x.im, y );
3052 : #ifdef WMOPS
3053 : multiCounter[currCounter].Mpy_32_16_1--;
3054 : multiCounter[currCounter].Mpy_32_16_1--;
3055 : multiCounter[currCounter].CL_scale++;
3056 : #endif
3057 2029753370 : return ( result );
3058 : }
3059 :
3060 7455930 : cmplx CL_dscale_t( cmplx x, Word16 y1, Word16 y2 )
3061 : {
3062 : cmplx result;
3063 7455930 : result.re = Mpy_32_16_1( x.re, y1 );
3064 7455930 : result.im = Mpy_32_16_1( x.im, y2 );
3065 : #ifdef WMOPS
3066 : multiCounter[currCounter].Mpy_32_16_1--;
3067 : multiCounter[currCounter].Mpy_32_16_1--;
3068 : multiCounter[currCounter].CL_dscale++;
3069 : #endif /* #ifdef WMOPS */
3070 7455930 : return ( result );
3071 : }
3072 :
3073 52725560 : cmplx CL_mult_32x16( cmplx input, cmplx_s coeff )
3074 : {
3075 : cmplx result;
3076 52725560 : result.re = L_sub( Mpy_32_16_1( input.re, coeff.re ), Mpy_32_16_1( input.im, coeff.im ) );
3077 52725560 : result.im = L_add( Mpy_32_16_1( input.re, coeff.im ), Mpy_32_16_1( input.im, coeff.re ) );
3078 : #ifdef WMOPS
3079 : multiCounter[currCounter].CL_multr_32x16++;
3080 : multiCounter[currCounter].Mpy_32_16_1--;
3081 : multiCounter[currCounter].Mpy_32_16_1--;
3082 : multiCounter[currCounter].Mpy_32_16_1--;
3083 : multiCounter[currCounter].Mpy_32_16_1--;
3084 : multiCounter[currCounter].L_sub--;
3085 : multiCounter[currCounter].L_add--;
3086 : #endif
3087 52725560 : return result;
3088 : }
|
