diff options
Diffstat (limited to 'src/math')
-rw-r--r-- | src/math/acoshl.c | 10 | ||||
-rw-r--r-- | src/math/fma.c | 2 | ||||
-rw-r--r-- | src/math/logf.c | 2 | ||||
-rw-r--r-- | src/math/powl.c | 34 | ||||
-rw-r--r-- | src/math/riscv32/copysign.c | 15 | ||||
-rw-r--r-- | src/math/riscv32/copysignf.c | 15 | ||||
-rw-r--r-- | src/math/riscv32/fabs.c | 15 | ||||
-rw-r--r-- | src/math/riscv32/fabsf.c | 15 | ||||
-rw-r--r-- | src/math/riscv32/fma.c | 15 | ||||
-rw-r--r-- | src/math/riscv32/fmaf.c | 15 | ||||
-rw-r--r-- | src/math/riscv32/fmax.c | 15 | ||||
-rw-r--r-- | src/math/riscv32/fmaxf.c | 15 | ||||
-rw-r--r-- | src/math/riscv32/fmin.c | 15 | ||||
-rw-r--r-- | src/math/riscv32/fminf.c | 15 | ||||
-rw-r--r-- | src/math/riscv32/sqrt.c | 15 | ||||
-rw-r--r-- | src/math/riscv32/sqrtf.c | 15 | ||||
-rw-r--r-- | src/math/sqrtl.c | 4 |
17 files changed, 212 insertions, 20 deletions
diff --git a/src/math/acoshl.c b/src/math/acoshl.c index 8d4b43f6..943cec17 100644 --- a/src/math/acoshl.c +++ b/src/math/acoshl.c @@ -10,14 +10,18 @@ long double acoshl(long double x) long double acoshl(long double x) { union ldshape u = {x}; - int e = u.i.se & 0x7fff; + int e = u.i.se; if (e < 0x3fff + 1) - /* |x| < 2, invalid if x < 1 or nan */ + /* 0 <= x < 2, invalid if x < 1 */ return log1pl(x-1 + sqrtl((x-1)*(x-1)+2*(x-1))); if (e < 0x3fff + 32) - /* |x| < 0x1p32 */ + /* 2 <= x < 0x1p32 */ return logl(2*x - 1/(x+sqrtl(x*x-1))); + if (e & 0x8000) + /* x < 0 or x = -0, invalid */ + return (x - x) / (x - x); + /* 0x1p32 <= x or nan */ return logl(x) + 0.693147180559945309417232121458176568L; } #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 diff --git a/src/math/fma.c b/src/math/fma.c index 0c6f90c9..adfadca8 100644 --- a/src/math/fma.c +++ b/src/math/fma.c @@ -53,7 +53,7 @@ double fma(double x, double y, double z) return x*y + z; if (nz.e >= ZEROINFNAN) { if (nz.e > ZEROINFNAN) /* z==0 */ - return x*y + z; + return x*y; return z; } diff --git a/src/math/logf.c b/src/math/logf.c index 7ee5d7fe..e4c2237c 100644 --- a/src/math/logf.c +++ b/src/math/logf.c @@ -53,7 +53,7 @@ float logf(float x) tmp = ix - OFF; i = (tmp >> (23 - LOGF_TABLE_BITS)) % N; k = (int32_t)tmp >> 23; /* arithmetic shift */ - iz = ix - (tmp & 0x1ff << 23); + iz = ix - (tmp & 0xff800000); invc = T[i].invc; logc = T[i].logc; z = (double_t)asfloat(iz); diff --git a/src/math/powl.c b/src/math/powl.c index 5b6da07b..6f64ea71 100644 --- a/src/math/powl.c +++ b/src/math/powl.c @@ -212,25 +212,33 @@ long double powl(long double x, long double y) } if (x == 1.0) return 1.0; /* 1**y = 1, even if y is nan */ - if (x == -1.0 && !isfinite(y)) - return 1.0; /* -1**inf = 1 */ if (y == 0.0) return 1.0; /* x**0 = 1, even if x is nan */ if (y == 1.0) return x; - if (y >= LDBL_MAX) { - if (x > 1.0 || x < -1.0) - return INFINITY; - if (x != 0.0) - return 0.0; - } - if (y <= -LDBL_MAX) { - if (x > 1.0 || x < -1.0) + /* if y*log2(x) < log2(LDBL_TRUE_MIN)-1 then x^y uflows to 0 + if y*log2(x) > -log2(LDBL_TRUE_MIN)+1 > LDBL_MAX_EXP then x^y oflows + if |x|!=1 then |log2(x)| > |log(x)| > LDBL_EPSILON/2 so + x^y oflows/uflows if |y|*LDBL_EPSILON/2 > -log2(LDBL_TRUE_MIN)+1 */ + if (fabsl(y) > 2*(-LDBL_MIN_EXP+LDBL_MANT_DIG+1)/LDBL_EPSILON) { + /* y is not an odd int */ + if (x == -1.0) + return 1.0; + if (y == INFINITY) { + if (x > 1.0 || x < -1.0) + return INFINITY; return 0.0; - if (x != 0.0 || y == -INFINITY) + } + if (y == -INFINITY) { + if (x > 1.0 || x < -1.0) + return 0.0; return INFINITY; + } + if ((x > 1.0 || x < -1.0) == (y > 0)) + return huge * huge; + return twom10000 * twom10000; } - if (x >= LDBL_MAX) { + if (x == INFINITY) { if (y > 0.0) return INFINITY; return 0.0; @@ -253,7 +261,7 @@ long double powl(long double x, long double y) yoddint = 1; } - if (x <= -LDBL_MAX) { + if (x == -INFINITY) { if (y > 0.0) { if (yoddint) return -INFINITY; diff --git a/src/math/riscv32/copysign.c b/src/math/riscv32/copysign.c new file mode 100644 index 00000000..c7854178 --- /dev/null +++ b/src/math/riscv32/copysign.c @@ -0,0 +1,15 @@ +#include <math.h> + +#if __riscv_flen >= 64 + +double copysign(double x, double y) +{ + __asm__ ("fsgnj.d %0, %1, %2" : "=f"(x) : "f"(x), "f"(y)); + return x; +} + +#else + +#include "../copysign.c" + +#endif diff --git a/src/math/riscv32/copysignf.c b/src/math/riscv32/copysignf.c new file mode 100644 index 00000000..a125611a --- /dev/null +++ b/src/math/riscv32/copysignf.c @@ -0,0 +1,15 @@ +#include <math.h> + +#if __riscv_flen >= 32 + +float copysignf(float x, float y) +{ + __asm__ ("fsgnj.s %0, %1, %2" : "=f"(x) : "f"(x), "f"(y)); + return x; +} + +#else + +#include "../copysignf.c" + +#endif diff --git a/src/math/riscv32/fabs.c b/src/math/riscv32/fabs.c new file mode 100644 index 00000000..5290b6f0 --- /dev/null +++ b/src/math/riscv32/fabs.c @@ -0,0 +1,15 @@ +#include <math.h> + +#if __riscv_flen >= 64 + +double fabs(double x) +{ + __asm__ ("fabs.d %0, %1" : "=f"(x) : "f"(x)); + return x; +} + +#else + +#include "../fabs.c" + +#endif diff --git a/src/math/riscv32/fabsf.c b/src/math/riscv32/fabsf.c new file mode 100644 index 00000000..f5032e35 --- /dev/null +++ b/src/math/riscv32/fabsf.c @@ -0,0 +1,15 @@ +#include <math.h> + +#if __riscv_flen >= 32 + +float fabsf(float x) +{ + __asm__ ("fabs.s %0, %1" : "=f"(x) : "f"(x)); + return x; +} + +#else + +#include "../fabsf.c" + +#endif diff --git a/src/math/riscv32/fma.c b/src/math/riscv32/fma.c new file mode 100644 index 00000000..99b05713 --- /dev/null +++ b/src/math/riscv32/fma.c @@ -0,0 +1,15 @@ +#include <math.h> + +#if __riscv_flen >= 64 + +double fma(double x, double y, double z) +{ + __asm__ ("fmadd.d %0, %1, %2, %3" : "=f"(x) : "f"(x), "f"(y), "f"(z)); + return x; +} + +#else + +#include "../fma.c" + +#endif diff --git a/src/math/riscv32/fmaf.c b/src/math/riscv32/fmaf.c new file mode 100644 index 00000000..f9dc47ed --- /dev/null +++ b/src/math/riscv32/fmaf.c @@ -0,0 +1,15 @@ +#include <math.h> + +#if __riscv_flen >= 32 + +float fmaf(float x, float y, float z) +{ + __asm__ ("fmadd.s %0, %1, %2, %3" : "=f"(x) : "f"(x), "f"(y), "f"(z)); + return x; +} + +#else + +#include "../fmaf.c" + +#endif diff --git a/src/math/riscv32/fmax.c b/src/math/riscv32/fmax.c new file mode 100644 index 00000000..023709cd --- /dev/null +++ b/src/math/riscv32/fmax.c @@ -0,0 +1,15 @@ +#include <math.h> + +#if __riscv_flen >= 64 + +double fmax(double x, double y) +{ + __asm__ ("fmax.d %0, %1, %2" : "=f"(x) : "f"(x), "f"(y)); + return x; +} + +#else + +#include "../fmax.c" + +#endif diff --git a/src/math/riscv32/fmaxf.c b/src/math/riscv32/fmaxf.c new file mode 100644 index 00000000..863d2bd1 --- /dev/null +++ b/src/math/riscv32/fmaxf.c @@ -0,0 +1,15 @@ +#include <math.h> + +#if __riscv_flen >= 32 + +float fmaxf(float x, float y) +{ + __asm__ ("fmax.s %0, %1, %2" : "=f"(x) : "f"(x), "f"(y)); + return x; +} + +#else + +#include "../fmaxf.c" + +#endif diff --git a/src/math/riscv32/fmin.c b/src/math/riscv32/fmin.c new file mode 100644 index 00000000..a4e3b067 --- /dev/null +++ b/src/math/riscv32/fmin.c @@ -0,0 +1,15 @@ +#include <math.h> + +#if __riscv_flen >= 64 + +double fmin(double x, double y) +{ + __asm__ ("fmin.d %0, %1, %2" : "=f"(x) : "f"(x), "f"(y)); + return x; +} + +#else + +#include "../fmin.c" + +#endif diff --git a/src/math/riscv32/fminf.c b/src/math/riscv32/fminf.c new file mode 100644 index 00000000..32156e80 --- /dev/null +++ b/src/math/riscv32/fminf.c @@ -0,0 +1,15 @@ +#include <math.h> + +#if __riscv_flen >= 32 + +float fminf(float x, float y) +{ + __asm__ ("fmin.s %0, %1, %2" : "=f"(x) : "f"(x), "f"(y)); + return x; +} + +#else + +#include "../fminf.c" + +#endif diff --git a/src/math/riscv32/sqrt.c b/src/math/riscv32/sqrt.c new file mode 100644 index 00000000..867a504c --- /dev/null +++ b/src/math/riscv32/sqrt.c @@ -0,0 +1,15 @@ +#include <math.h> + +#if __riscv_flen >= 64 + +double sqrt(double x) +{ + __asm__ ("fsqrt.d %0, %1" : "=f"(x) : "f"(x)); + return x; +} + +#else + +#include "../sqrt.c" + +#endif diff --git a/src/math/riscv32/sqrtf.c b/src/math/riscv32/sqrtf.c new file mode 100644 index 00000000..610c2cf8 --- /dev/null +++ b/src/math/riscv32/sqrtf.c @@ -0,0 +1,15 @@ +#include <math.h> + +#if __riscv_flen >= 32 + +float sqrtf(float x) +{ + __asm__ ("fsqrt.s %0, %1" : "=f"(x) : "f"(x)); + return x; +} + +#else + +#include "../sqrtf.c" + +#endif diff --git a/src/math/sqrtl.c b/src/math/sqrtl.c index 1b9f19c7..a231b3f2 100644 --- a/src/math/sqrtl.c +++ b/src/math/sqrtl.c @@ -205,7 +205,7 @@ long double sqrtl(long double x) top = (top + 0x3fff) >> 1; /* r ~ 1/sqrt(m) */ - static const uint64_t three = 0xc0000000; + const uint64_t three = 0xc0000000; uint64_t r, s, d, u, i; i = (ix.hi >> 42) % 128; r = (uint32_t)__rsqrt_tab[i] << 16; @@ -227,7 +227,7 @@ long double sqrtl(long double x) r = mul64(u, r) << 1; /* |r sqrt(m) - 1| < 0x1.c001p-59, switch to 128bit */ - static const u128 threel = {.hi=three<<32, .lo=0}; + const u128 threel = {.hi=three<<32, .lo=0}; u128 rl, sl, dl, ul; rl.hi = r; rl.lo = 0; |