oper_32b.c revision e2e838afcf03e603a41a0455846eaf9614537c16 
1/*
2 ** Copyright 2003-2010, VisualOn, Inc.
3 **
4 ** Licensed under the Apache License, Version 2.0 (the "License");
5 ** you may not use this file except in compliance with the License.
6 ** You may obtain a copy of the License at
7 **
8 **     http://www.apache.org/licenses/LICENSE-2.0
9 **
10 ** Unless required by applicable law or agreed to in writing, software
11 ** distributed under the License is distributed on an "AS IS" BASIS,
12 ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 ** See the License for the specific language governing permissions and
14 ** limitations under the License.
15 */
16/*******************************************************************************
17	File:		oper_32b.c
18
19	Content:	  This file contains operations in double precision.
20
21*******************************************************************************/
22
23#include "typedef.h"
24#include "basic_op.h"
25#include "oper_32b.h"
26
27/*****************************************************************************
28 *                                                                           *
29 *  Function L_Extract()                                                     *
30 *                                                                           *
31 *  Extract from a 32 bit integer two 16 bit DPF.                            *
32 *                                                                           *
33 *  Arguments:                                                               *
34 *                                                                           *
35 *   L_32      : 32 bit integer.                                             *
36 *               0x8000 0000 <= L_32 <= 0x7fff ffff.                         *
37 *   hi        : b16 to b31 of L_32                                          *
38 *   lo        : (L_32 - hi<<16)>>1                                          *
39 *****************************************************************************
40*/
41
42void L_Extract (Word32 L_32, Word16 *hi, Word16 *lo)
43{
44    *hi = extract_h (L_32);
45    *lo = extract_l (L_msu (L_shr (L_32, 1), *hi, 16384));
46    return;
47}
48
49/*****************************************************************************
50 *                                                                           *
51 *  Function L_Comp()                                                        *
52 *                                                                           *
53 *  Compose from two 16 bit DPF a 32 bit integer.                            *
54 *                                                                           *
55 *     L_32 = hi<<16 + lo<<1                                                 *
56 *                                                                           *
57 *  Arguments:                                                               *
58 *                                                                           *
59 *   hi        msb                                                           *
60 *   lo        lsf (with sign)                                               *
61 *                                                                           *
62 *   Return Value :                                                          *
63 *                                                                           *
64 *             32 bit long signed integer (Word32) whose value falls in the  *
65 *             range : 0x8000 0000 <= L_32 <= 0x7fff fff0.                   *
66 *                                                                           *
67 *****************************************************************************
68*/
69
70Word32 L_Comp (Word16 hi, Word16 lo)
71{
72    Word32 L_32;
73
74    L_32 = L_deposit_h (hi);
75    return (L_mac (L_32, lo, 1));       /* = hi<<16 + lo<<1 */
76}
77
78/*****************************************************************************
79 * Function Mpy_32()                                                         *
80 *                                                                           *
81 *   Multiply two 32 bit integers (DPF). The result is divided by 2**31      *
82 *                                                                           *
83 *   L_32 = (hi1*hi2)<<1 + ( (hi1*lo2)>>15 + (lo1*hi2)>>15 )<<1              *
84 *                                                                           *
85 *   This operation can also be viewed as the multiplication of two Q31      *
86 *   number and the result is also in Q31.                                   *
87 *                                                                           *
88 * Arguments:                                                                *
89 *                                                                           *
90 *  hi1         hi part of first number                                      *
91 *  lo1         lo part of first number                                      *
92 *  hi2         hi part of second number                                     *
93 *  lo2         lo part of second number                                     *
94 *                                                                           *
95 *****************************************************************************
96*/
97
98Word32 Mpy_32 (Word16 hi1, Word16 lo1, Word16 hi2, Word16 lo2)
99{
100    Word32 L_32;
101
102    L_32 = L_mult (hi1, hi2);
103    L_32 = L_mac (L_32, mult (hi1, lo2), 1);
104    L_32 = L_mac (L_32, mult (lo1, hi2), 1);
105
106    return (L_32);
107}
108
109/*****************************************************************************
110 * Function Mpy_32_16()                                                      *
111 *                                                                           *
112 *   Multiply a 16 bit integer by a 32 bit (DPF). The result is divided      *
113 *   by 2**15                                                                *
114 *                                                                           *
115 *                                                                           *
116 *   L_32 = (hi1*lo2)<<1 + ((lo1*lo2)>>15)<<1                                *
117 *                                                                           *
118 * Arguments:                                                                *
119 *                                                                           *
120 *  hi          hi part of 32 bit number.                                    *
121 *  lo          lo part of 32 bit number.                                    *
122 *  n           16 bit number.                                               *
123 *                                                                           *
124 *****************************************************************************
125*/
126
127Word32 Mpy_32_16 (Word16 hi, Word16 lo, Word16 n)
128{
129    Word32 L_32;
130
131    L_32 = L_mult (hi, n);
132    L_32 = L_mac (L_32, mult (lo, n), 1);
133
134    return (L_32);
135}
136
137/*****************************************************************************
138 *                                                                           *
139 *   Function Name : Div_32                                                  *
140 *                                                                           *
141 *   Purpose :                                                               *
142 *             Fractional integer division of two 32 bit numbers.            *
143 *             L_num / L_denom.                                              *
144 *             L_num and L_denom must be positive and L_num < L_denom.       *
145 *             L_denom = denom_hi<<16 + denom_lo<<1                          *
146 *             denom_hi is a normalize number.                               *
147 *                                                                           *
148 *   Inputs :                                                                *
149 *                                                                           *
150 *    L_num                                                                  *
151 *             32 bit long signed integer (Word32) whose value falls in the  *
152 *             range : 0x0000 0000 < L_num < L_denom                         *
153 *                                                                           *
154 *    L_denom = denom_hi<<16 + denom_lo<<1      (DPF)                        *
155 *                                                                           *
156 *       denom_hi                                                            *
157 *             16 bit positive normalized integer whose value falls in the   *
158 *             range : 0x4000 < hi < 0x7fff                                  *
159 *       denom_lo                                                            *
160 *             16 bit positive integer whose value falls in the              *
161 *             range : 0 < lo < 0x7fff                                       *
162 *                                                                           *
163 *   Return Value :                                                          *
164 *                                                                           *
165 *    L_div                                                                  *
166 *             32 bit long signed integer (Word32) whose value falls in the  *
167 *             range : 0x0000 0000 <= L_div <= 0x7fff ffff.                  *
168 *                                                                           *
169 *  Algorithm:                                                               *
170 *                                                                           *
171 *  - find = 1/L_denom.                                                      *
172 *      First approximation: approx = 1 / denom_hi                           *
173 *      1/L_denom = approx * (2.0 - L_denom * approx )                       *
174 *                                                                           *
175 *  -  result = L_num * (1/L_denom)                                          *
176 *****************************************************************************
177*/
178
179Word32 Div_32 (Word32 L_num, Word32 denom)
180{
181    Word16 approx;
182    Word32 L_32;
183    /* First approximation: 1 / L_denom = 1/denom_hi */
184
185    approx = div_s ((Word16) 0x3fff, denom >> 16);
186
187    /* 1/L_denom = approx * (2.0 - L_denom * approx) */
188
189    L_32 = L_mpy_ls (denom, approx);
190
191    L_32 = L_sub ((Word32) 0x7fffffffL, L_32);
192
193	L_32 = L_mpy_ls (L_32, approx);
194    /* L_num * (1/L_denom) */
195
196	L_32 = MULHIGH(L_32, L_num);
197    L_32 = L_shl (L_32, 3);
198
199    return (L_32);
200}
201
202/*!
203
204  \brief  calculates the log dualis times 4 of argument
205          iLog4(x) = (Word32)(4 * log(value)/log(2.0))
206
207  \return ilog4 value
208
209*/
210Word16 iLog4(Word32 value)
211{
212  Word16 iLog4;
213
214  if(value != 0){
215    Word32 tmp;
216    Word16 tmp16;
217    iLog4 = norm_l(value);
218    tmp = (value << iLog4);
219    tmp16 = round16(tmp);
220    tmp = L_mult(tmp16, tmp16);
221    tmp16 = round16(tmp);
222    tmp = L_mult(tmp16, tmp16);
223    tmp16 = round16(tmp);
224
225    iLog4 = (-(iLog4 << 2) - norm_s(tmp16)) - 1;
226  }
227  else {
228    iLog4 = -128; /* -(INT_BITS*4); */
229  }
230
231  return iLog4;
232}
233
234#define step(shift) \
235    if ((0x40000000l >> shift) + root <= value)       \
236    {                                                 \
237        value -= (0x40000000l >> shift) + root;       \
238        root = (root >> 1) | (0x40000000l >> shift);  \
239    } else {                                          \
240        root = root >> 1;                             \
241    }
242
243Word32 rsqrt(Word32 value,     /*!< Operand to square root (0.0 ... 1) */
244             Word32 accuracy)  /*!< Number of valid bits that will be calculated */
245{
246    Word32 root = 0;
247	Word32 scale;
248
249	if(value < 0)
250		return 0;
251
252	scale = norm_l(value);
253	if(scale & 1) scale--;
254
255	value <<= scale;
256
257	step( 0); step( 2); step( 4); step( 6);
258    step( 8); step(10); step(12); step(14);
259    step(16); step(18); step(20); step(22);
260    step(24); step(26); step(28); step(30);
261
262    scale >>= 1;
263	if (root < value)
264        ++root;
265
266	root >>= scale;
267    return root* 46334;
268}
269
270static const Word32 pow2Table[POW2_TABLE_SIZE] = {
2710x7fffffff, 0x7fa765ad, 0x7f4f08ae, 0x7ef6e8da,
2720x7e9f0606, 0x7e476009, 0x7deff6b6, 0x7d98c9e6,
2730x7d41d96e, 0x7ceb2523, 0x7c94acde, 0x7c3e7073,
2740x7be86fb9, 0x7b92aa88, 0x7b3d20b6, 0x7ae7d21a,
2750x7a92be8b, 0x7a3de5df, 0x79e947ef, 0x7994e492,
2760x7940bb9e, 0x78ecccec, 0x78991854, 0x78459dac,
2770x77f25cce, 0x779f5591, 0x774c87cc, 0x76f9f359,
2780x76a7980f, 0x765575c8, 0x76038c5b, 0x75b1dba2,
2790x75606374, 0x750f23ab, 0x74be1c20, 0x746d4cac,
2800x741cb528, 0x73cc556d, 0x737c2d55, 0x732c3cba,
2810x72dc8374, 0x728d015d, 0x723db650, 0x71eea226,
2820x719fc4b9, 0x71511de4, 0x7102ad80, 0x70b47368,
2830x70666f76, 0x7018a185, 0x6fcb096f, 0x6f7da710,
2840x6f307a41, 0x6ee382de, 0x6e96c0c3, 0x6e4a33c9,
2850x6dfddbcc, 0x6db1b8a8, 0x6d65ca38, 0x6d1a1057,
2860x6cce8ae1, 0x6c8339b2, 0x6c381ca6, 0x6bed3398,
2870x6ba27e66, 0x6b57fce9, 0x6b0daeff, 0x6ac39485,
2880x6a79ad56, 0x6a2ff94f, 0x69e6784d, 0x699d2a2c,
2890x69540ec9, 0x690b2601, 0x68c26fb1, 0x6879ebb6,
2900x683199ed, 0x67e97a34, 0x67a18c68, 0x6759d065,
2910x6712460b, 0x66caed35, 0x6683c5c3, 0x663ccf92,
2920x65f60a80, 0x65af766a, 0x6569132f, 0x6522e0ad,
2930x64dcdec3, 0x64970d4f, 0x64516c2e, 0x640bfb41,
2940x63c6ba64, 0x6381a978, 0x633cc85b, 0x62f816eb,
2950x62b39509, 0x626f4292, 0x622b1f66, 0x61e72b65,
2960x61a3666d, 0x615fd05f, 0x611c6919, 0x60d9307b,
2970x60962665, 0x60534ab7, 0x60109d51, 0x5fce1e12,
2980x5f8bccdb, 0x5f49a98c, 0x5f07b405, 0x5ec5ec26,
2990x5e8451d0, 0x5e42e4e3, 0x5e01a540, 0x5dc092c7,
3000x5d7fad59, 0x5d3ef4d7, 0x5cfe6923, 0x5cbe0a1c,
3010x5c7dd7a4, 0x5c3dd19c, 0x5bfdf7e5, 0x5bbe4a61,
3020x5b7ec8f2, 0x5b3f7377, 0x5b0049d4, 0x5ac14bea,
3030x5a82799a, 0x5a43d2c6, 0x5a055751, 0x59c7071c,
3040x5988e209, 0x594ae7fb, 0x590d18d3, 0x58cf7474,
3050x5891fac1, 0x5854ab9b, 0x581786e6, 0x57da8c83,
3060x579dbc57, 0x57611642, 0x57249a29, 0x56e847ef,
3070x56ac1f75, 0x567020a0, 0x56344b52, 0x55f89f70,
3080x55bd1cdb, 0x5581c378, 0x55469329, 0x550b8bd4,
3090x54d0ad5b, 0x5495f7a1, 0x545b6a8b, 0x542105fd,
3100x53e6c9db, 0x53acb607, 0x5372ca68, 0x533906e0,
3110x52ff6b55, 0x52c5f7aa, 0x528cabc3, 0x52538786,
3120x521a8ad7, 0x51e1b59a, 0x51a907b4, 0x5170810b,
3130x51382182, 0x50ffe8fe, 0x50c7d765, 0x508fec9c,
3140x50582888, 0x50208b0e, 0x4fe91413, 0x4fb1c37c,
3150x4f7a9930, 0x4f439514, 0x4f0cb70c, 0x4ed5ff00,
3160x4e9f6cd4, 0x4e69006e, 0x4e32b9b4, 0x4dfc988c,
3170x4dc69cdd, 0x4d90c68b, 0x4d5b157e, 0x4d25899c,
3180x4cf022ca, 0x4cbae0ef, 0x4c85c3f1, 0x4c50cbb8,
3190x4c1bf829, 0x4be7492b, 0x4bb2bea5, 0x4b7e587d,
3200x4b4a169c, 0x4b15f8e6, 0x4ae1ff43, 0x4aae299b,
3210x4a7a77d5, 0x4a46e9d6, 0x4a137f88, 0x49e038d0,
3220x49ad1598, 0x497a15c4, 0x4947393f, 0x49147fee,
3230x48e1e9ba, 0x48af768a, 0x487d2646, 0x484af8d6,
3240x4818ee22, 0x47e70611, 0x47b5408c, 0x47839d7b,
3250x47521cc6, 0x4720be55, 0x46ef8210, 0x46be67e0,
3260x468d6fae, 0x465c9961, 0x462be4e2, 0x45fb521a,
3270x45cae0f2, 0x459a9152, 0x456a6323, 0x453a564d,
3280x450a6abb, 0x44daa054, 0x44aaf702, 0x447b6ead,
3290x444c0740, 0x441cc0a3, 0x43ed9ac0, 0x43be9580,
3300x438fb0cb, 0x4360ec8d, 0x433248ae, 0x4303c517,
3310x42d561b4, 0x42a71e6c, 0x4278fb2b, 0x424af7da,
3320x421d1462, 0x41ef50ae, 0x41c1aca8, 0x41942839,
3330x4166c34c, 0x41397dcc, 0x410c57a2, 0x40df50b8,
3340x40b268fa, 0x4085a051, 0x4058f6a8, 0x402c6be9
335};
336
337/*!
338
339  \brief calculates 2 ^ (x/y) for x<=0, y > 0, x <= 32768 * y
340
341  avoids integer division
342
343  \return
344*/
345Word32 pow2_xy(Word32 x, Word32 y)
346{
347  Word32 iPart;
348  Word32 fPart;
349  Word32 res;
350  Word32 tmp, tmp2;
351  Word32 shift, shift2;
352
353  tmp2 = -x;
354  iPart = tmp2 / y;
355  fPart = tmp2 - iPart*y;
356  iPart = min(iPart,INT_BITS-1);
357
358  res = pow2Table[(POW2_TABLE_SIZE*fPart)/y] >> iPart;
359
360  return(res);
361}