mp_flipper_polyfill.c 8.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398
  1. #include <limits.h>
  2. #include <stdint.h>
  3. #include <float.h>
  4. #include <math.h>
  5. #include <string.h>
  6. #define POLYFILL_FUN_1(name, ret, arg1) ret name(arg1)
  7. #define POLYFILL_FUN_2(name, ret, arg1, arg2) ret name(arg1, arg2)
  8. #define POLYFILL_FUN_3(name, ret, arg1, arg2, arg3) ret name(arg1, arg2, arg3)
  9. #define POLYFILL_FUN_4(name, ret, arg1, arg2, arg3, arg4) ret name(arg1, arg2, arg3, arg4)
  10. #ifndef __aeabi_dcmple
  11. POLYFILL_FUN_2(__aeabi_dcmple, double, double, double) {
  12. }
  13. #endif
  14. #ifndef __aeabi_dcmplt
  15. POLYFILL_FUN_2(__aeabi_dcmplt, double, double, double) {
  16. }
  17. #endif
  18. #ifndef __aeabi_dcmpun
  19. POLYFILL_FUN_2(__aeabi_dcmpun, double, double, double) {
  20. }
  21. #endif
  22. #ifndef __aeabi_dcmpeq
  23. POLYFILL_FUN_2(__aeabi_dcmpeq, double, double, double) {
  24. }
  25. #endif
  26. #ifndef __aeabi_dmul
  27. double __aeabi_dmul(double x, double y) {
  28. return x * y;
  29. }
  30. #endif
  31. #ifndef __aeabi_dadd
  32. POLYFILL_FUN_2(__aeabi_dadd, double, double x, double y) {
  33. return x + y;
  34. }
  35. #endif
  36. #ifndef __aeabi_ddiv
  37. POLYFILL_FUN_2(__aeabi_ddiv, double, double x, double y) {
  38. return x / y;
  39. }
  40. #endif
  41. #ifndef __aeabi_l2d
  42. POLYFILL_FUN_1(__aeabi_l2d, double, long x) {
  43. return x;
  44. }
  45. #endif
  46. #ifndef __aeabi_f2d
  47. POLYFILL_FUN_1(__aeabi_f2d, double, float x) {
  48. return x;
  49. }
  50. #endif
  51. #ifndef __aeabi_dsub
  52. POLYFILL_FUN_2(__aeabi_dsub, double, double x, double y) {
  53. return x - y;
  54. }
  55. #endif
  56. #ifndef __aeabi_dcmpge
  57. POLYFILL_FUN_2(__aeabi_dcmpge, double, double, double) {
  58. }
  59. #endif
  60. #ifndef __aeabi_i2d
  61. POLYFILL_FUN_1(__aeabi_i2d, double, int x) {
  62. return x;
  63. }
  64. #endif
  65. #ifndef __aeabi_dcmpgt
  66. POLYFILL_FUN_1(__aeabi_dcmpgt, double, double) {
  67. }
  68. #endif
  69. #ifndef __aeabi_d2iz
  70. POLYFILL_FUN_1(__aeabi_d2iz, long int, double x) {
  71. return x;
  72. }
  73. #endif
  74. #ifndef __aeabi_d2lz
  75. POLYFILL_FUN_1(__aeabi_d2lz, long, double x) {
  76. return x;
  77. }
  78. #endif
  79. #ifndef __aeabi_d2uiz
  80. POLYFILL_FUN_1(__aeabi_d2uiz, unsigned long int, double x) {
  81. return x;
  82. }
  83. #endif
  84. #ifndef __aeabi_d2f
  85. POLYFILL_FUN_1(__aeabi_d2f, float, double x) {
  86. return x;
  87. }
  88. #endif
  89. #ifndef __aeabi_ldivmod
  90. POLYFILL_FUN_2(__aeabi_ldivmod, long, long, long) {
  91. }
  92. #endif
  93. #ifndef strtox
  94. POLYFILL_FUN_4(strtox, long long unsigned int, const char*, char**, int, long long unsigned int) {
  95. }
  96. #endif
  97. #if FLT_EVAL_METHOD == 0 || FLT_EVAL_METHOD == 1
  98. #define EPS DBL_EPSILON
  99. #elif FLT_EVAL_METHOD == 2
  100. #define EPS LDBL_EPSILON
  101. #endif
  102. static const double_t toint = 1 / EPS;
  103. #ifndef FORCE_EVAL
  104. #define FORCE_EVAL(x) \
  105. do { \
  106. if(sizeof(x) == sizeof(float)) { \
  107. volatile float __x; \
  108. __x = (x); \
  109. (void)__x; \
  110. } else if(sizeof(x) == sizeof(double)) { \
  111. volatile double __x; \
  112. __x = (x); \
  113. (void)__x; \
  114. } else { \
  115. volatile long double __x; \
  116. __x = (x); \
  117. (void)__x; \
  118. } \
  119. } while(0)
  120. #endif
  121. #ifndef floorf
  122. float floorf(float x) {
  123. union {
  124. float f;
  125. uint32_t i;
  126. } u = {x};
  127. int e = (int)(u.i >> 23 & 0xff) - 0x7f;
  128. uint32_t m;
  129. if(e >= 23) return x;
  130. if(e >= 0) {
  131. m = 0x007fffff >> e;
  132. if((u.i & m) == 0) return x;
  133. FORCE_EVAL(x + 0x1p120f);
  134. if(u.i >> 31) u.i += m;
  135. u.i &= ~m;
  136. } else {
  137. FORCE_EVAL(x + 0x1p120f);
  138. if(u.i >> 31 == 0)
  139. u.i = 0;
  140. else if(u.i << 1)
  141. u.f = -1.0;
  142. }
  143. return u.f;
  144. }
  145. #endif
  146. #ifndef floor
  147. double floor(double x) {
  148. union {
  149. double f;
  150. uint64_t i;
  151. } u = {x};
  152. int e = u.i >> 52 & 0x7ff;
  153. double_t y;
  154. if(e >= 0x3ff + 52 || x == 0) return x;
  155. /* y = int(x) - x, where int(x) is an integer neighbor of x */
  156. if(u.i >> 63)
  157. y = x - toint + toint - x;
  158. else
  159. y = x + toint - toint - x;
  160. /* special case because of non-nearest rounding modes */
  161. if(e <= 0x3ff - 1) {
  162. FORCE_EVAL(y);
  163. return u.i >> 63 ? -1 : 0;
  164. }
  165. if(y > 0) return x + y - 1;
  166. return x + y;
  167. }
  168. #endif
  169. #ifndef fmodf
  170. float fmodf(float x, float y) {
  171. union {
  172. float f;
  173. uint32_t i;
  174. } ux = {x}, uy = {y};
  175. int ex = ux.i >> 23 & 0xff;
  176. int ey = uy.i >> 23 & 0xff;
  177. uint32_t sx = ux.i & 0x80000000;
  178. uint32_t i;
  179. uint32_t uxi = ux.i;
  180. if(uy.i << 1 == 0 || isnan(y) || ex == 0xff) return (x * y) / (x * y);
  181. if(uxi << 1 <= uy.i << 1) {
  182. if(uxi << 1 == uy.i << 1) return 0 * x;
  183. return x;
  184. }
  185. /* normalize x and y */
  186. if(!ex) {
  187. for(i = uxi << 9; i >> 31 == 0; ex--, i <<= 1)
  188. ;
  189. uxi <<= -ex + 1;
  190. } else {
  191. uxi &= -1U >> 9;
  192. uxi |= 1U << 23;
  193. }
  194. if(!ey) {
  195. for(i = uy.i << 9; i >> 31 == 0; ey--, i <<= 1)
  196. ;
  197. uy.i <<= -ey + 1;
  198. } else {
  199. uy.i &= -1U >> 9;
  200. uy.i |= 1U << 23;
  201. }
  202. /* x mod y */
  203. for(; ex > ey; ex--) {
  204. i = uxi - uy.i;
  205. if(i >> 31 == 0) {
  206. if(i == 0) return 0 * x;
  207. uxi = i;
  208. }
  209. uxi <<= 1;
  210. }
  211. i = uxi - uy.i;
  212. if(i >> 31 == 0) {
  213. if(i == 0) return 0 * x;
  214. uxi = i;
  215. }
  216. for(; uxi >> 23 == 0; uxi <<= 1, ex--)
  217. ;
  218. /* scale result up */
  219. if(ex > 0) {
  220. uxi -= 1U << 23;
  221. uxi |= (uint32_t)ex << 23;
  222. } else {
  223. uxi >>= -ex + 1;
  224. }
  225. uxi |= sx;
  226. ux.i = uxi;
  227. return ux.f;
  228. }
  229. #endif
  230. #ifndef fmod
  231. double fmod(double x, double y) {
  232. union {
  233. double f;
  234. uint64_t i;
  235. } ux = {x}, uy = {y};
  236. int ex = ux.i >> 52 & 0x7ff;
  237. int ey = uy.i >> 52 & 0x7ff;
  238. int sx = ux.i >> 63;
  239. uint64_t i;
  240. /* in the followings uxi should be ux.i, but then gcc wrongly adds */
  241. /* float load/store to inner loops ruining performance and code size */
  242. uint64_t uxi = ux.i;
  243. if(uy.i << 1 == 0 || isnan(y) || ex == 0x7ff) return (x * y) / (x * y);
  244. if(uxi << 1 <= uy.i << 1) {
  245. if(uxi << 1 == uy.i << 1) return 0 * x;
  246. return x;
  247. }
  248. /* normalize x and y */
  249. if(!ex) {
  250. for(i = uxi << 12; i >> 63 == 0; ex--, i <<= 1)
  251. ;
  252. uxi <<= -ex + 1;
  253. } else {
  254. uxi &= -1ULL >> 12;
  255. uxi |= 1ULL << 52;
  256. }
  257. if(!ey) {
  258. for(i = uy.i << 12; i >> 63 == 0; ey--, i <<= 1)
  259. ;
  260. uy.i <<= -ey + 1;
  261. } else {
  262. uy.i &= -1ULL >> 12;
  263. uy.i |= 1ULL << 52;
  264. }
  265. /* x mod y */
  266. for(; ex > ey; ex--) {
  267. i = uxi - uy.i;
  268. if(i >> 63 == 0) {
  269. if(i == 0) return 0 * x;
  270. uxi = i;
  271. }
  272. uxi <<= 1;
  273. }
  274. i = uxi - uy.i;
  275. if(i >> 63 == 0) {
  276. if(i == 0) return 0 * x;
  277. uxi = i;
  278. }
  279. for(; uxi >> 52 == 0; uxi <<= 1, ex--)
  280. ;
  281. /* scale result */
  282. if(ex > 0) {
  283. uxi -= 1ULL << 52;
  284. uxi |= (uint64_t)ex << 52;
  285. } else {
  286. uxi >>= -ex + 1;
  287. }
  288. uxi |= (uint64_t)sx << 63;
  289. ux.i = uxi;
  290. return ux.f;
  291. }
  292. #endif
  293. #ifndef nearbyintf
  294. float nearbyintf(float x) {
  295. union {
  296. float f;
  297. uint32_t i;
  298. } u = {x};
  299. int e = u.i >> 23 & 0xff;
  300. int s = u.i >> 31;
  301. float_t y;
  302. if(e >= 0x7f + 23) return x;
  303. if(s)
  304. y = x - 0x1p23f + 0x1p23f;
  305. else
  306. y = x + 0x1p23f - 0x1p23f;
  307. if(y == 0) return s ? -0.0f : 0.0f;
  308. return y;
  309. }
  310. #endif
  311. #ifndef stroll
  312. long long strtoll(const char* restrict s, char** restrict p, int base) {
  313. return strtox(s, p, base, LLONG_MIN);
  314. }
  315. #endif
  316. #ifndef pow
  317. double pow(double x, double y) {
  318. return powf(x, y);
  319. }
  320. #endif
  321. #ifndef rint
  322. double rint(double x) {
  323. union {
  324. double f;
  325. uint64_t i;
  326. } u = {x};
  327. int e = u.i >> 52 & 0x7ff;
  328. int s = u.i >> 63;
  329. double_t y;
  330. if(e >= 0x3ff + 52) return x;
  331. if(s)
  332. y = x - toint + toint;
  333. else
  334. y = x + toint - toint;
  335. if(y == 0) return s ? -0.0 : 0;
  336. return y;
  337. }
  338. #endif
  339. #ifndef nearbyint
  340. double nearbyint(double x) {
  341. #ifdef FE_INEXACT
  342. #pragma STDC FENV_ACCESS ON
  343. int e;
  344. e = fetestexcept(FE_INEXACT);
  345. #endif
  346. x = rint(x);
  347. #ifdef FE_INEXACT
  348. if(!e) feclearexcept(FE_INEXACT);
  349. #endif
  350. return x;
  351. }
  352. #endif