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 #define UNUSED(x) (void)(x) 28 29 /***************************************************************************** 30 * * 31 * Function L_Extract() * 32 * * 33 * Extract from a 32 bit integer two 16 bit DPF. * 34 * * 35 * Arguments: * 36 * * 37 * L_32 : 32 bit integer. * 38 * 0x8000 0000 <= L_32 <= 0x7fff ffff. * 39 * hi : b16 to b31 of L_32 * 40 * lo : (L_32 - hi<<16)>>1 * 41 ***************************************************************************** 42 */ 43 44 void L_Extract (Word32 L_32, Word16 *hi, Word16 *lo) 45 { 46 *hi = extract_h (L_32); 47 *lo = extract_l (L_msu (L_shr (L_32, 1), *hi, 16384)); 48 return; 49 } 50 51 /***************************************************************************** 52 * * 53 * Function L_Comp() * 54 * * 55 * Compose from two 16 bit DPF a 32 bit integer. * 56 * * 57 * L_32 = hi<<16 + lo<<1 * 58 * * 59 * Arguments: * 60 * * 61 * hi msb * 62 * lo lsf (with sign) * 63 * * 64 * Return Value : * 65 * * 66 * 32 bit long signed integer (Word32) whose value falls in the * 67 * range : 0x8000 0000 <= L_32 <= 0x7fff fff0. * 68 * * 69 ***************************************************************************** 70 */ 71 72 Word32 L_Comp (Word16 hi, Word16 lo) 73 { 74 Word32 L_32; 75 76 L_32 = L_deposit_h (hi); 77 return (L_mac (L_32, lo, 1)); /* = hi<<16 + lo<<1 */ 78 } 79 80 /***************************************************************************** 81 * Function Mpy_32() * 82 * * 83 * Multiply two 32 bit integers (DPF). The result is divided by 2**31 * 84 * * 85 * L_32 = (hi1*hi2)<<1 + ( (hi1*lo2)>>15 + (lo1*hi2)>>15 )<<1 * 86 * * 87 * This operation can also be viewed as the multiplication of two Q31 * 88 * number and the result is also in Q31. * 89 * * 90 * Arguments: * 91 * * 92 * hi1 hi part of first number * 93 * lo1 lo part of first number * 94 * hi2 hi part of second number * 95 * lo2 lo part of second number * 96 * * 97 ***************************************************************************** 98 */ 99 100 Word32 Mpy_32 (Word16 hi1, Word16 lo1, Word16 hi2, Word16 lo2) 101 { 102 Word32 L_32; 103 104 L_32 = L_mult (hi1, hi2); 105 L_32 = L_mac (L_32, mult (hi1, lo2), 1); 106 L_32 = L_mac (L_32, mult (lo1, hi2), 1); 107 108 return (L_32); 109 } 110 111 /***************************************************************************** 112 * Function Mpy_32_16() * 113 * * 114 * Multiply a 16 bit integer by a 32 bit (DPF). The result is divided * 115 * by 2**15 * 116 * * 117 * * 118 * L_32 = (hi1*lo2)<<1 + ((lo1*lo2)>>15)<<1 * 119 * * 120 * Arguments: * 121 * * 122 * hi hi part of 32 bit number. * 123 * lo lo part of 32 bit number. * 124 * n 16 bit number. * 125 * * 126 ***************************************************************************** 127 */ 128 129 Word32 Mpy_32_16 (Word16 hi, Word16 lo, Word16 n) 130 { 131 Word32 L_32; 132 133 L_32 = L_mult (hi, n); 134 L_32 = L_mac (L_32, mult (lo, n), 1); 135 136 return (L_32); 137 } 138 139 /***************************************************************************** 140 * * 141 * Function Name : Div_32 * 142 * * 143 * Purpose : * 144 * Fractional integer division of two 32 bit numbers. * 145 * L_num / L_denom. * 146 * L_num and L_denom must be positive and L_num < L_denom. * 147 * L_denom = denom_hi<<16 + denom_lo<<1 * 148 * denom_hi is a normalize number. * 149 * * 150 * Inputs : * 151 * * 152 * L_num * 153 * 32 bit long signed integer (Word32) whose value falls in the * 154 * range : 0x0000 0000 < L_num < L_denom * 155 * * 156 * L_denom = denom_hi<<16 + denom_lo<<1 (DPF) * 157 * * 158 * denom_hi * 159 * 16 bit positive normalized integer whose value falls in the * 160 * range : 0x4000 < hi < 0x7fff * 161 * denom_lo * 162 * 16 bit positive integer whose value falls in the * 163 * range : 0 < lo < 0x7fff * 164 * * 165 * Return Value : * 166 * * 167 * L_div * 168 * 32 bit long signed integer (Word32) whose value falls in the * 169 * range : 0x0000 0000 <= L_div <= 0x7fff ffff. * 170 * * 171 * Algorithm: * 172 * * 173 * - find = 1/L_denom. * 174 * First approximation: approx = 1 / denom_hi * 175 * 1/L_denom = approx * (2.0 - L_denom * approx ) * 176 * * 177 * - result = L_num * (1/L_denom) * 178 ***************************************************************************** 179 */ 180 181 Word32 Div_32 (Word32 L_num, Word32 denom) 182 { 183 Word16 approx; 184 Word32 L_32; 185 /* First approximation: 1 / L_denom = 1/denom_hi */ 186 187 approx = div_s ((Word16) 0x3fff, denom >> 16); 188 189 /* 1/L_denom = approx * (2.0 - L_denom * approx) */ 190 191 L_32 = L_mpy_ls (denom, approx); 192 193 L_32 = L_sub ((Word32) 0x7fffffffL, L_32); 194 195 L_32 = L_mpy_ls (L_32, approx); 196 /* L_num * (1/L_denom) */ 197 198 L_32 = MULHIGH(L_32, L_num); 199 L_32 = L_shl (L_32, 3); 200 201 return (L_32); 202 } 203 204 /*! 205 206 \brief calculates the log dualis times 4 of argument 207 iLog4(x) = (Word32)(4 * log(value)/log(2.0)) 208 209 \return ilog4 value 210 211 */ 212 Word16 iLog4(Word32 value) 213 { 214 Word16 iLog4; 215 216 if(value != 0){ 217 Word32 tmp; 218 Word16 tmp16; 219 iLog4 = norm_l(value); 220 tmp = (value << iLog4); 221 tmp16 = round16(tmp); 222 tmp = L_mult(tmp16, tmp16); 223 tmp16 = round16(tmp); 224 tmp = L_mult(tmp16, tmp16); 225 tmp16 = round16(tmp); 226 227 iLog4 = (-(iLog4 << 2) - norm_s(tmp16)) - 1; 228 } 229 else { 230 iLog4 = -128; /* -(INT_BITS*4); */ 231 } 232 233 return iLog4; 234 } 235 236 #define step(shift) \ 237 if ((0x40000000l >> shift) + root <= value) \ 238 { \ 239 value -= (0x40000000l >> shift) + root; \ 240 root = (root >> 1) | (0x40000000l >> shift); \ 241 } else { \ 242 root = root >> 1; \ 243 } 244 245 Word32 rsqrt(Word32 value, /*!< Operand to square root (0.0 ... 1) */ 246 Word32 accuracy) /*!< Number of valid bits that will be calculated */ 247 { 248 UNUSED(accuracy); 249 250 Word32 root = 0; 251 Word32 scale; 252 253 if(value < 0) 254 return 0; 255 256 scale = norm_l(value); 257 if(scale & 1) scale--; 258 259 value <<= scale; 260 261 step( 0); step( 2); step( 4); step( 6); 262 step( 8); step(10); step(12); step(14); 263 step(16); step(18); step(20); step(22); 264 step(24); step(26); step(28); step(30); 265 266 scale >>= 1; 267 if (root < value) 268 ++root; 269 270 root >>= scale; 271 return root* 46334; 272 } 273 274 static const Word32 pow2Table[POW2_TABLE_SIZE] = { 275 0x7fffffff, 0x7fa765ad, 0x7f4f08ae, 0x7ef6e8da, 276 0x7e9f0606, 0x7e476009, 0x7deff6b6, 0x7d98c9e6, 277 0x7d41d96e, 0x7ceb2523, 0x7c94acde, 0x7c3e7073, 278 0x7be86fb9, 0x7b92aa88, 0x7b3d20b6, 0x7ae7d21a, 279 0x7a92be8b, 0x7a3de5df, 0x79e947ef, 0x7994e492, 280 0x7940bb9e, 0x78ecccec, 0x78991854, 0x78459dac, 281 0x77f25cce, 0x779f5591, 0x774c87cc, 0x76f9f359, 282 0x76a7980f, 0x765575c8, 0x76038c5b, 0x75b1dba2, 283 0x75606374, 0x750f23ab, 0x74be1c20, 0x746d4cac, 284 0x741cb528, 0x73cc556d, 0x737c2d55, 0x732c3cba, 285 0x72dc8374, 0x728d015d, 0x723db650, 0x71eea226, 286 0x719fc4b9, 0x71511de4, 0x7102ad80, 0x70b47368, 287 0x70666f76, 0x7018a185, 0x6fcb096f, 0x6f7da710, 288 0x6f307a41, 0x6ee382de, 0x6e96c0c3, 0x6e4a33c9, 289 0x6dfddbcc, 0x6db1b8a8, 0x6d65ca38, 0x6d1a1057, 290 0x6cce8ae1, 0x6c8339b2, 0x6c381ca6, 0x6bed3398, 291 0x6ba27e66, 0x6b57fce9, 0x6b0daeff, 0x6ac39485, 292 0x6a79ad56, 0x6a2ff94f, 0x69e6784d, 0x699d2a2c, 293 0x69540ec9, 0x690b2601, 0x68c26fb1, 0x6879ebb6, 294 0x683199ed, 0x67e97a34, 0x67a18c68, 0x6759d065, 295 0x6712460b, 0x66caed35, 0x6683c5c3, 0x663ccf92, 296 0x65f60a80, 0x65af766a, 0x6569132f, 0x6522e0ad, 297 0x64dcdec3, 0x64970d4f, 0x64516c2e, 0x640bfb41, 298 0x63c6ba64, 0x6381a978, 0x633cc85b, 0x62f816eb, 299 0x62b39509, 0x626f4292, 0x622b1f66, 0x61e72b65, 300 0x61a3666d, 0x615fd05f, 0x611c6919, 0x60d9307b, 301 0x60962665, 0x60534ab7, 0x60109d51, 0x5fce1e12, 302 0x5f8bccdb, 0x5f49a98c, 0x5f07b405, 0x5ec5ec26, 303 0x5e8451d0, 0x5e42e4e3, 0x5e01a540, 0x5dc092c7, 304 0x5d7fad59, 0x5d3ef4d7, 0x5cfe6923, 0x5cbe0a1c, 305 0x5c7dd7a4, 0x5c3dd19c, 0x5bfdf7e5, 0x5bbe4a61, 306 0x5b7ec8f2, 0x5b3f7377, 0x5b0049d4, 0x5ac14bea, 307 0x5a82799a, 0x5a43d2c6, 0x5a055751, 0x59c7071c, 308 0x5988e209, 0x594ae7fb, 0x590d18d3, 0x58cf7474, 309 0x5891fac1, 0x5854ab9b, 0x581786e6, 0x57da8c83, 310 0x579dbc57, 0x57611642, 0x57249a29, 0x56e847ef, 311 0x56ac1f75, 0x567020a0, 0x56344b52, 0x55f89f70, 312 0x55bd1cdb, 0x5581c378, 0x55469329, 0x550b8bd4, 313 0x54d0ad5b, 0x5495f7a1, 0x545b6a8b, 0x542105fd, 314 0x53e6c9db, 0x53acb607, 0x5372ca68, 0x533906e0, 315 0x52ff6b55, 0x52c5f7aa, 0x528cabc3, 0x52538786, 316 0x521a8ad7, 0x51e1b59a, 0x51a907b4, 0x5170810b, 317 0x51382182, 0x50ffe8fe, 0x50c7d765, 0x508fec9c, 318 0x50582888, 0x50208b0e, 0x4fe91413, 0x4fb1c37c, 319 0x4f7a9930, 0x4f439514, 0x4f0cb70c, 0x4ed5ff00, 320 0x4e9f6cd4, 0x4e69006e, 0x4e32b9b4, 0x4dfc988c, 321 0x4dc69cdd, 0x4d90c68b, 0x4d5b157e, 0x4d25899c, 322 0x4cf022ca, 0x4cbae0ef, 0x4c85c3f1, 0x4c50cbb8, 323 0x4c1bf829, 0x4be7492b, 0x4bb2bea5, 0x4b7e587d, 324 0x4b4a169c, 0x4b15f8e6, 0x4ae1ff43, 0x4aae299b, 325 0x4a7a77d5, 0x4a46e9d6, 0x4a137f88, 0x49e038d0, 326 0x49ad1598, 0x497a15c4, 0x4947393f, 0x49147fee, 327 0x48e1e9ba, 0x48af768a, 0x487d2646, 0x484af8d6, 328 0x4818ee22, 0x47e70611, 0x47b5408c, 0x47839d7b, 329 0x47521cc6, 0x4720be55, 0x46ef8210, 0x46be67e0, 330 0x468d6fae, 0x465c9961, 0x462be4e2, 0x45fb521a, 331 0x45cae0f2, 0x459a9152, 0x456a6323, 0x453a564d, 332 0x450a6abb, 0x44daa054, 0x44aaf702, 0x447b6ead, 333 0x444c0740, 0x441cc0a3, 0x43ed9ac0, 0x43be9580, 334 0x438fb0cb, 0x4360ec8d, 0x433248ae, 0x4303c517, 335 0x42d561b4, 0x42a71e6c, 0x4278fb2b, 0x424af7da, 336 0x421d1462, 0x41ef50ae, 0x41c1aca8, 0x41942839, 337 0x4166c34c, 0x41397dcc, 0x410c57a2, 0x40df50b8, 338 0x40b268fa, 0x4085a051, 0x4058f6a8, 0x402c6be9 339 }; 340 341 /*! 342 343 \brief calculates 2 ^ (x/y) for x<=0, y > 0, x <= 32768 * y 344 345 avoids integer division 346 347 \return 348 */ 349 Word32 pow2_xy(Word32 x, Word32 y) 350 { 351 UWord32 iPart; 352 UWord32 fPart; 353 Word32 res; 354 Word32 tmp, tmp2; 355 Word32 shift, shift2; 356 357 tmp2 = -x; 358 iPart = tmp2 / y; 359 fPart = tmp2 - iPart*y; 360 iPart = min(iPart,INT_BITS-1); 361 362 res = pow2Table[(POW2_TABLE_SIZE*fPart)/y] >> iPart; 363 364 return(res); 365 } 366
