|
@@ -4,89 +4,115 @@
|
|
|
#include <math.h>
|
|
#include <math.h>
|
|
|
#include <string.h>
|
|
#include <string.h>
|
|
|
|
|
|
|
|
-#define POLYFILL_FUN_1(name, ret, arg1) \
|
|
|
|
|
- ret name(arg1) {}
|
|
|
|
|
-#define POLYFILL_FUN_2(name, ret, arg1, arg2) \
|
|
|
|
|
- ret name(arg1, arg2) {}
|
|
|
|
|
-#define POLYFILL_FUN_3(name, ret, arg1, arg2, arg3) \
|
|
|
|
|
- ret name(arg1, arg2, arg3) {}
|
|
|
|
|
-#define POLYFILL_FUN_4(name, ret, arg1, arg2, arg3, arg4) \
|
|
|
|
|
- ret name(arg1, arg2, arg3, arg4) {}
|
|
|
|
|
|
|
+#define POLYFILL_FUN_1(name, ret, arg1) ret name(arg1)
|
|
|
|
|
+#define POLYFILL_FUN_2(name, ret, arg1, arg2) ret name(arg1, arg2)
|
|
|
|
|
+#define POLYFILL_FUN_3(name, ret, arg1, arg2, arg3) ret name(arg1, arg2, arg3)
|
|
|
|
|
+#define POLYFILL_FUN_4(name, ret, arg1, arg2, arg3, arg4) ret name(arg1, arg2, arg3, arg4)
|
|
|
|
|
|
|
|
#ifndef __aeabi_dcmple
|
|
#ifndef __aeabi_dcmple
|
|
|
-POLYFILL_FUN_2(__aeabi_dcmple, double, double, double)
|
|
|
|
|
|
|
+POLYFILL_FUN_2(__aeabi_dcmple, double, double, double) {
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_dcmplt
|
|
#ifndef __aeabi_dcmplt
|
|
|
-POLYFILL_FUN_2(__aeabi_dcmplt, double, double, double)
|
|
|
|
|
|
|
+POLYFILL_FUN_2(__aeabi_dcmplt, double, double, double) {
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_dcmpun
|
|
#ifndef __aeabi_dcmpun
|
|
|
-POLYFILL_FUN_2(__aeabi_dcmpun, double, double, double)
|
|
|
|
|
|
|
+POLYFILL_FUN_2(__aeabi_dcmpun, double, double, double) {
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_dcmpeq
|
|
#ifndef __aeabi_dcmpeq
|
|
|
-POLYFILL_FUN_2(__aeabi_dcmpeq, double, double, double)
|
|
|
|
|
|
|
+POLYFILL_FUN_2(__aeabi_dcmpeq, double, double, double) {
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_dmul
|
|
#ifndef __aeabi_dmul
|
|
|
-POLYFILL_FUN_2(__aeabi_dmul, double, double, double)
|
|
|
|
|
|
|
+double __aeabi_dmul(double x, double y) {
|
|
|
|
|
+ return x * y;
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_dadd
|
|
#ifndef __aeabi_dadd
|
|
|
-POLYFILL_FUN_2(__aeabi_dadd, double, double, double)
|
|
|
|
|
|
|
+POLYFILL_FUN_2(__aeabi_dadd, double, double x, double y) {
|
|
|
|
|
+ return x + y;
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_ddiv
|
|
#ifndef __aeabi_ddiv
|
|
|
-POLYFILL_FUN_2(__aeabi_ddiv, double, double, double)
|
|
|
|
|
|
|
+POLYFILL_FUN_2(__aeabi_ddiv, double, double x, double y) {
|
|
|
|
|
+ return x / y;
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_l2d
|
|
#ifndef __aeabi_l2d
|
|
|
-POLYFILL_FUN_1(__aeabi_l2d, double, long)
|
|
|
|
|
|
|
+POLYFILL_FUN_1(__aeabi_l2d, double, long x) {
|
|
|
|
|
+ return x;
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_f2d
|
|
#ifndef __aeabi_f2d
|
|
|
-POLYFILL_FUN_1(__aeabi_f2d, double, float)
|
|
|
|
|
|
|
+POLYFILL_FUN_1(__aeabi_f2d, double, float x) {
|
|
|
|
|
+ return x;
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_dsub
|
|
#ifndef __aeabi_dsub
|
|
|
-POLYFILL_FUN_2(__aeabi_dsub, double, double, double)
|
|
|
|
|
|
|
+POLYFILL_FUN_2(__aeabi_dsub, double, double x, double y) {
|
|
|
|
|
+ return x - y;
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_dcmpge
|
|
#ifndef __aeabi_dcmpge
|
|
|
-POLYFILL_FUN_2(__aeabi_dcmpge, double, double, double)
|
|
|
|
|
|
|
+POLYFILL_FUN_2(__aeabi_dcmpge, double, double, double) {
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_i2d
|
|
#ifndef __aeabi_i2d
|
|
|
-POLYFILL_FUN_1(__aeabi_i2d, double, long)
|
|
|
|
|
|
|
+POLYFILL_FUN_1(__aeabi_i2d, double, int x) {
|
|
|
|
|
+ return x;
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_dcmpgt
|
|
#ifndef __aeabi_dcmpgt
|
|
|
-POLYFILL_FUN_1(__aeabi_dcmpgt, double, double)
|
|
|
|
|
|
|
+POLYFILL_FUN_1(__aeabi_dcmpgt, double, double) {
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_d2iz
|
|
#ifndef __aeabi_d2iz
|
|
|
-POLYFILL_FUN_1(__aeabi_d2iz, long int, double)
|
|
|
|
|
|
|
+POLYFILL_FUN_1(__aeabi_d2iz, long int, double x) {
|
|
|
|
|
+ return x;
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_d2lz
|
|
#ifndef __aeabi_d2lz
|
|
|
-POLYFILL_FUN_1(__aeabi_d2lz, long, double)
|
|
|
|
|
|
|
+POLYFILL_FUN_1(__aeabi_d2lz, long, double x) {
|
|
|
|
|
+ return x;
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_d2uiz
|
|
#ifndef __aeabi_d2uiz
|
|
|
-POLYFILL_FUN_1(__aeabi_d2uiz, unsigned long int, double)
|
|
|
|
|
|
|
+POLYFILL_FUN_1(__aeabi_d2uiz, unsigned long int, double x) {
|
|
|
|
|
+ return x;
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_d2f
|
|
#ifndef __aeabi_d2f
|
|
|
-POLYFILL_FUN_1(__aeabi_d2f, float, double)
|
|
|
|
|
|
|
+POLYFILL_FUN_1(__aeabi_d2f, float, double x) {
|
|
|
|
|
+ return x;
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef __aeabi_ldivmod
|
|
#ifndef __aeabi_ldivmod
|
|
|
-POLYFILL_FUN_2(__aeabi_ldivmod, long, long, long)
|
|
|
|
|
|
|
+POLYFILL_FUN_2(__aeabi_ldivmod, long, long, long) {
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef strtox
|
|
#ifndef strtox
|
|
|
-POLYFILL_FUN_4(strtox, long long unsigned int, const char *, char **, int, long long unsigned int)
|
|
|
|
|
|
|
+POLYFILL_FUN_4(strtox, long long unsigned int, const char*, char**, int, long long unsigned int) {
|
|
|
|
|
+}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#if FLT_EVAL_METHOD == 0 || FLT_EVAL_METHOD == 1
|
|
#if FLT_EVAL_METHOD == 0 || FLT_EVAL_METHOD == 1
|
|
@@ -97,345 +123,276 @@ POLYFILL_FUN_4(strtox, long long unsigned int, const char *, char **, int, long
|
|
|
static const double_t toint = 1 / EPS;
|
|
static const double_t toint = 1 / EPS;
|
|
|
|
|
|
|
|
#ifndef FORCE_EVAL
|
|
#ifndef FORCE_EVAL
|
|
|
-#define FORCE_EVAL(x) \
|
|
|
|
|
- do \
|
|
|
|
|
- { \
|
|
|
|
|
- if (sizeof(x) == sizeof(float)) \
|
|
|
|
|
- { \
|
|
|
|
|
- volatile float __x; \
|
|
|
|
|
- __x = (x); \
|
|
|
|
|
- (void)__x; \
|
|
|
|
|
- } \
|
|
|
|
|
- else if (sizeof(x) == sizeof(double)) \
|
|
|
|
|
- { \
|
|
|
|
|
- volatile double __x; \
|
|
|
|
|
- __x = (x); \
|
|
|
|
|
- (void)__x; \
|
|
|
|
|
- } \
|
|
|
|
|
- else \
|
|
|
|
|
- { \
|
|
|
|
|
- volatile long double __x; \
|
|
|
|
|
- __x = (x); \
|
|
|
|
|
- (void)__x; \
|
|
|
|
|
- } \
|
|
|
|
|
- } while (0)
|
|
|
|
|
|
|
+#define FORCE_EVAL(x) \
|
|
|
|
|
+ do { \
|
|
|
|
|
+ if(sizeof(x) == sizeof(float)) { \
|
|
|
|
|
+ volatile float __x; \
|
|
|
|
|
+ __x = (x); \
|
|
|
|
|
+ (void)__x; \
|
|
|
|
|
+ } else if(sizeof(x) == sizeof(double)) { \
|
|
|
|
|
+ volatile double __x; \
|
|
|
|
|
+ __x = (x); \
|
|
|
|
|
+ (void)__x; \
|
|
|
|
|
+ } else { \
|
|
|
|
|
+ volatile long double __x; \
|
|
|
|
|
+ __x = (x); \
|
|
|
|
|
+ (void)__x; \
|
|
|
|
|
+ } \
|
|
|
|
|
+ } while(0)
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef floorf
|
|
#ifndef floorf
|
|
|
-float floorf(float x)
|
|
|
|
|
-{
|
|
|
|
|
- union
|
|
|
|
|
- {
|
|
|
|
|
- float f;
|
|
|
|
|
- uint32_t i;
|
|
|
|
|
- } u = {x};
|
|
|
|
|
- int e = (int)(u.i >> 23 & 0xff) - 0x7f;
|
|
|
|
|
- uint32_t m;
|
|
|
|
|
-
|
|
|
|
|
- if (e >= 23)
|
|
|
|
|
- return x;
|
|
|
|
|
- if (e >= 0)
|
|
|
|
|
- {
|
|
|
|
|
- m = 0x007fffff >> e;
|
|
|
|
|
- if ((u.i & m) == 0)
|
|
|
|
|
- return x;
|
|
|
|
|
- FORCE_EVAL(x + 0x1p120f);
|
|
|
|
|
- if (u.i >> 31)
|
|
|
|
|
- u.i += m;
|
|
|
|
|
- u.i &= ~m;
|
|
|
|
|
- }
|
|
|
|
|
- else
|
|
|
|
|
- {
|
|
|
|
|
- FORCE_EVAL(x + 0x1p120f);
|
|
|
|
|
- if (u.i >> 31 == 0)
|
|
|
|
|
- u.i = 0;
|
|
|
|
|
- else if (u.i << 1)
|
|
|
|
|
- u.f = -1.0;
|
|
|
|
|
- }
|
|
|
|
|
- return u.f;
|
|
|
|
|
|
|
+float floorf(float x) {
|
|
|
|
|
+ union {
|
|
|
|
|
+ float f;
|
|
|
|
|
+ uint32_t i;
|
|
|
|
|
+ } u = {x};
|
|
|
|
|
+ int e = (int)(u.i >> 23 & 0xff) - 0x7f;
|
|
|
|
|
+ uint32_t m;
|
|
|
|
|
+
|
|
|
|
|
+ if(e >= 23) return x;
|
|
|
|
|
+ if(e >= 0) {
|
|
|
|
|
+ m = 0x007fffff >> e;
|
|
|
|
|
+ if((u.i & m) == 0) return x;
|
|
|
|
|
+ FORCE_EVAL(x + 0x1p120f);
|
|
|
|
|
+ if(u.i >> 31) u.i += m;
|
|
|
|
|
+ u.i &= ~m;
|
|
|
|
|
+ } else {
|
|
|
|
|
+ FORCE_EVAL(x + 0x1p120f);
|
|
|
|
|
+ if(u.i >> 31 == 0)
|
|
|
|
|
+ u.i = 0;
|
|
|
|
|
+ else if(u.i << 1)
|
|
|
|
|
+ u.f = -1.0;
|
|
|
|
|
+ }
|
|
|
|
|
+ return u.f;
|
|
|
}
|
|
}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef floor
|
|
#ifndef floor
|
|
|
-double floor(double x)
|
|
|
|
|
-{
|
|
|
|
|
- union
|
|
|
|
|
- {
|
|
|
|
|
- double f;
|
|
|
|
|
- uint64_t i;
|
|
|
|
|
- } u = {x};
|
|
|
|
|
- int e = u.i >> 52 & 0x7ff;
|
|
|
|
|
- double_t y;
|
|
|
|
|
-
|
|
|
|
|
- if (e >= 0x3ff + 52 || x == 0)
|
|
|
|
|
- return x;
|
|
|
|
|
- /* y = int(x) - x, where int(x) is an integer neighbor of x */
|
|
|
|
|
- if (u.i >> 63)
|
|
|
|
|
- y = x - toint + toint - x;
|
|
|
|
|
- else
|
|
|
|
|
- y = x + toint - toint - x;
|
|
|
|
|
- /* special case because of non-nearest rounding modes */
|
|
|
|
|
- if (e <= 0x3ff - 1)
|
|
|
|
|
- {
|
|
|
|
|
- FORCE_EVAL(y);
|
|
|
|
|
- return u.i >> 63 ? -1 : 0;
|
|
|
|
|
- }
|
|
|
|
|
- if (y > 0)
|
|
|
|
|
- return x + y - 1;
|
|
|
|
|
- return x + y;
|
|
|
|
|
|
|
+double floor(double x) {
|
|
|
|
|
+ union {
|
|
|
|
|
+ double f;
|
|
|
|
|
+ uint64_t i;
|
|
|
|
|
+ } u = {x};
|
|
|
|
|
+ int e = u.i >> 52 & 0x7ff;
|
|
|
|
|
+ double_t y;
|
|
|
|
|
+
|
|
|
|
|
+ if(e >= 0x3ff + 52 || x == 0) return x;
|
|
|
|
|
+ /* y = int(x) - x, where int(x) is an integer neighbor of x */
|
|
|
|
|
+ if(u.i >> 63)
|
|
|
|
|
+ y = x - toint + toint - x;
|
|
|
|
|
+ else
|
|
|
|
|
+ y = x + toint - toint - x;
|
|
|
|
|
+ /* special case because of non-nearest rounding modes */
|
|
|
|
|
+ if(e <= 0x3ff - 1) {
|
|
|
|
|
+ FORCE_EVAL(y);
|
|
|
|
|
+ return u.i >> 63 ? -1 : 0;
|
|
|
|
|
+ }
|
|
|
|
|
+ if(y > 0) return x + y - 1;
|
|
|
|
|
+ return x + y;
|
|
|
}
|
|
}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef fmodf
|
|
#ifndef fmodf
|
|
|
-float fmodf(float x, float y)
|
|
|
|
|
-{
|
|
|
|
|
- union
|
|
|
|
|
- {
|
|
|
|
|
- float f;
|
|
|
|
|
|
|
+float fmodf(float x, float y) {
|
|
|
|
|
+ union {
|
|
|
|
|
+ float f;
|
|
|
|
|
+ uint32_t i;
|
|
|
|
|
+ } ux = {x}, uy = {y};
|
|
|
|
|
+ int ex = ux.i >> 23 & 0xff;
|
|
|
|
|
+ int ey = uy.i >> 23 & 0xff;
|
|
|
|
|
+ uint32_t sx = ux.i & 0x80000000;
|
|
|
uint32_t i;
|
|
uint32_t i;
|
|
|
- } ux = {x}, uy = {y};
|
|
|
|
|
- int ex = ux.i >> 23 & 0xff;
|
|
|
|
|
- int ey = uy.i >> 23 & 0xff;
|
|
|
|
|
- uint32_t sx = ux.i & 0x80000000;
|
|
|
|
|
- uint32_t i;
|
|
|
|
|
- uint32_t uxi = ux.i;
|
|
|
|
|
-
|
|
|
|
|
- if (uy.i << 1 == 0 || isnan(y) || ex == 0xff)
|
|
|
|
|
- return (x * y) / (x * y);
|
|
|
|
|
- if (uxi << 1 <= uy.i << 1)
|
|
|
|
|
- {
|
|
|
|
|
- if (uxi << 1 == uy.i << 1)
|
|
|
|
|
- return 0 * x;
|
|
|
|
|
- return x;
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- /* normalize x and y */
|
|
|
|
|
- if (!ex)
|
|
|
|
|
- {
|
|
|
|
|
- for (i = uxi << 9; i >> 31 == 0; ex--, i <<= 1)
|
|
|
|
|
- ;
|
|
|
|
|
- uxi <<= -ex + 1;
|
|
|
|
|
- }
|
|
|
|
|
- else
|
|
|
|
|
- {
|
|
|
|
|
- uxi &= -1U >> 9;
|
|
|
|
|
- uxi |= 1U << 23;
|
|
|
|
|
- }
|
|
|
|
|
- if (!ey)
|
|
|
|
|
- {
|
|
|
|
|
- for (i = uy.i << 9; i >> 31 == 0; ey--, i <<= 1)
|
|
|
|
|
- ;
|
|
|
|
|
- uy.i <<= -ey + 1;
|
|
|
|
|
- }
|
|
|
|
|
- else
|
|
|
|
|
- {
|
|
|
|
|
- uy.i &= -1U >> 9;
|
|
|
|
|
- uy.i |= 1U << 23;
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- /* x mod y */
|
|
|
|
|
- for (; ex > ey; ex--)
|
|
|
|
|
- {
|
|
|
|
|
|
|
+ uint32_t uxi = ux.i;
|
|
|
|
|
+
|
|
|
|
|
+ if(uy.i << 1 == 0 || isnan(y) || ex == 0xff) return (x * y) / (x * y);
|
|
|
|
|
+ if(uxi << 1 <= uy.i << 1) {
|
|
|
|
|
+ if(uxi << 1 == uy.i << 1) return 0 * x;
|
|
|
|
|
+ return x;
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ /* normalize x and y */
|
|
|
|
|
+ if(!ex) {
|
|
|
|
|
+ for(i = uxi << 9; i >> 31 == 0; ex--, i <<= 1)
|
|
|
|
|
+ ;
|
|
|
|
|
+ uxi <<= -ex + 1;
|
|
|
|
|
+ } else {
|
|
|
|
|
+ uxi &= -1U >> 9;
|
|
|
|
|
+ uxi |= 1U << 23;
|
|
|
|
|
+ }
|
|
|
|
|
+ if(!ey) {
|
|
|
|
|
+ for(i = uy.i << 9; i >> 31 == 0; ey--, i <<= 1)
|
|
|
|
|
+ ;
|
|
|
|
|
+ uy.i <<= -ey + 1;
|
|
|
|
|
+ } else {
|
|
|
|
|
+ uy.i &= -1U >> 9;
|
|
|
|
|
+ uy.i |= 1U << 23;
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ /* x mod y */
|
|
|
|
|
+ for(; ex > ey; ex--) {
|
|
|
|
|
+ i = uxi - uy.i;
|
|
|
|
|
+ if(i >> 31 == 0) {
|
|
|
|
|
+ if(i == 0) return 0 * x;
|
|
|
|
|
+ uxi = i;
|
|
|
|
|
+ }
|
|
|
|
|
+ uxi <<= 1;
|
|
|
|
|
+ }
|
|
|
i = uxi - uy.i;
|
|
i = uxi - uy.i;
|
|
|
- if (i >> 31 == 0)
|
|
|
|
|
- {
|
|
|
|
|
- if (i == 0)
|
|
|
|
|
- return 0 * x;
|
|
|
|
|
- uxi = i;
|
|
|
|
|
|
|
+ if(i >> 31 == 0) {
|
|
|
|
|
+ if(i == 0) return 0 * x;
|
|
|
|
|
+ uxi = i;
|
|
|
}
|
|
}
|
|
|
- uxi <<= 1;
|
|
|
|
|
- }
|
|
|
|
|
- i = uxi - uy.i;
|
|
|
|
|
- if (i >> 31 == 0)
|
|
|
|
|
- {
|
|
|
|
|
- if (i == 0)
|
|
|
|
|
- return 0 * x;
|
|
|
|
|
- uxi = i;
|
|
|
|
|
- }
|
|
|
|
|
- for (; uxi >> 23 == 0; uxi <<= 1, ex--)
|
|
|
|
|
- ;
|
|
|
|
|
-
|
|
|
|
|
- /* scale result up */
|
|
|
|
|
- if (ex > 0)
|
|
|
|
|
- {
|
|
|
|
|
- uxi -= 1U << 23;
|
|
|
|
|
- uxi |= (uint32_t)ex << 23;
|
|
|
|
|
- }
|
|
|
|
|
- else
|
|
|
|
|
- {
|
|
|
|
|
- uxi >>= -ex + 1;
|
|
|
|
|
- }
|
|
|
|
|
- uxi |= sx;
|
|
|
|
|
- ux.i = uxi;
|
|
|
|
|
- return ux.f;
|
|
|
|
|
|
|
+ for(; uxi >> 23 == 0; uxi <<= 1, ex--)
|
|
|
|
|
+ ;
|
|
|
|
|
+
|
|
|
|
|
+ /* scale result up */
|
|
|
|
|
+ if(ex > 0) {
|
|
|
|
|
+ uxi -= 1U << 23;
|
|
|
|
|
+ uxi |= (uint32_t)ex << 23;
|
|
|
|
|
+ } else {
|
|
|
|
|
+ uxi >>= -ex + 1;
|
|
|
|
|
+ }
|
|
|
|
|
+ uxi |= sx;
|
|
|
|
|
+ ux.i = uxi;
|
|
|
|
|
+ return ux.f;
|
|
|
}
|
|
}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef fmod
|
|
#ifndef fmod
|
|
|
-double fmod(double x, double y)
|
|
|
|
|
-{
|
|
|
|
|
- union
|
|
|
|
|
- {
|
|
|
|
|
- double f;
|
|
|
|
|
|
|
+double fmod(double x, double y) {
|
|
|
|
|
+ union {
|
|
|
|
|
+ double f;
|
|
|
|
|
+ uint64_t i;
|
|
|
|
|
+ } ux = {x}, uy = {y};
|
|
|
|
|
+ int ex = ux.i >> 52 & 0x7ff;
|
|
|
|
|
+ int ey = uy.i >> 52 & 0x7ff;
|
|
|
|
|
+ int sx = ux.i >> 63;
|
|
|
uint64_t i;
|
|
uint64_t i;
|
|
|
- } ux = {x}, uy = {y};
|
|
|
|
|
- int ex = ux.i >> 52 & 0x7ff;
|
|
|
|
|
- int ey = uy.i >> 52 & 0x7ff;
|
|
|
|
|
- int sx = ux.i >> 63;
|
|
|
|
|
- uint64_t i;
|
|
|
|
|
-
|
|
|
|
|
- /* in the followings uxi should be ux.i, but then gcc wrongly adds */
|
|
|
|
|
- /* float load/store to inner loops ruining performance and code size */
|
|
|
|
|
- uint64_t uxi = ux.i;
|
|
|
|
|
-
|
|
|
|
|
- if (uy.i << 1 == 0 || isnan(y) || ex == 0x7ff)
|
|
|
|
|
- return (x * y) / (x * y);
|
|
|
|
|
- if (uxi << 1 <= uy.i << 1)
|
|
|
|
|
- {
|
|
|
|
|
- if (uxi << 1 == uy.i << 1)
|
|
|
|
|
- return 0 * x;
|
|
|
|
|
- return x;
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- /* normalize x and y */
|
|
|
|
|
- if (!ex)
|
|
|
|
|
- {
|
|
|
|
|
- for (i = uxi << 12; i >> 63 == 0; ex--, i <<= 1)
|
|
|
|
|
- ;
|
|
|
|
|
- uxi <<= -ex + 1;
|
|
|
|
|
- }
|
|
|
|
|
- else
|
|
|
|
|
- {
|
|
|
|
|
- uxi &= -1ULL >> 12;
|
|
|
|
|
- uxi |= 1ULL << 52;
|
|
|
|
|
- }
|
|
|
|
|
- if (!ey)
|
|
|
|
|
- {
|
|
|
|
|
- for (i = uy.i << 12; i >> 63 == 0; ey--, i <<= 1)
|
|
|
|
|
- ;
|
|
|
|
|
- uy.i <<= -ey + 1;
|
|
|
|
|
- }
|
|
|
|
|
- else
|
|
|
|
|
- {
|
|
|
|
|
- uy.i &= -1ULL >> 12;
|
|
|
|
|
- uy.i |= 1ULL << 52;
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- /* x mod y */
|
|
|
|
|
- for (; ex > ey; ex--)
|
|
|
|
|
- {
|
|
|
|
|
|
|
+
|
|
|
|
|
+ /* in the followings uxi should be ux.i, but then gcc wrongly adds */
|
|
|
|
|
+ /* float load/store to inner loops ruining performance and code size */
|
|
|
|
|
+ uint64_t uxi = ux.i;
|
|
|
|
|
+
|
|
|
|
|
+ if(uy.i << 1 == 0 || isnan(y) || ex == 0x7ff) return (x * y) / (x * y);
|
|
|
|
|
+ if(uxi << 1 <= uy.i << 1) {
|
|
|
|
|
+ if(uxi << 1 == uy.i << 1) return 0 * x;
|
|
|
|
|
+ return x;
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ /* normalize x and y */
|
|
|
|
|
+ if(!ex) {
|
|
|
|
|
+ for(i = uxi << 12; i >> 63 == 0; ex--, i <<= 1)
|
|
|
|
|
+ ;
|
|
|
|
|
+ uxi <<= -ex + 1;
|
|
|
|
|
+ } else {
|
|
|
|
|
+ uxi &= -1ULL >> 12;
|
|
|
|
|
+ uxi |= 1ULL << 52;
|
|
|
|
|
+ }
|
|
|
|
|
+ if(!ey) {
|
|
|
|
|
+ for(i = uy.i << 12; i >> 63 == 0; ey--, i <<= 1)
|
|
|
|
|
+ ;
|
|
|
|
|
+ uy.i <<= -ey + 1;
|
|
|
|
|
+ } else {
|
|
|
|
|
+ uy.i &= -1ULL >> 12;
|
|
|
|
|
+ uy.i |= 1ULL << 52;
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ /* x mod y */
|
|
|
|
|
+ for(; ex > ey; ex--) {
|
|
|
|
|
+ i = uxi - uy.i;
|
|
|
|
|
+ if(i >> 63 == 0) {
|
|
|
|
|
+ if(i == 0) return 0 * x;
|
|
|
|
|
+ uxi = i;
|
|
|
|
|
+ }
|
|
|
|
|
+ uxi <<= 1;
|
|
|
|
|
+ }
|
|
|
i = uxi - uy.i;
|
|
i = uxi - uy.i;
|
|
|
- if (i >> 63 == 0)
|
|
|
|
|
- {
|
|
|
|
|
- if (i == 0)
|
|
|
|
|
- return 0 * x;
|
|
|
|
|
- uxi = i;
|
|
|
|
|
|
|
+ if(i >> 63 == 0) {
|
|
|
|
|
+ if(i == 0) return 0 * x;
|
|
|
|
|
+ uxi = i;
|
|
|
|
|
+ }
|
|
|
|
|
+ for(; uxi >> 52 == 0; uxi <<= 1, ex--)
|
|
|
|
|
+ ;
|
|
|
|
|
+
|
|
|
|
|
+ /* scale result */
|
|
|
|
|
+ if(ex > 0) {
|
|
|
|
|
+ uxi -= 1ULL << 52;
|
|
|
|
|
+ uxi |= (uint64_t)ex << 52;
|
|
|
|
|
+ } else {
|
|
|
|
|
+ uxi >>= -ex + 1;
|
|
|
}
|
|
}
|
|
|
- uxi <<= 1;
|
|
|
|
|
- }
|
|
|
|
|
- i = uxi - uy.i;
|
|
|
|
|
- if (i >> 63 == 0)
|
|
|
|
|
- {
|
|
|
|
|
- if (i == 0)
|
|
|
|
|
- return 0 * x;
|
|
|
|
|
- uxi = i;
|
|
|
|
|
- }
|
|
|
|
|
- for (; uxi >> 52 == 0; uxi <<= 1, ex--)
|
|
|
|
|
- ;
|
|
|
|
|
-
|
|
|
|
|
- /* scale result */
|
|
|
|
|
- if (ex > 0)
|
|
|
|
|
- {
|
|
|
|
|
- uxi -= 1ULL << 52;
|
|
|
|
|
- uxi |= (uint64_t)ex << 52;
|
|
|
|
|
- }
|
|
|
|
|
- else
|
|
|
|
|
- {
|
|
|
|
|
- uxi >>= -ex + 1;
|
|
|
|
|
- }
|
|
|
|
|
- uxi |= (uint64_t)sx << 63;
|
|
|
|
|
- ux.i = uxi;
|
|
|
|
|
- return ux.f;
|
|
|
|
|
|
|
+ uxi |= (uint64_t)sx << 63;
|
|
|
|
|
+ ux.i = uxi;
|
|
|
|
|
+ return ux.f;
|
|
|
}
|
|
}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef nearbyintf
|
|
#ifndef nearbyintf
|
|
|
-float nearbyintf(float x)
|
|
|
|
|
-{
|
|
|
|
|
- union
|
|
|
|
|
- {
|
|
|
|
|
- float f;
|
|
|
|
|
- uint32_t i;
|
|
|
|
|
- } u = {x};
|
|
|
|
|
- int e = u.i >> 23 & 0xff;
|
|
|
|
|
- int s = u.i >> 31;
|
|
|
|
|
- float_t y;
|
|
|
|
|
-
|
|
|
|
|
- if (e >= 0x7f + 23)
|
|
|
|
|
- return x;
|
|
|
|
|
- if (s)
|
|
|
|
|
- y = x - 0x1p23f + 0x1p23f;
|
|
|
|
|
- else
|
|
|
|
|
- y = x + 0x1p23f - 0x1p23f;
|
|
|
|
|
- if (y == 0)
|
|
|
|
|
- return s ? -0.0f : 0.0f;
|
|
|
|
|
- return y;
|
|
|
|
|
|
|
+float nearbyintf(float x) {
|
|
|
|
|
+ union {
|
|
|
|
|
+ float f;
|
|
|
|
|
+ uint32_t i;
|
|
|
|
|
+ } u = {x};
|
|
|
|
|
+ int e = u.i >> 23 & 0xff;
|
|
|
|
|
+ int s = u.i >> 31;
|
|
|
|
|
+ float_t y;
|
|
|
|
|
+
|
|
|
|
|
+ if(e >= 0x7f + 23) return x;
|
|
|
|
|
+ if(s)
|
|
|
|
|
+ y = x - 0x1p23f + 0x1p23f;
|
|
|
|
|
+ else
|
|
|
|
|
+ y = x + 0x1p23f - 0x1p23f;
|
|
|
|
|
+ if(y == 0) return s ? -0.0f : 0.0f;
|
|
|
|
|
+ return y;
|
|
|
}
|
|
}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef stroll
|
|
#ifndef stroll
|
|
|
-long long strtoll(const char *restrict s, char **restrict p, int base)
|
|
|
|
|
-{
|
|
|
|
|
- return strtox(s, p, base, LLONG_MIN);
|
|
|
|
|
|
|
+long long strtoll(const char* restrict s, char** restrict p, int base) {
|
|
|
|
|
+ return strtox(s, p, base, LLONG_MIN);
|
|
|
}
|
|
}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef pow
|
|
#ifndef pow
|
|
|
-double pow(double x, double y)
|
|
|
|
|
-{
|
|
|
|
|
- return powf(x, y);
|
|
|
|
|
|
|
+double pow(double x, double y) {
|
|
|
|
|
+ return powf(x, y);
|
|
|
}
|
|
}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef rint
|
|
#ifndef rint
|
|
|
|
|
|
|
|
-double rint(double x)
|
|
|
|
|
-{
|
|
|
|
|
- union
|
|
|
|
|
- {
|
|
|
|
|
- double f;
|
|
|
|
|
- uint64_t i;
|
|
|
|
|
- } u = {x};
|
|
|
|
|
- int e = u.i >> 52 & 0x7ff;
|
|
|
|
|
- int s = u.i >> 63;
|
|
|
|
|
- double_t y;
|
|
|
|
|
-
|
|
|
|
|
- if (e >= 0x3ff + 52)
|
|
|
|
|
- return x;
|
|
|
|
|
- if (s)
|
|
|
|
|
- y = x - toint + toint;
|
|
|
|
|
- else
|
|
|
|
|
- y = x + toint - toint;
|
|
|
|
|
- if (y == 0)
|
|
|
|
|
- return s ? -0.0 : 0;
|
|
|
|
|
- return y;
|
|
|
|
|
|
|
+double rint(double x) {
|
|
|
|
|
+ union {
|
|
|
|
|
+ double f;
|
|
|
|
|
+ uint64_t i;
|
|
|
|
|
+ } u = {x};
|
|
|
|
|
+ int e = u.i >> 52 & 0x7ff;
|
|
|
|
|
+ int s = u.i >> 63;
|
|
|
|
|
+ double_t y;
|
|
|
|
|
+
|
|
|
|
|
+ if(e >= 0x3ff + 52) return x;
|
|
|
|
|
+ if(s)
|
|
|
|
|
+ y = x - toint + toint;
|
|
|
|
|
+ else
|
|
|
|
|
+ y = x + toint - toint;
|
|
|
|
|
+ if(y == 0) return s ? -0.0 : 0;
|
|
|
|
|
+ return y;
|
|
|
}
|
|
}
|
|
|
#endif
|
|
#endif
|
|
|
|
|
|
|
|
#ifndef nearbyint
|
|
#ifndef nearbyint
|
|
|
-double nearbyint(double x)
|
|
|
|
|
-{
|
|
|
|
|
|
|
+double nearbyint(double x) {
|
|
|
#ifdef FE_INEXACT
|
|
#ifdef FE_INEXACT
|
|
|
#pragma STDC FENV_ACCESS ON
|
|
#pragma STDC FENV_ACCESS ON
|
|
|
- int e;
|
|
|
|
|
|
|
+ int e;
|
|
|
|
|
|
|
|
- e = fetestexcept(FE_INEXACT);
|
|
|
|
|
|
|
+ e = fetestexcept(FE_INEXACT);
|
|
|
#endif
|
|
#endif
|
|
|
- x = rint(x);
|
|
|
|
|
|
|
+ x = rint(x);
|
|
|
#ifdef FE_INEXACT
|
|
#ifdef FE_INEXACT
|
|
|
- if (!e)
|
|
|
|
|
- feclearexcept(FE_INEXACT);
|
|
|
|
|
|
|
+ if(!e) feclearexcept(FE_INEXACT);
|
|
|
#endif
|
|
#endif
|
|
|
- return x;
|
|
|
|
|
|
|
+ return x;
|
|
|
}
|
|
}
|
|
|
#endif
|
|
#endif
|