fixpt31_32.c 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483
  1. /*
  2. * Copyright 2012-15 Advanced Micro Devices, Inc.
  3. *
  4. * Permission is hereby granted, free of charge, to any person obtaining a
  5. * copy of this software and associated documentation files (the "Software"),
  6. * to deal in the Software without restriction, including without limitation
  7. * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  8. * and/or sell copies of the Software, and to permit persons to whom the
  9. * Software is furnished to do so, subject to the following conditions:
  10. *
  11. * The above copyright notice and this permission notice shall be included in
  12. * all copies or substantial portions of the Software.
  13. *
  14. * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  15. * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  16. * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  17. * THE COPYRIGHT HOLDER(S) OR AUTHOR(S) BE LIABLE FOR ANY CLAIM, DAMAGES OR
  18. * OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
  19. * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  20. * OTHER DEALINGS IN THE SOFTWARE.
  21. *
  22. * Authors: AMD
  23. *
  24. */
  25. #include "dm_services.h"
  26. #include "include/fixed31_32.h"
  27. static inline unsigned long long abs_i64(
  28. long long arg)
  29. {
  30. if (arg > 0)
  31. return (unsigned long long)arg;
  32. else
  33. return (unsigned long long)(-arg);
  34. }
  35. /*
  36. * @brief
  37. * result = dividend / divisor
  38. * *remainder = dividend % divisor
  39. */
  40. static inline unsigned long long complete_integer_division_u64(
  41. unsigned long long dividend,
  42. unsigned long long divisor,
  43. unsigned long long *remainder)
  44. {
  45. unsigned long long result;
  46. ASSERT(divisor);
  47. result = div64_u64_rem(dividend, divisor, remainder);
  48. return result;
  49. }
  50. #define FRACTIONAL_PART_MASK \
  51. ((1ULL << FIXED31_32_BITS_PER_FRACTIONAL_PART) - 1)
  52. #define GET_INTEGER_PART(x) \
  53. ((x) >> FIXED31_32_BITS_PER_FRACTIONAL_PART)
  54. #define GET_FRACTIONAL_PART(x) \
  55. (FRACTIONAL_PART_MASK & (x))
  56. struct fixed31_32 dc_fixpt_from_fraction(long long numerator, long long denominator)
  57. {
  58. struct fixed31_32 res;
  59. bool arg1_negative = numerator < 0;
  60. bool arg2_negative = denominator < 0;
  61. unsigned long long arg1_value = arg1_negative ? -numerator : numerator;
  62. unsigned long long arg2_value = arg2_negative ? -denominator : denominator;
  63. unsigned long long remainder;
  64. /* determine integer part */
  65. unsigned long long res_value = complete_integer_division_u64(
  66. arg1_value, arg2_value, &remainder);
  67. ASSERT(res_value <= LONG_MAX);
  68. /* determine fractional part */
  69. {
  70. unsigned int i = FIXED31_32_BITS_PER_FRACTIONAL_PART;
  71. do {
  72. remainder <<= 1;
  73. res_value <<= 1;
  74. if (remainder >= arg2_value) {
  75. res_value |= 1;
  76. remainder -= arg2_value;
  77. }
  78. } while (--i != 0);
  79. }
  80. /* round up LSB */
  81. {
  82. unsigned long long summand = (remainder << 1) >= arg2_value;
  83. ASSERT(res_value <= LLONG_MAX - summand);
  84. res_value += summand;
  85. }
  86. res.value = (long long)res_value;
  87. if (arg1_negative ^ arg2_negative)
  88. res.value = -res.value;
  89. return res;
  90. }
  91. struct fixed31_32 dc_fixpt_mul(struct fixed31_32 arg1, struct fixed31_32 arg2)
  92. {
  93. struct fixed31_32 res;
  94. bool arg1_negative = arg1.value < 0;
  95. bool arg2_negative = arg2.value < 0;
  96. unsigned long long arg1_value = arg1_negative ? -arg1.value : arg1.value;
  97. unsigned long long arg2_value = arg2_negative ? -arg2.value : arg2.value;
  98. unsigned long long arg1_int = GET_INTEGER_PART(arg1_value);
  99. unsigned long long arg2_int = GET_INTEGER_PART(arg2_value);
  100. unsigned long long arg1_fra = GET_FRACTIONAL_PART(arg1_value);
  101. unsigned long long arg2_fra = GET_FRACTIONAL_PART(arg2_value);
  102. unsigned long long tmp;
  103. res.value = arg1_int * arg2_int;
  104. ASSERT(res.value <= LONG_MAX);
  105. res.value <<= FIXED31_32_BITS_PER_FRACTIONAL_PART;
  106. tmp = arg1_int * arg2_fra;
  107. ASSERT(tmp <= (unsigned long long)(LLONG_MAX - res.value));
  108. res.value += tmp;
  109. tmp = arg2_int * arg1_fra;
  110. ASSERT(tmp <= (unsigned long long)(LLONG_MAX - res.value));
  111. res.value += tmp;
  112. tmp = arg1_fra * arg2_fra;
  113. tmp = (tmp >> FIXED31_32_BITS_PER_FRACTIONAL_PART) +
  114. (tmp >= (unsigned long long)dc_fixpt_half.value);
  115. ASSERT(tmp <= (unsigned long long)(LLONG_MAX - res.value));
  116. res.value += tmp;
  117. if (arg1_negative ^ arg2_negative)
  118. res.value = -res.value;
  119. return res;
  120. }
  121. struct fixed31_32 dc_fixpt_sqr(struct fixed31_32 arg)
  122. {
  123. struct fixed31_32 res;
  124. unsigned long long arg_value = abs_i64(arg.value);
  125. unsigned long long arg_int = GET_INTEGER_PART(arg_value);
  126. unsigned long long arg_fra = GET_FRACTIONAL_PART(arg_value);
  127. unsigned long long tmp;
  128. res.value = arg_int * arg_int;
  129. ASSERT(res.value <= LONG_MAX);
  130. res.value <<= FIXED31_32_BITS_PER_FRACTIONAL_PART;
  131. tmp = arg_int * arg_fra;
  132. ASSERT(tmp <= (unsigned long long)(LLONG_MAX - res.value));
  133. res.value += tmp;
  134. ASSERT(tmp <= (unsigned long long)(LLONG_MAX - res.value));
  135. res.value += tmp;
  136. tmp = arg_fra * arg_fra;
  137. tmp = (tmp >> FIXED31_32_BITS_PER_FRACTIONAL_PART) +
  138. (tmp >= (unsigned long long)dc_fixpt_half.value);
  139. ASSERT(tmp <= (unsigned long long)(LLONG_MAX - res.value));
  140. res.value += tmp;
  141. return res;
  142. }
  143. struct fixed31_32 dc_fixpt_recip(struct fixed31_32 arg)
  144. {
  145. /*
  146. * @note
  147. * Good idea to use Newton's method
  148. */
  149. ASSERT(arg.value);
  150. return dc_fixpt_from_fraction(
  151. dc_fixpt_one.value,
  152. arg.value);
  153. }
  154. struct fixed31_32 dc_fixpt_sinc(struct fixed31_32 arg)
  155. {
  156. struct fixed31_32 square;
  157. struct fixed31_32 res = dc_fixpt_one;
  158. int n = 27;
  159. struct fixed31_32 arg_norm = arg;
  160. if (dc_fixpt_le(
  161. dc_fixpt_two_pi,
  162. dc_fixpt_abs(arg))) {
  163. arg_norm = dc_fixpt_sub(
  164. arg_norm,
  165. dc_fixpt_mul_int(
  166. dc_fixpt_two_pi,
  167. (int)div64_s64(
  168. arg_norm.value,
  169. dc_fixpt_two_pi.value)));
  170. }
  171. square = dc_fixpt_sqr(arg_norm);
  172. do {
  173. res = dc_fixpt_sub(
  174. dc_fixpt_one,
  175. dc_fixpt_div_int(
  176. dc_fixpt_mul(
  177. square,
  178. res),
  179. n * (n - 1)));
  180. n -= 2;
  181. } while (n > 2);
  182. if (arg.value != arg_norm.value)
  183. res = dc_fixpt_div(
  184. dc_fixpt_mul(res, arg_norm),
  185. arg);
  186. return res;
  187. }
  188. struct fixed31_32 dc_fixpt_sin(struct fixed31_32 arg)
  189. {
  190. return dc_fixpt_mul(
  191. arg,
  192. dc_fixpt_sinc(arg));
  193. }
  194. struct fixed31_32 dc_fixpt_cos(struct fixed31_32 arg)
  195. {
  196. /* TODO implement argument normalization */
  197. const struct fixed31_32 square = dc_fixpt_sqr(arg);
  198. struct fixed31_32 res = dc_fixpt_one;
  199. int n = 26;
  200. do {
  201. res = dc_fixpt_sub(
  202. dc_fixpt_one,
  203. dc_fixpt_div_int(
  204. dc_fixpt_mul(
  205. square,
  206. res),
  207. n * (n - 1)));
  208. n -= 2;
  209. } while (n != 0);
  210. return res;
  211. }
  212. /*
  213. * @brief
  214. * result = exp(arg),
  215. * where abs(arg) < 1
  216. *
  217. * Calculated as Taylor series.
  218. */
  219. static struct fixed31_32 fixed31_32_exp_from_taylor_series(struct fixed31_32 arg)
  220. {
  221. unsigned int n = 9;
  222. struct fixed31_32 res = dc_fixpt_from_fraction(
  223. n + 2,
  224. n + 1);
  225. /* TODO find correct res */
  226. ASSERT(dc_fixpt_lt(arg, dc_fixpt_one));
  227. do
  228. res = dc_fixpt_add(
  229. dc_fixpt_one,
  230. dc_fixpt_div_int(
  231. dc_fixpt_mul(
  232. arg,
  233. res),
  234. n));
  235. while (--n != 1);
  236. return dc_fixpt_add(
  237. dc_fixpt_one,
  238. dc_fixpt_mul(
  239. arg,
  240. res));
  241. }
  242. struct fixed31_32 dc_fixpt_exp(struct fixed31_32 arg)
  243. {
  244. /*
  245. * @brief
  246. * Main equation is:
  247. * exp(x) = exp(r + m * ln(2)) = (1 << m) * exp(r),
  248. * where m = round(x / ln(2)), r = x - m * ln(2)
  249. */
  250. if (dc_fixpt_le(
  251. dc_fixpt_ln2_div_2,
  252. dc_fixpt_abs(arg))) {
  253. int m = dc_fixpt_round(
  254. dc_fixpt_div(
  255. arg,
  256. dc_fixpt_ln2));
  257. struct fixed31_32 r = dc_fixpt_sub(
  258. arg,
  259. dc_fixpt_mul_int(
  260. dc_fixpt_ln2,
  261. m));
  262. ASSERT(m != 0);
  263. ASSERT(dc_fixpt_lt(
  264. dc_fixpt_abs(r),
  265. dc_fixpt_one));
  266. if (m > 0)
  267. return dc_fixpt_shl(
  268. fixed31_32_exp_from_taylor_series(r),
  269. (unsigned char)m);
  270. else
  271. return dc_fixpt_div_int(
  272. fixed31_32_exp_from_taylor_series(r),
  273. 1LL << -m);
  274. } else if (arg.value != 0)
  275. return fixed31_32_exp_from_taylor_series(arg);
  276. else
  277. return dc_fixpt_one;
  278. }
  279. struct fixed31_32 dc_fixpt_log(struct fixed31_32 arg)
  280. {
  281. struct fixed31_32 res = dc_fixpt_neg(dc_fixpt_one);
  282. /* TODO improve 1st estimation */
  283. struct fixed31_32 error;
  284. ASSERT(arg.value > 0);
  285. /* TODO if arg is negative, return NaN */
  286. /* TODO if arg is zero, return -INF */
  287. do {
  288. struct fixed31_32 res1 = dc_fixpt_add(
  289. dc_fixpt_sub(
  290. res,
  291. dc_fixpt_one),
  292. dc_fixpt_div(
  293. arg,
  294. dc_fixpt_exp(res)));
  295. error = dc_fixpt_sub(
  296. res,
  297. res1);
  298. res = res1;
  299. /* TODO determine max_allowed_error based on quality of exp() */
  300. } while (abs_i64(error.value) > 100ULL);
  301. return res;
  302. }
  303. /* this function is a generic helper to translate fixed point value to
  304. * specified integer format that will consist of integer_bits integer part and
  305. * fractional_bits fractional part. For example it is used in
  306. * dc_fixpt_u2d19 to receive 2 bits integer part and 19 bits fractional
  307. * part in 32 bits. It is used in hw programming (scaler)
  308. */
  309. static inline unsigned int ux_dy(
  310. long long value,
  311. unsigned int integer_bits,
  312. unsigned int fractional_bits)
  313. {
  314. /* 1. create mask of integer part */
  315. unsigned int result = (1 << integer_bits) - 1;
  316. /* 2. mask out fractional part */
  317. unsigned int fractional_part = FRACTIONAL_PART_MASK & value;
  318. /* 3. shrink fixed point integer part to be of integer_bits width*/
  319. result &= GET_INTEGER_PART(value);
  320. /* 4. make space for fractional part to be filled in after integer */
  321. result <<= fractional_bits;
  322. /* 5. shrink fixed point fractional part to of fractional_bits width*/
  323. fractional_part >>= FIXED31_32_BITS_PER_FRACTIONAL_PART - fractional_bits;
  324. /* 6. merge the result */
  325. return result | fractional_part;
  326. }
  327. static inline unsigned int clamp_ux_dy(
  328. long long value,
  329. unsigned int integer_bits,
  330. unsigned int fractional_bits,
  331. unsigned int min_clamp)
  332. {
  333. unsigned int truncated_val = ux_dy(value, integer_bits, fractional_bits);
  334. if (value >= (1LL << (integer_bits + FIXED31_32_BITS_PER_FRACTIONAL_PART)))
  335. return (1 << (integer_bits + fractional_bits)) - 1;
  336. else if (truncated_val > min_clamp)
  337. return truncated_val;
  338. else
  339. return min_clamp;
  340. }
  341. unsigned int dc_fixpt_u3d19(struct fixed31_32 arg)
  342. {
  343. return ux_dy(arg.value, 3, 19);
  344. }
  345. unsigned int dc_fixpt_u2d19(struct fixed31_32 arg)
  346. {
  347. return ux_dy(arg.value, 2, 19);
  348. }
  349. unsigned int dc_fixpt_u0d19(struct fixed31_32 arg)
  350. {
  351. return ux_dy(arg.value, 0, 19);
  352. }
  353. unsigned int dc_fixpt_clamp_u0d14(struct fixed31_32 arg)
  354. {
  355. return clamp_ux_dy(arg.value, 0, 14, 1);
  356. }
  357. unsigned int dc_fixpt_clamp_u0d10(struct fixed31_32 arg)
  358. {
  359. return clamp_ux_dy(arg.value, 0, 10, 1);
  360. }
  361. int dc_fixpt_s4d19(struct fixed31_32 arg)
  362. {
  363. if (arg.value < 0)
  364. return -(int)ux_dy(dc_fixpt_abs(arg).value, 4, 19);
  365. else
  366. return ux_dy(arg.value, 4, 19);
  367. }