2 ** Copyright 2003-2010, VisualOn, Inc.
\r
4 ** Licensed under the Apache License, Version 2.0 (the "License");
\r
5 ** you may not use this file except in compliance with the License.
\r
6 ** You may obtain a copy of the License at
\r
8 ** http://www.apache.org/licenses/LICENSE-2.0
\r
10 ** Unless required by applicable law or agreed to in writing, software
\r
11 ** distributed under the License is distributed on an "AS IS" BASIS,
\r
12 ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
\r
13 ** See the License for the specific language governing permissions and
\r
14 ** limitations under the License.
\r
17 /*___________________________________________________________________________
\r
19 | This file contains mathematic operations in fixed point. |
\r
21 | Isqrt() : inverse square root (16 bits precision). |
\r
22 | Pow2() : 2^x (16 bits precision). |
\r
23 | Log2() : log2 (16 bits precision). |
\r
24 | Dot_product() : scalar product of <x[],y[]> |
\r
26 | These operations are not standard double precision operations. |
\r
27 | They are used where low complexity is important and the full 32 bits |
\r
28 | precision is not necessary. For example, the function Div_32() has a |
\r
29 | 24 bits precision which is enough for our purposes. |
\r
31 | In this file, the values use theses representations: |
\r
33 | Word32 L_32 : standard signed 32 bits format |
\r
34 | Word16 hi, lo : L_32 = hi<<16 + lo<<1 (DPF - Double Precision Format) |
\r
35 | Word32 frac, Word16 exp : L_32 = frac << exp-31 (normalised format) |
\r
36 | Word16 int, frac : L_32 = int.frac (fractional format) |
\r
37 |___________________________________________________________________________|
\r
39 #include "typedef.h"
\r
40 #include "basic_op.h"
\r
41 #include "math_op.h"
\r
43 /*___________________________________________________________________________
\r
45 | Function Name : Isqrt |
\r
47 | Compute 1/sqrt(L_x). |
\r
48 | if L_x is negative or zero, result is 1 (7fffffff). |
\r
49 |---------------------------------------------------------------------------|
\r
52 | 1- Normalization of L_x. |
\r
53 | 2- call Isqrt_n(L_x, exponant) |
\r
54 | 3- L_y = L_x << exponant |
\r
55 |___________________________________________________________________________|
\r
57 Word32 Isqrt( /* (o) Q31 : output value (range: 0<=val<1) */
\r
58 Word32 L_x /* (i) Q0 : input value (range: 0<=val<=7fffffff) */
\r
64 L_x = (L_x << exp); /* L_x is normalized */
\r
66 Isqrt_n(&L_x, &exp);
\r
67 L_y = (L_x << exp); /* denormalization */
\r
71 /*___________________________________________________________________________
\r
73 | Function Name : Isqrt_n |
\r
75 | Compute 1/sqrt(value). |
\r
76 | if value is negative or zero, result is 1 (frac=7fffffff, exp=0). |
\r
77 |---------------------------------------------------------------------------|
\r
80 | The function 1/sqrt(value) is approximated by a table and linear |
\r
83 | 1- If exponant is odd then shift fraction right once. |
\r
84 | 2- exponant = -((exponant-1)>>1) |
\r
85 | 3- i = bit25-b30 of fraction, 16 <= i <= 63 ->because of normalization. |
\r
86 | 4- a = bit10-b24 |
\r
88 | 6- fraction = table[i]<<16 - (table[i] - table[i+1]) * a * 2 |
\r
89 |___________________________________________________________________________|
\r
91 static Word16 table_isqrt[49] =
\r
93 32767, 31790, 30894, 30070, 29309, 28602, 27945, 27330, 26755, 26214,
\r
94 25705, 25225, 24770, 24339, 23930, 23541, 23170, 22817, 22479, 22155,
\r
95 21845, 21548, 21263, 20988, 20724, 20470, 20225, 19988, 19760, 19539,
\r
96 19326, 19119, 18919, 18725, 18536, 18354, 18176, 18004, 17837, 17674,
\r
97 17515, 17361, 17211, 17064, 16921, 16782, 16646, 16514, 16384
\r
101 Word32 * frac, /* (i/o) Q31: normalized value (1.0 < frac <= 0.5) */
\r
102 Word16 * exp /* (i/o) : exponent (value = frac x 2^exponent) */
\r
107 if (*frac <= (Word32) 0)
\r
110 *frac = 0x7fffffffL;
\r
114 if((*exp & 1) == 1) /*If exponant odd -> shift right */
\r
115 *frac = (*frac) >> 1;
\r
117 *exp = negate((*exp - 1) >> 1);
\r
119 *frac = (*frac >> 9);
\r
120 i = extract_h(*frac); /* Extract b25-b31 */
\r
121 *frac = (*frac >> 1);
\r
122 a = (Word16)(*frac); /* Extract b10-b24 */
\r
123 a = (Word16) (a & (Word16) 0x7fff);
\r
125 *frac = L_deposit_h(table_isqrt[i]); /* table[i] << 16 */
\r
126 tmp = vo_sub(table_isqrt[i], table_isqrt[i + 1]); /* table[i] - table[i+1]) */
\r
127 *frac = vo_L_msu(*frac, tmp, a); /* frac -= tmp*a*2 */
\r
132 /*___________________________________________________________________________
\r
134 | Function Name : Pow2() |
\r
136 | L_x = pow(2.0, exponant.fraction) (exponant = interger part) |
\r
137 | = pow(2.0, 0.fraction) << exponant |
\r
138 |---------------------------------------------------------------------------|
\r
141 | The function Pow2(L_x) is approximated by a table and linear |
\r
144 | 1- i = bit10-b15 of fraction, 0 <= i <= 31 |
\r
145 | 2- a = bit0-b9 of fraction |
\r
146 | 3- L_x = table[i]<<16 - (table[i] - table[i+1]) * a * 2 |
\r
147 | 4- L_x = L_x >> (30-exponant) (with rounding) |
\r
148 |___________________________________________________________________________|
\r
150 static Word16 table_pow2[33] =
\r
152 16384, 16743, 17109, 17484, 17867, 18258, 18658, 19066, 19484, 19911,
\r
153 20347, 20792, 21247, 21713, 22188, 22674, 23170, 23678, 24196, 24726,
\r
154 25268, 25821, 26386, 26964, 27554, 28158, 28774, 29405, 30048, 30706,
\r
155 31379, 32066, 32767
\r
158 Word32 Pow2( /* (o) Q0 : result (range: 0<=val<=0x7fffffff) */
\r
159 Word16 exponant, /* (i) Q0 : Integer part. (range: 0<=val<=30) */
\r
160 Word16 fraction /* (i) Q15 : Fractionnal part. (range: 0.0<=val<1.0) */
\r
163 Word16 exp, i, a, tmp;
\r
166 L_x = vo_L_mult(fraction, 32); /* L_x = fraction<<6 */
\r
167 i = extract_h(L_x); /* Extract b10-b16 of fraction */
\r
169 a = (Word16)(L_x); /* Extract b0-b9 of fraction */
\r
170 a = (Word16) (a & (Word16) 0x7fff);
\r
172 L_x = L_deposit_h(table_pow2[i]); /* table[i] << 16 */
\r
173 tmp = vo_sub(table_pow2[i], table_pow2[i + 1]); /* table[i] - table[i+1] */
\r
174 L_x -= (tmp * a)<<1; /* L_x -= tmp*a*2 */
\r
176 exp = vo_sub(30, exponant);
\r
177 L_x = vo_L_shr_r(L_x, exp);
\r
182 /*___________________________________________________________________________
\r
184 | Function Name : Dot_product12() |
\r
186 | Compute scalar product of <x[],y[]> using accumulator. |
\r
188 | The result is normalized (in Q31) with exponent (0..30). |
\r
189 |---------------------------------------------------------------------------|
\r
192 | dot_product = sum(x[i]*y[i]) i=0..N-1 |
\r
193 |___________________________________________________________________________|
\r
196 Word32 Dot_product12( /* (o) Q31: normalized result (1 < val <= -1) */
\r
197 Word16 x[], /* (i) 12bits: x vector */
\r
198 Word16 y[], /* (i) 12bits: y vector */
\r
199 Word16 lg, /* (i) : vector length */
\r
200 Word16 * exp /* (o) : exponent of result (0..+30) */
\r
206 for (i = 0; i < lg; i++)
\r
208 L_sum += x[i] * y[i];
\r
210 L_sum = (L_sum << 1) + 1;
\r
211 /* Normalize acc in Q31 */
\r
212 sft = norm_l(L_sum);
\r
213 L_sum = L_sum << sft;
\r
214 *exp = 30 - sft; /* exponent = 0..30 */
\r