FreeCalypso > hg > gsm-codec-lib

comparison libgsmefr/inv_sqrt.c @ 53:49dd1ac8e75b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
libgsmefr: import most *.c files from ETSI source
author Mychaela Falconia <falcon@freecalypso.org>
date Fri, 25 Nov 2022 16:18:21 +0000
parents
children 2f0828ba0725
comparison
equal deleted inserted replaced
52:988fd7ff514f 53:49dd1ac8e75b
1 /*************************************************************************
2 *
3 * FUNCTION: Inv_sqrt
4 *
5 * PURPOSE: Computes 1/sqrt(L_x), where L_x is positive.
6 * If L_x is negative or zero, the result is 1 (3fff ffff).
7 *
8 * DESCRIPTION:
9 * The function 1/sqrt(L_x) is approximated by a table and linear
10 * interpolation. The inverse square root is computed using the
11 * following steps:
12 * 1- Normalization of L_x.
13 * 2- If (30-exponent) is even then shift right once.
14 * 3- exponent = (30-exponent)/2 +1
15 * 4- i = bit25-b31 of L_x; 16<=i<=63 because of normalization.
16 * 5- a = bit10-b24
17 * 6- i -=16
18 * 7- L_y = table[i]<<16 - (table[i] - table[i+1]) * a * 2
19 * 8- L_y >>= exponent
20 *
21 *************************************************************************/
22
23 #include "typedef.h"
24 #include "basic_op.h"
25 #include "count.h"
26
27 #include "inv_sqrt.tab" /* Table for inv_sqrt() */
28
29 Word32 Inv_sqrt ( /* (o) : output value */
30 Word32 L_x /* (i) : input value */
31 )
32 {
33 Word16 exp, i, a, tmp;
34 Word32 L_y;
35
36 test ();
37 if (L_x <= (Word32) 0)
38 return ((Word32) 0x3fffffffL);
39
40 exp = norm_l (L_x);
41 L_x = L_shl (L_x, exp); /* L_x is normalize */
42
43 exp = sub (30, exp);
44 test (); logic16 ();
45 if ((exp & 1) == 0) /* If exponent even -> shift right */
46 {
47 L_x = L_shr (L_x, 1);
48 }
49 exp = shr (exp, 1);
50 exp = add (exp, 1);
51
52 L_x = L_shr (L_x, 9);
53 i = extract_h (L_x); /* Extract b25-b31 */
54 L_x = L_shr (L_x, 1);
55 a = extract_l (L_x); /* Extract b10-b24 */
56 a = a & (Word16) 0x7fff; logic16 ();
57
58 i = sub (i, 16);
59
60 L_y = L_deposit_h (table[i]); /* table[i] << 16 */
61 tmp = sub (table[i], table[i + 1]); /* table[i] - table[i+1]) */
62 L_y = L_msu (L_y, tmp, a); /* L_y -= tmp*a*2 */
63
64 L_y = L_shr (L_y, exp); /* denormalization */
65
66 return (L_y);
67 }