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