phy_qmath.c 7.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310
  1. /*
  2. * Copyright (c) 2010 Broadcom Corporation
  3. *
  4. * Permission to use, copy, modify, and/or distribute this software for any
  5. * purpose with or without fee is hereby granted, provided that the above
  6. * copyright notice and this permission notice appear in all copies.
  7. *
  8. * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  9. * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
  10. * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
  11. * SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  12. * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION
  13. * OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
  14. * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  15. */
  16. #include "phy_qmath.h"
  17. /*
  18. * Description: This function make 16 bit unsigned multiplication.
  19. * To fit the output into 16 bits the 32 bit multiplication result is right
  20. * shifted by 16 bits.
  21. */
  22. u16 qm_mulu16(u16 op1, u16 op2)
  23. {
  24. return (u16) (((u32) op1 * (u32) op2) >> 16);
  25. }
  26. /*
  27. * Description: This function make 16 bit multiplication and return the result
  28. * in 16 bits. To fit the multiplication result into 16 bits the multiplication
  29. * result is right shifted by 15 bits. Right shifting 15 bits instead of 16 bits
  30. * is done to remove the extra sign bit formed due to the multiplication.
  31. * When both the 16bit inputs are 0x8000 then the output is saturated to
  32. * 0x7fffffff.
  33. */
  34. s16 qm_muls16(s16 op1, s16 op2)
  35. {
  36. s32 result;
  37. if (op1 == (s16) 0x8000 && op2 == (s16) 0x8000)
  38. result = 0x7fffffff;
  39. else
  40. result = ((s32) (op1) * (s32) (op2));
  41. return (s16) (result >> 15);
  42. }
  43. /*
  44. * Description: This function add two 32 bit numbers and return the 32bit
  45. * result. If the result overflow 32 bits, the output will be saturated to
  46. * 32bits.
  47. */
  48. s32 qm_add32(s32 op1, s32 op2)
  49. {
  50. s32 result;
  51. result = op1 + op2;
  52. if (op1 < 0 && op2 < 0 && result > 0)
  53. result = 0x80000000;
  54. else if (op1 > 0 && op2 > 0 && result < 0)
  55. result = 0x7fffffff;
  56. return result;
  57. }
  58. /*
  59. * Description: This function add two 16 bit numbers and return the 16bit
  60. * result. If the result overflow 16 bits, the output will be saturated to
  61. * 16bits.
  62. */
  63. s16 qm_add16(s16 op1, s16 op2)
  64. {
  65. s16 result;
  66. s32 temp = (s32) op1 + (s32) op2;
  67. if (temp > (s32) 0x7fff)
  68. result = (s16) 0x7fff;
  69. else if (temp < (s32) 0xffff8000)
  70. result = (s16) 0xffff8000;
  71. else
  72. result = (s16) temp;
  73. return result;
  74. }
  75. /*
  76. * Description: This function make 16 bit subtraction and return the 16bit
  77. * result. If the result overflow 16 bits, the output will be saturated to
  78. * 16bits.
  79. */
  80. s16 qm_sub16(s16 op1, s16 op2)
  81. {
  82. s16 result;
  83. s32 temp = (s32) op1 - (s32) op2;
  84. if (temp > (s32) 0x7fff)
  85. result = (s16) 0x7fff;
  86. else if (temp < (s32) 0xffff8000)
  87. result = (s16) 0xffff8000;
  88. else
  89. result = (s16) temp;
  90. return result;
  91. }
  92. /*
  93. * Description: This function make a 32 bit saturated left shift when the
  94. * specified shift is +ve. This function will make a 32 bit right shift when
  95. * the specified shift is -ve. This function return the result after shifting
  96. * operation.
  97. */
  98. s32 qm_shl32(s32 op, int shift)
  99. {
  100. int i;
  101. s32 result;
  102. result = op;
  103. if (shift > 31)
  104. shift = 31;
  105. else if (shift < -31)
  106. shift = -31;
  107. if (shift >= 0) {
  108. for (i = 0; i < shift; i++)
  109. result = qm_add32(result, result);
  110. } else {
  111. result = result >> (-shift);
  112. }
  113. return result;
  114. }
  115. /*
  116. * Description: This function make a 16 bit saturated left shift when the
  117. * specified shift is +ve. This function will make a 16 bit right shift when
  118. * the specified shift is -ve. This function return the result after shifting
  119. * operation.
  120. */
  121. s16 qm_shl16(s16 op, int shift)
  122. {
  123. int i;
  124. s16 result;
  125. result = op;
  126. if (shift > 15)
  127. shift = 15;
  128. else if (shift < -15)
  129. shift = -15;
  130. if (shift > 0) {
  131. for (i = 0; i < shift; i++)
  132. result = qm_add16(result, result);
  133. } else {
  134. result = result >> (-shift);
  135. }
  136. return result;
  137. }
  138. /*
  139. * Description: This function make a 16 bit right shift when shift is +ve.
  140. * This function make a 16 bit saturated left shift when shift is -ve. This
  141. * function return the result of the shift operation.
  142. */
  143. s16 qm_shr16(s16 op, int shift)
  144. {
  145. return qm_shl16(op, -shift);
  146. }
  147. /*
  148. * Description: This function return the number of redundant sign bits in a
  149. * 32 bit number. Example: qm_norm32(0x00000080) = 23
  150. */
  151. s16 qm_norm32(s32 op)
  152. {
  153. u16 u16extraSignBits;
  154. if (op == 0) {
  155. return 31;
  156. } else {
  157. u16extraSignBits = 0;
  158. while ((op >> 31) == (op >> 30)) {
  159. u16extraSignBits++;
  160. op = op << 1;
  161. }
  162. }
  163. return u16extraSignBits;
  164. }
  165. /* This table is log2(1+(i/32)) where i=[0:1:32], in q.15 format */
  166. static const s16 log_table[] = {
  167. 0,
  168. 1455,
  169. 2866,
  170. 4236,
  171. 5568,
  172. 6863,
  173. 8124,
  174. 9352,
  175. 10549,
  176. 11716,
  177. 12855,
  178. 13968,
  179. 15055,
  180. 16117,
  181. 17156,
  182. 18173,
  183. 19168,
  184. 20143,
  185. 21098,
  186. 22034,
  187. 22952,
  188. 23852,
  189. 24736,
  190. 25604,
  191. 26455,
  192. 27292,
  193. 28114,
  194. 28922,
  195. 29717,
  196. 30498,
  197. 31267,
  198. 32024,
  199. 32767
  200. };
  201. #define LOG_TABLE_SIZE 32 /* log_table size */
  202. #define LOG2_LOG_TABLE_SIZE 5 /* log2(log_table size) */
  203. #define Q_LOG_TABLE 15 /* qformat of log_table */
  204. #define LOG10_2 19728 /* log10(2) in q.16 */
  205. /*
  206. * Description:
  207. * This routine takes the input number N and its q format qN and compute
  208. * the log10(N). This routine first normalizes the input no N. Then N is in
  209. * mag*(2^x) format. mag is any number in the range 2^30-(2^31 - 1).
  210. * Then log2(mag * 2^x) = log2(mag) + x is computed. From that
  211. * log10(mag * 2^x) = log2(mag * 2^x) * log10(2) is computed.
  212. * This routine looks the log2 value in the table considering
  213. * LOG2_LOG_TABLE_SIZE+1 MSBs. As the MSB is always 1, only next
  214. * LOG2_OF_LOG_TABLE_SIZE MSBs are used for table lookup. Next 16 MSBs are used
  215. * for interpolation.
  216. * Inputs:
  217. * N - number to which log10 has to be found.
  218. * qN - q format of N
  219. * log10N - address where log10(N) will be written.
  220. * qLog10N - address where log10N qformat will be written.
  221. * Note/Problem:
  222. * For accurate results input should be in normalized or near normalized form.
  223. */
  224. void qm_log10(s32 N, s16 qN, s16 *log10N, s16 *qLog10N)
  225. {
  226. s16 s16norm, s16tableIndex, s16errorApproximation;
  227. u16 u16offset;
  228. s32 s32log;
  229. /* normalize the N. */
  230. s16norm = qm_norm32(N);
  231. N = N << s16norm;
  232. /* The qformat of N after normalization.
  233. * -30 is added to treat the no as between 1.0 to 2.0
  234. * i.e. after adding the -30 to the qformat the decimal point will be
  235. * just rigtht of the MSB. (i.e. after sign bit and 1st MSB). i.e.
  236. * at the right side of 30th bit.
  237. */
  238. qN = qN + s16norm - 30;
  239. /* take the table index as the LOG2_OF_LOG_TABLE_SIZE bits right of the
  240. * MSB */
  241. s16tableIndex = (s16) (N >> (32 - (2 + LOG2_LOG_TABLE_SIZE)));
  242. /* remove the MSB. the MSB is always 1 after normalization. */
  243. s16tableIndex =
  244. s16tableIndex & (s16) ((1 << LOG2_LOG_TABLE_SIZE) - 1);
  245. /* remove the (1+LOG2_OF_LOG_TABLE_SIZE) MSBs in the N. */
  246. N = N & ((1 << (32 - (2 + LOG2_LOG_TABLE_SIZE))) - 1);
  247. /* take the offset as the 16 MSBS after table index.
  248. */
  249. u16offset = (u16) (N >> (32 - (2 + LOG2_LOG_TABLE_SIZE + 16)));
  250. /* look the log value in the table. */
  251. s32log = log_table[s16tableIndex]; /* q.15 format */
  252. /* interpolate using the offset. q.15 format. */
  253. s16errorApproximation = (s16) qm_mulu16(u16offset,
  254. (u16) (log_table[s16tableIndex + 1] -
  255. log_table[s16tableIndex]));
  256. /* q.15 format */
  257. s32log = qm_add16((s16) s32log, s16errorApproximation);
  258. /* adjust for the qformat of the N as
  259. * log2(mag * 2^x) = log2(mag) + x
  260. */
  261. s32log = qm_add32(s32log, ((s32) -qN) << 15); /* q.15 format */
  262. /* normalize the result. */
  263. s16norm = qm_norm32(s32log);
  264. /* bring all the important bits into lower 16 bits */
  265. /* q.15+s16norm-16 format */
  266. s32log = qm_shl32(s32log, s16norm - 16);
  267. /* compute the log10(N) by multiplying log2(N) with log10(2).
  268. * as log10(mag * 2^x) = log2(mag * 2^x) * log10(2)
  269. * log10N in q.15+s16norm-16+1 (LOG10_2 is in q.16)
  270. */
  271. *log10N = qm_muls16((s16) s32log, (s16) LOG10_2);
  272. /* write the q format of the result. */
  273. *qLog10N = 15 + s16norm - 16 + 1;
  274. return;
  275. }