arm_sqrt_q31.c 4.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143
  1. /* ----------------------------------------------------------------------
  2. * Project: CMSIS DSP Library
  3. * Title: arm_sqrt_q31.c
  4. * Description: Q31 square root function
  5. *
  6. * $Date: 27. January 2017
  7. * $Revision: V.1.5.1
  8. *
  9. * Target Processor: Cortex-M cores
  10. * -------------------------------------------------------------------- */
  11. /*
  12. * Copyright (C) 2010-2017 ARM Limited or its affiliates. All rights reserved.
  13. *
  14. * SPDX-License-Identifier: Apache-2.0
  15. *
  16. * Licensed under the Apache License, Version 2.0 (the License); you may
  17. * not use this file except in compliance with the License.
  18. * You may obtain a copy of the License at
  19. *
  20. * www.apache.org/licenses/LICENSE-2.0
  21. *
  22. * Unless required by applicable law or agreed to in writing, software
  23. * distributed under the License is distributed on an AS IS BASIS, WITHOUT
  24. * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  25. * See the License for the specific language governing permissions and
  26. * limitations under the License.
  27. */
  28. #include "arm_math.h"
  29. #include "arm_common_tables.h"
  30. /**
  31. * @ingroup groupFastMath
  32. */
  33. /**
  34. * @addtogroup SQRT
  35. * @{
  36. */
  37. /**
  38. * @brief Q31 square root function.
  39. * @param[in] in input value. The range of the input value is [0 +1) or 0x00000000 to 0x7FFFFFFF.
  40. * @param[out] *pOut square root of input value.
  41. * @return The function returns ARM_MATH_SUCCESS if the input value is positive
  42. * and ARM_MATH_ARGUMENT_ERROR if the input is negative. For
  43. * negative inputs, the function returns *pOut = 0.
  44. */
  45. arm_status arm_sqrt_q31(
  46. q31_t in,
  47. q31_t * pOut)
  48. {
  49. q31_t number, temp1, bits_val1, var1, signBits1, half;
  50. float32_t temp_float1;
  51. union
  52. {
  53. q31_t fracval;
  54. float32_t floatval;
  55. } tempconv;
  56. number = in;
  57. /* If the input is a positive number then compute the signBits. */
  58. if (number > 0)
  59. {
  60. signBits1 = __CLZ(number) - 1;
  61. /* Shift by the number of signBits1 */
  62. if ((signBits1 % 2) == 0)
  63. {
  64. number = number << signBits1;
  65. }
  66. else
  67. {
  68. number = number << (signBits1 - 1);
  69. }
  70. /* Calculate half value of the number */
  71. half = number >> 1;
  72. /* Store the number for later use */
  73. temp1 = number;
  74. /*Convert to float */
  75. temp_float1 = number * 4.6566128731e-010f;
  76. /*Store as integer */
  77. tempconv.floatval = temp_float1;
  78. bits_val1 = tempconv.fracval;
  79. /* Subtract the shifted value from the magic number to give intial guess */
  80. bits_val1 = 0x5f3759df - (bits_val1 >> 1); /* gives initial guess */
  81. /* Store as float */
  82. tempconv.fracval = bits_val1;
  83. temp_float1 = tempconv.floatval;
  84. /* Convert to integer format */
  85. var1 = (q31_t) (temp_float1 * 1073741824);
  86. /* 1st iteration */
  87. var1 = ((q31_t) ((q63_t) var1 * (0x30000000 -
  88. ((q31_t)
  89. ((((q31_t)
  90. (((q63_t) var1 * var1) >> 31)) *
  91. (q63_t) half) >> 31))) >> 31)) << 2;
  92. /* 2nd iteration */
  93. var1 = ((q31_t) ((q63_t) var1 * (0x30000000 -
  94. ((q31_t)
  95. ((((q31_t)
  96. (((q63_t) var1 * var1) >> 31)) *
  97. (q63_t) half) >> 31))) >> 31)) << 2;
  98. /* 3rd iteration */
  99. var1 = ((q31_t) ((q63_t) var1 * (0x30000000 -
  100. ((q31_t)
  101. ((((q31_t)
  102. (((q63_t) var1 * var1) >> 31)) *
  103. (q63_t) half) >> 31))) >> 31)) << 2;
  104. /* Multiply the inverse square root with the original value */
  105. var1 = ((q31_t) (((q63_t) temp1 * var1) >> 31)) << 1;
  106. /* Shift the output down accordingly */
  107. if ((signBits1 % 2) == 0)
  108. {
  109. var1 = var1 >> (signBits1 / 2);
  110. }
  111. else
  112. {
  113. var1 = var1 >> ((signBits1 - 1) / 2);
  114. }
  115. *pOut = var1;
  116. return (ARM_MATH_SUCCESS);
  117. }
  118. /* If the number is a negative number then store zero as its square root value */
  119. else
  120. {
  121. *pOut = 0;
  122. return (ARM_MATH_ARGUMENT_ERROR);
  123. }
  124. }
  125. /**
  126. * @} end of SQRT group
  127. */