release/src-rt-6.x.4708/shared/qmath.c

   1 /*
   2  * qmath functions used in arithmatic and DSP operations where
   3  * fractional operations, saturation support is needed.
   4  *
   5  * Copyright (C) 2012, Broadcom Corporation
   6  * All Rights Reserved.
   7  *
   8  * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Broadcom Corporation;
   9  * the contents of this file may not be disclosed to third parties, copied
  10  * or duplicated in any form, in whole or in part, without the prior
  11  * written permission of Broadcom Corporation.
  12  *
  13  * $Id: qmath.c 300516 2011-12-04 17:39:44Z $
  14  */
  15
  16 #include <bcm_cfg.h>
  17 #include "qmath.h"
  18
  19 /*
  20 Description: This function saturate input 32 bit number into a 16 bit number.
  21 If input number is greater than 0x7fff then output is saturated to 0x7fff.
  22 else if input number is less than 0xffff8000 then output is saturated to 0xffff8000
  23 else output is same as input.
  24 */
  25 int16
  26 qm_sat32(int32 op)
  27 {
  28         int16 result;
  29         if (op > (int32)0x7fff) {
  30                 result = 0x7fff;
  31         }
  32         else if (op < (int32)0xffff8000) {
  33                 result = (int16)(0x8000);
  34         }
  35         else {
  36                 result = (int16)op;
  37         }
  38         return result;
  39 }
  40
  41 /*
  42 Description: This function multiply two input 16 bit numbers and return the 32 bit result.
  43 This multiplication is similar to compiler multiplication. This operation is defined if
  44 16 bit multiplication on the processor platform is cheaper than 32 bit multiplication (as
  45 the most of qmath functions can be replaced with processor intrinsic instructions).
  46 */
  47 int32
  48 qm_mul321616(int16 op1, int16 op2)
  49 {
  50         return ((int32)(op1)*(int32)(op2));
  51 }
  52
  53 /*
  54 Description: This function make 16 bit multiplication and return the result in 16 bits.
  55 To fit the result into 16 bits the 32 bit multiplication result is right
  56 shifted by 16 bits.
  57 */
  58 int16
  59 qm_mul16(int16 op1, int16 op2)
  60 {
  61         int32 result;
  62         result = ((int32)(op1)*(int32)(op2));
  63         return ((int16)(result >> 16));
  64 }
  65
  66 /*
  67 Description: This function multiply two 16 bit numbers and return the result in 32 bits.
  68 This function remove the extra sign bit created by the multiplication by leftshifting the
  69 32 bit multiplication result by 1 bit before returning the result. So the output is
  70 twice that of compiler multiplication. (i.e. qm_muls321616(2,3)=12).
  71 When both input 16 bit numbers are 0x8000, then the result is saturated to 0x7fffffff.
  72 */
  73 int32
  74 qm_muls321616(int16 op1, int16 op2)
  75 {
  76         int32 result;
  77         if (op1 == (int16)(0x8000) && op2 == (int16)(0x8000)) {
  78                 result = 0x7fffffff;
  79         }
  80         else {
  81                 result = ((int32)(op1)*(int32)(op2));
  82                 result = result << 1;
  83         }
  84         return result;
  85 }
  86
  87 /*
  88 Description: This function make 16 bit unsigned multiplication. To fit the output into
  89 16 bits the 32 bit multiplication result is right shifted by 16 bits.
  90 */
  91 uint16
  92 qm_mulu16(uint16 op1, uint16 op2)
  93 {
  94         return (uint16)(((uint32)op1 * (uint32)op2)>>16);
  95 }
  96
  97 /*
  98 Description: This function make 16 bit multiplication and return the result in 16 bits.
  99 To fit the multiplication result into 16 bits the multiplication result is right shifted by
 100 15 bits. Right shifting 15 bits instead of 16 bits is done to remove the extra sign bit formed
 101 due to the multiplication.
 102 When both the 16bit inputs are 0x8000 then the output is saturated to 0x7fffffff.
 103 */
 104 int16
 105 qm_muls16(int16 op1, int16 op2)
 106 {
 107         int32 result;
 108         if (op1 == (int16)0x8000 && op2 == (int16)0x8000) {
 109                 result = 0x7fffffff;
 110         }
 111         else {
 112                 result = ((int32)(op1)*(int32)(op2));
 113         }
 114         return ((int16)(result >> 15));
 115 }
 116
 117 /*
 118 Description: This function add two 32 bit numbers and return the 32bit result.
 119 If the result overflow 32 bits, the output will be saturated to 32bits.
 120 */
 121 int32
 122 qm_add32(int32 op1, int32 op2)
 123 {
 124         int32 result;
 125         result = op1 + op2;
 126         if (op1 < 0 && op2 < 0 && result > 0) {
 127                 result = 0x80000000;
 128         } else if (op1 > 0 && op2 > 0 && result < 0) {
 129                 result = 0x7fffffff;
 130         }
 131         return result;
 132 }
 133
 134 /*
 135 Description: This function add two 16 bit numbers and return the 16bit result.
 136 If the result overflow 16 bits, the output will be saturated to 16bits.
 137 */
 138 int16
 139 qm_add16(int16 op1, int16 op2)
 140 {
 141         int16 result;
 142         int32 temp = (int32)op1 + (int32)op2;
 143         if (temp > (int32)0x7fff) {
 144                 result = (int16)0x7fff;
 145         } else if (temp < (int32)0xffff8000) {
 146                 result = (int16)0xffff8000;
 147         } else {
 148                 result = (int16) temp;
 149         }
 150         return result;
 151 }
 152
 153 /*
 154 Description: This function make 16 bit subtraction and return the 16bit result.
 155 If the result overflow 16 bits, the output will be saturated to 16bits.
 156 */
 157 int16
 158 qm_sub16(int16 op1, int16 op2)
 159 {
 160         int16 result;
 161         int32 temp = (int32)op1 - (int32)op2;
 162         if (temp > (int32)0x7fff) {
 163                 result = (int16)0x7fff;
 164         } else if (temp < (int32)0xffff8000) {
 165                 result = (int16)0xffff8000;
 166         } else {
 167                 result = (int16) temp;
 168         }
 169         return result;
 170 }
 171
 172 /*
 173 Description: This function make 32 bit subtraction and return the 32bit result.
 174 If the result overflow 32 bits, the output will be saturated to 32bits.
 175 */
 176 int32
 177 qm_sub32(int32 op1, int32 op2)
 178 {
 179         int32 result;
 180         result = op1 - op2;
 181         if (op1 >= 0 && op2 < 0 && result < 0) {
 182                 result = 0x7fffffff;
 183         }
 184         else if (op1 < 0 && op2 > 0 && result > 0)      {
 185                 result = 0x80000000;
 186         }
 187         return result;
 188 }
 189
 190 /*
 191 Description: This function multiply input 16 bit numbers and accumulate the result
 192 into the input 32 bit number and return the 32 bit accumulated result.
 193 If the accumulation result in overflow, then the output will be saturated.
 194 */
 195 int32
 196 qm_mac321616(int32 acc, int16 op1, int16 op2)
 197 {
 198         int32 result;
 199         result = qm_add32(acc, qm_mul321616(op1, op2));
 200         return result;
 201 }
 202
 203 /*
 204 Description: This function make a 32 bit saturated left shift when the specified shift
 205 is +ve. This function will make a 32 bit right shift when the specified shift is -ve.
 206 This function return the result after shifting operation.
 207 */
 208 int32
 209 qm_shl32(int32 op, int shift)
 210 {
 211         int i;
 212         int32 result;
 213         result = op;
 214         if (shift > 31) shift = 31;
 215         else if (shift < -31) shift = -31;
 216         if (shift >= 0) {
 217                 for (i = 0; i < shift; i++) {
 218                         result = qm_add32(result, result);
 219                 }
 220         }
 221         else {
 222                 result = result >> (-shift);
 223         }
 224         return result;
 225 }
 226
 227 /*
 228 Description: This function make a 32 bit right shift when shift is +ve.
 229 This function make a 32 bit saturated left shift when shift is -ve. This function
 230 return the result of the shift operation.
 231 */
 232 int32
 233 qm_shr32(int32 op, int shift)
 234 {
 235         return qm_shl32(op, -shift);
 236 }
 237
 238 /*
 239 Description: This function make a 16 bit saturated left shift when the specified shift
 240 is +ve. This function will make a 16 bit right shift when the specified shift is -ve.
 241 This function return the result after shifting operation.
 242 */
 243 int16
 244 qm_shl16(int16 op, int shift)
 245 {
 246         int i;
 247         int16 result;
 248         result = op;
 249         if (shift > 15) shift = 15;
 250         else if (shift < -15) shift = -15;
 251         if (shift > 0) {
 252                 for (i = 0; i < shift; i++) {
 253                         result = qm_add16(result, result);
 254                 }
 255         }
 256         else {
 257                 result = result >> (-shift);
 258         }
 259         return result;
 260 }
 261
 262 /*
 263 Description: This function make a 16 bit right shift when shift is +ve.
 264 This function make a 16 bit saturated left shift when shift is -ve. This function
 265 return the result of the shift operation.
 266 */
 267 int16
 268 qm_shr16(int16 op, int shift)
 269 {
 270         return qm_shl16(op, -shift);
 271 }
 272
 273 /*
 274 Description: This function return the number of redundant sign bits in a 16 bit number.
 275 Example: qm_norm16(0x0080) = 7.
 276 */
 277 int16
 278 qm_norm16(int16 op)
 279 {
 280         uint16 u16extraSignBits;
 281         if (op == 0) {
 282                 return 15;
 283         }
 284         else {
 285                 u16extraSignBits = 0;
 286                 while ((op >> 15) == (op >> 14)) {
 287                         u16extraSignBits++;
 288                         op = op << 1;
 289                 }
 290         }
 291         return u16extraSignBits;
 292 }
 293
 294 /*
 295 Description: This function return the number of redundant sign bits in a 32 bit number.
 296 Example: qm_norm32(0x00000080) = 23
 297 */
 298 int16
 299 qm_norm32(int32 op)
 300 {
 301         uint16 u16extraSignBits;
 302         if (op == 0) {
 303                 return 31;
 304         }
 305         else {
 306                 u16extraSignBits = 0;
 307                 while ((op >> 31) == (op >> 30)) {
 308                         u16extraSignBits++;
 309                         op = op << 1;
 310                 }
 311         }
 312         return u16extraSignBits;
 313 }
 314
 315 /*
 316 Description: This function divide two 16 bit unsigned numbers.
 317 The numerator should be less than denominator. So the quotient is always less than 1.
 318 This function return the quotient in q.15 format.
 319 */
 320 int16
 321 qm_div_s(int16 num, int16 denom)
 322 {
 323         int16 var_out;
 324         int16 iteration;
 325         int32 L_num;
 326         int32 L_denom;
 327         L_num = (num) << 15;
 328         L_denom = (denom) << 15;
 329         for (iteration = 0; iteration < 15; iteration++) {
 330                 L_num <<= 1;
 331                 if (L_num >= L_denom) {
 332                         L_num = qm_sub32(L_num, L_denom);
 333                         L_num = qm_add32(L_num, 1);
 334                 }
 335         }
 336         var_out = (int16)(L_num & 0x7fff);
 337         return (var_out);
 338 }
 339
 340 /*
 341 Description: This function compute the absolute value of a 16 bit number.
 342 */
 343 int16
 344 qm_abs16(int16 op)
 345 {
 346         if (op < 0) {
 347                 if (op == (int16)0xffff8000) {
 348                         return 0x7fff;
 349                 } else {
 350                         return -op;
 351                 }
 352         }
 353         else {
 354                 return op;
 355         }
 356 }
 357
 358 /*
 359 Description: This function divide two 16 bit numbers.
 360 The quotient is returned through return value.
 361 The qformat of the quotient is returned through the pointer (qQuotient) passed
 362 to this function. The qformat of quotient is adjusted appropriately such that
 363 the quotient occupies all 16 bits.
 364 */
 365 int16
 366 qm_div16(int16 num, int16 denom, int16 *qQuotient)
 367 {
 368         int16 sign;
 369         int16 nNum, nDenom;
 370         sign = num ^ denom;
 371         num = qm_abs16(num);
 372         denom = qm_abs16(denom);
 373         nNum = qm_norm16(num);
 374         nDenom = qm_norm16(denom);
 375         num = qm_shl16(num, nNum-1);
 376         denom = qm_shl16(denom, nDenom);
 377         *qQuotient = nNum-1-nDenom+15;
 378         if (sign >= 0)  {
 379                 return qm_div_s(num, denom);
 380         }
 381         else {
 382                 return -qm_div_s(num, denom);
 383         }
 384 }
 385
 386 /*
 387 Description: This function compute absolute value of a 32 bit number.
 388 */
 389 int32
 390 qm_abs32(int32 op)
 391 {
 392         if (op < 0) {
 393                 if (op == (int32)0x80000000) {
 394                         return 0x7fffffff;
 395                 } else {
 396                         return -op;
 397                 }
 398         }
 399         else {
 400                 return op;
 401         }
 402 }
 403
 404 /*
 405 Description: This function divide two 32 bit numbers. The division is performed
 406 by considering only important 16 bits in 32 bit numbers.
 407 The quotient is returned through return value.
 408 The qformat of the quotient is returned through the pointer (qquotient) passed
 409 to this function. The qformat of quotient is adjusted appropriately such that
 410 the quotient occupies all 16 bits.
 411 */
 412 int16
 413 qm_div163232(int32 num, int32 denom, int16 *qquotient)
 414 {
 415         int32 sign;
 416         int16 nNum, nDenom;
 417         sign = num ^ denom;
 418         num = qm_abs32(num);
 419         denom = qm_abs32(denom);
 420         nNum = qm_norm32(num);
 421         nDenom = qm_norm32(denom);
 422         num = qm_shl32(num, nNum - 1);
 423         denom = qm_shl32(denom, nDenom);
 424         *qquotient = nNum - 1 - nDenom + 15;
 425         if (sign >= 0)  {
 426                 return qm_div_s((int16)(num >> 16), (int16)(denom >> 16));
 427         }
 428         else {
 429                 return -qm_div_s((int16)(num >> 16), (int16)(denom >> 16));
 430         }
 431 }
 432
 433 /*
 434 Description: This function multiply a 32 bit number with a 16 bit number.
 435 The multiplicaton result is right shifted by 16 bits to fit the result
 436 into 32 bit output.
 437 */
 438 int32
 439 qm_mul323216(int32 op1, int16 op2)
 440 {
 441         int16 hi;
 442         uint16 lo;
 443         int32 result;
 444         hi = op1 >> 16;
 445         lo = (int16)(op1 & 0xffff);
 446         result = qm_mul321616(hi, op2);
 447         result = result + (qm_mulsu321616(op2, lo) >> 16);
 448         return result;
 449 }
 450
 451 /*
 452 Description: This function multiply signed 16 bit number with unsigned 16 bit number and return
 453 the result in 32 bits.
 454 */
 455 int32
 456 qm_mulsu321616(int16 op1, uint16 op2)
 457 {
 458         return (int32)(op1) * op2;
 459 }
 460
 461 /*
 462 Description: This function multiply 32 bit number with 16 bit number. The multiplication result is
 463 right shifted by 15 bits to fit the result into 32 bits. Right shifting by only 15 bits instead of
 464 16 bits is done to remove the extra sign bit formed by multiplication from the return value.
 465 When the input numbers are 0x80000000, 0x8000 the return value is saturated to 0x7fffffff.
 466 */
 467 int32
 468 qm_muls323216(int32 op1, int16 op2)
 469 {
 470         int16 hi;
 471         uint16 lo;
 472         int32 result;
 473         hi = op1 >> 16;
 474         lo = (int16)(op1 & 0xffff);
 475         result = qm_muls321616(hi, op2);
 476         result = qm_add32(result, (qm_mulsu321616(op2, lo) >> 15));
 477         return result;
 478 }
 479
 480 /*
 481 Description: This function multiply two 32 bit numbers. The multiplication result is right
 482 shifted by 32 bits to fit the multiplication result into 32 bits. The right shifted
 483 multiplication result is returned as output.
 484 */
 485 int32
 486 qm_mul32(int32 a, int32 b)
 487 {
 488         int16 hi1, hi2;
 489         uint16 lo1, lo2;
 490         int32 result;
 491         hi1 = a >> 16;
 492         hi2 = b >> 16;
 493         lo1 = (uint16)(a & 0xffff);
 494         lo2 = (uint16)(b & 0xffff);
 495         result = qm_mul321616(hi1, hi2);
 496         result = result + (qm_mulsu321616(hi1, lo2) >> 16);
 497         result = result + (qm_mulsu321616(hi2, lo1) >> 16);
 498         return result;
 499 }
 500
 501 /*
 502 Description: This function multiply two 32 bit numbers. The multiplication result is
 503 right shifted by 31 bits to fit the multiplication result into 32 bits. The right
 504 shifted multiplication result is returned as output. Right shifting by only 31 bits
 505 instead of 32 bits is done to remove the extra sign bit formed by multiplication.
 506 When the input numbers are 0x80000000, 0x80000000 the return value is saturated to
 507 0x7fffffff.
 508 */
 509 int32
 510 qm_muls32(int32 a, int32 b)
 511 {
 512         int16 hi1, hi2;
 513         uint16 lo1, lo2;
 514         int32 result;
 515         hi1 = a >> 16;
 516         hi2 = b >> 16;
 517         lo1 = (uint16)(a & 0xffff);
 518         lo2 = (uint16)(b & 0xffff);
 519         result = qm_muls321616(hi1, hi2);
 520         result = qm_add32(result, (qm_mulsu321616(hi1, lo2) >> 15));
 521         result = qm_add32(result, (qm_mulsu321616(hi2, lo1) >> 15));
 522         result = qm_add32(result, (qm_mulu16(lo1, lo2) >> 15));
 523         return result;
 524 }
 525
 526 /* This table is log2(1+(i/32)) where i=[0:1:31], in q.15 format */
 527 static const int16 log_table[] = {
 528 0,
 529 1455,
 530 2866,
 531 4236,
 532 5568,
 533 6863,
 534 8124,
 535 9352,
 536 10549,
 537 11716,
 538 12855,
 539 13968,
 540 15055,
 541 16117,
 542 17156,
 543 18173,
 544 19168,
 545 20143,
 546 21098,
 547 22034,
 548 22952,
 549 23852,
 550 24736,
 551 25604,
 552 26455,
 553 27292,
 554 28114,
 555 28922,
 556 29717,
 557 30498,
 558 31267,
 559 32024
 560 };
 561 #define LOG_TABLE_SIZE 32               /* log_table size */
 562 #define LOG2_LOG_TABLE_SIZE 5   /* log2(log_table size) */
 563 #define Q_LOG_TABLE 15                  /* qformat of log_table */
 564 #define LOG10_2         19728           /* log10(2) in q.16 */
 565
 566 /*
 567 Description:
 568 This routine takes the input number N and its q format qN and compute
 569 the log10(N). This routine first normalizes the input no N.     Then N is in mag*(2^x) format.
 570 mag is any number in the range 2^30-(2^31 - 1). Then log2(mag * 2^x) = log2(mag) + x is computed.
 571 From that log10(mag * 2^x) = log2(mag * 2^x) * log10(2) is computed.
 572 This routine looks the log2 value in the table considering LOG2_LOG_TABLE_SIZE+1 MSBs.
 573 As the MSB is always 1, only next LOG2_OF_LOG_TABLE_SIZE MSBs are used for table lookup.
 574 Next 16 MSBs are used for interpolation.
 575 Inputs:
 576 N - number to which log10 has to be found.
 577 qN - q format of N
 578 log10N - address where log10(N) will be written.
 579 qLog10N - address where log10N qformat will be written.
 580 Note/Problem:
 581 For accurate results input should be in normalized or near normalized form.
 582 */
 583 void
 584 qm_log10(int32 N, int16 qN, int16 *log10N, int16 *qLog10N)
 585 {
 586         int16 s16norm, s16tableIndex, s16errorApproximation;
 587         uint16 u16offset;
 588         int32 s32log;
 589
 590         /* Logerithm of negative values is undefined.
 591         * assert N is greater than 0.
 592         */
 593         /* ASSERT(N > 0); */
 594
 595         /* normalize the N. */
 596         s16norm = qm_norm32(N);
 597         N = N << s16norm;
 598
 599         /* The qformat of N after normalization.
 600         * -30 is added to treat the no as between 1.0 to 2.0
 601         * i.e. after adding the -30 to the qformat the decimal point will be
 602         * just rigtht of the MSB. (i.e. after sign bit and 1st MSB). i.e.
 603         * at the right side of 30th bit.
 604         */
 605         qN = qN + s16norm - 30;
 606
 607         /* take the table index as the LOG2_OF_LOG_TABLE_SIZE bits right of the MSB */
 608         s16tableIndex = (int16)(N >> (32-(2+LOG2_LOG_TABLE_SIZE)));
 609
 610         /* remove the MSB. the MSB is always 1 after normalization. */
 611         s16tableIndex = s16tableIndex & (int16)((1<<LOG2_LOG_TABLE_SIZE)-1);
 612
 613         /* remove the (1+LOG2_OF_LOG_TABLE_SIZE) MSBs in the N. */
 614         N = N & ((1<<(32-(2+LOG2_LOG_TABLE_SIZE))) - 1);
 615
 616         /* take the offset as the 16 MSBS after table index.
 617         */
 618         u16offset = (uint16)(N >> (32-(2+LOG2_LOG_TABLE_SIZE+16)));
 619
 620         /* look the log value in the table. */
 621         s32log = log_table[s16tableIndex];              /* q.15 format */
 622
 623         /* interpolate using the offset. */
 624         s16errorApproximation = (int16)qm_mulu16(u16offset,
 625                 (uint16)(log_table[s16tableIndex+1]-log_table[s16tableIndex])); /* q.15 */
 626
 627         s32log = qm_add16((int16)s32log, s16errorApproximation);        /* q.15 format */
 628
 629         /* adjust for the qformat of the N as
 630         * log2(mag * 2^x) = log2(mag) + x
 631         */
 632         s32log = qm_add32(s32log, ((int32)-qN)<<15);            /* q.15 format */
 633
 634         /* normalize the result. */
 635         s16norm = qm_norm32(s32log);
 636
 637         /* bring all the important bits into lower 16 bits */
 638         s32log = qm_shl32(s32log, s16norm-16);          /* q.15+s16norm-16 format */
 639
 640         /* compute the log10(N) by multiplying log2(N) with log10(2).
 641         * as log10(mag * 2^x) = log2(mag * 2^x) * log10(2)
 642         * log10N in q.15+s16norm-16+1 (LOG10_2 is in q.16)
 643         */
 644         *log10N = qm_muls16((int16)s32log, (int16)LOG10_2);
 645
 646         /* write the q format of the result. */
 647         *qLog10N = 15 + s16norm - 16 + 1;
 648
 649         return;
 650 }
 651
 652 /*
 653 Description:
 654 This routine compute 1/N.
 655 This routine reformates the given no N as N * 2^qN where N is in between 0.5 and 1.0
 656 in q.15 format in 16 bits. So the problem now boils down to finding the inverse of a
 657 q.15 no in 16 bits which is in the range of 0.5 to 1.0. The output is always between
 658 2.0 to 1. So the output is 2.0 to 1.0 in q.30 format. Once the final output format is found
 659 by taking the qN into account. Inverse is found with newton rapson method. Initially
 660 inverse (x) is guessed as 1/0.75 (with appropriate sign). The new guess is calculated
 661 using the formula x' = 2*x - N*x*x. After 4 or 5 iterations the inverse is very close to
 662 inverse of N.
 663 Inputs:
 664 N - number to which 1/N has to be found.
 665 qn - q format of N.
 666 sqrtN - address where 1/N has to be written.
 667 qsqrtN - address where q format of 1/N has to be written.
 668 */
 669 #define qx 29
 670 void
 671 qm_1byN(int32 N, int16 qN, int32 *result, int16 *qResult)
 672 {
 673         int16 normN;
 674         int32 s32firstTerm, s32secondTerm, x;
 675         int i;
 676
 677         normN = qm_norm32(N);
 678
 679         /* limit N to least significant 16 bits. 15th bit is the sign bit. */
 680         N = qm_shl32(N, normN-16);
 681         qN = qN + normN - 16 - 15;
 682         /* -15 is added to treat N as 16 bit q.15 number in the range from 0.5 to 1 */
 683
 684         /* Take the initial guess as 1/0.75 in qx format with appropriate sign. */
 685         if (N >= 0)
 686         {
 687                 x = (int32)((1/0.75)*(1<<qx));
 688                 /* input no is in the range 0.5 to 1. So 1/0.75 is taken as initial guess. */
 689         }
 690         else
 691         {
 692                 x = (int32)((1/-0.75)*(1<<qx));
 693                 /* input no is in the range -0.5 to -1. So 1/-0.75 is taken as initial guess. */
 694         }
 695
 696         /* iterate the equation x = 2*x - N*x*x for 4 times. */
 697         for (i = 0; i < 4; i++)
 698         {
 699                 s32firstTerm = qm_shl32(x, 1); /* s32firstTerm = 2*x in q.29 */
 700                 s32secondTerm = qm_muls321616((int16)(s32firstTerm>>16), (int16)(s32firstTerm>>16));
 701                 /* s32secondTerm = x*x in q.(29+1-16)*2+1 */
 702                 s32secondTerm = qm_muls321616((int16)(s32secondTerm>>16), (int16)N);
 703                 /* s32secondTerm = N*x*x in q.((29+1-16)*2+1)-16+15+1 i.e. in q.29 */
 704                 x = qm_sub32(s32firstTerm, s32secondTerm);
 705                 /* can be added directly as both are in q.29 */
 706         }
 707
 708         /* Bring the x to q.30 format. */
 709         *result = qm_shl32(x, 1);
 710         /* giving the output in q.30 format for q.15 input in 16 bits. */
 711
 712         /* compute the final q format of the result. */
 713         *qResult = -qN + 30;    /* adjusting the q format of actual output */
 714
 715         return;
 716 }
 717 #undef qx