summaryrefslogtreecommitdiff
path: root/src/math
diff options
context:
space:
mode:
Diffstat (limited to 'src/math')
-rw-r--r--src/math/acoshf.c8
-rw-r--r--src/math/acoshl.c10
-rw-r--r--src/math/arm/fabs.c2
-rw-r--r--src/math/arm/sqrt.c2
-rw-r--r--src/math/expm1f.c3
-rw-r--r--src/math/fma.c2
-rw-r--r--src/math/fmaf.c21
-rw-r--r--src/math/logf.c2
-rw-r--r--src/math/powerpc/fabs.c2
-rw-r--r--src/math/powerpc/fabsf.c2
-rw-r--r--src/math/powerpc/fma.c2
-rw-r--r--src/math/powerpc/fmaf.c2
-rw-r--r--src/math/powl.c34
-rw-r--r--src/math/riscv32/copysign.c15
-rw-r--r--src/math/riscv32/copysignf.c15
-rw-r--r--src/math/riscv32/fabs.c15
-rw-r--r--src/math/riscv32/fabsf.c15
-rw-r--r--src/math/riscv32/fma.c15
-rw-r--r--src/math/riscv32/fmaf.c15
-rw-r--r--src/math/riscv32/fmax.c15
-rw-r--r--src/math/riscv32/fmaxf.c15
-rw-r--r--src/math/riscv32/fmin.c15
-rw-r--r--src/math/riscv32/fminf.c15
-rw-r--r--src/math/riscv32/sqrt.c15
-rw-r--r--src/math/riscv32/sqrtf.c15
-rw-r--r--src/math/sqrtl.c4
26 files changed, 233 insertions, 43 deletions
diff --git a/src/math/acoshf.c b/src/math/acoshf.c
index 8a4ec4d5..b773d48e 100644
--- a/src/math/acoshf.c
+++ b/src/math/acoshf.c
@@ -15,12 +15,12 @@ float acoshf(float x)
uint32_t a = u.i & 0x7fffffff;
if (a < 0x3f800000+(1<<23))
- /* |x| < 2, invalid if x < 1 or nan */
+ /* |x| < 2, invalid if x < 1 */
/* up to 2ulp error in [1,1.125] */
return log1pf(x-1 + sqrtf((x-1)*(x-1)+2*(x-1)));
- if (a < 0x3f800000+(12<<23))
- /* |x| < 0x1p12 */
+ if (u.i < 0x3f800000+(12<<23))
+ /* 2 <= x < 0x1p12 */
return logf(2*x - 1/(x+sqrtf(x*x-1)));
- /* x >= 0x1p12 */
+ /* x >= 0x1p12 or x <= -2 or nan */
return logf(x) + 0.693147180559945309417232121458176568f;
}
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/arm/fabs.c b/src/math/arm/fabs.c
index f890520a..6e1d367d 100644
--- a/src/math/arm/fabs.c
+++ b/src/math/arm/fabs.c
@@ -1,6 +1,6 @@
#include <math.h>
-#if __ARM_PCS_VFP
+#if __ARM_PCS_VFP && __ARM_FP&8
double fabs(double x)
{
diff --git a/src/math/arm/sqrt.c b/src/math/arm/sqrt.c
index 874af960..567e2e91 100644
--- a/src/math/arm/sqrt.c
+++ b/src/math/arm/sqrt.c
@@ -1,6 +1,6 @@
#include <math.h>
-#if __ARM_PCS_VFP || (__VFP_FP__ && !__SOFTFP__)
+#if (__ARM_PCS_VFP || (__VFP_FP__ && !__SOFTFP__)) && (__ARM_FP&8)
double sqrt(double x)
{
diff --git a/src/math/expm1f.c b/src/math/expm1f.c
index 297e0b44..09a41afe 100644
--- a/src/math/expm1f.c
+++ b/src/math/expm1f.c
@@ -16,7 +16,6 @@
#include "libm.h"
static const float
-o_threshold = 8.8721679688e+01, /* 0x42b17180 */
ln2_hi = 6.9313812256e-01, /* 0x3f317180 */
ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */
invln2 = 1.4426950216e+00, /* 0x3fb8aa3b */
@@ -41,7 +40,7 @@ float expm1f(float x)
return x;
if (sign)
return -1;
- if (x > o_threshold) {
+ if (hx > 0x42b17217) { /* x > log(FLT_MAX) */
x *= 0x1p127f;
return x;
}
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/fmaf.c b/src/math/fmaf.c
index 80f5cd8a..7c65acf1 100644
--- a/src/math/fmaf.c
+++ b/src/math/fmaf.c
@@ -77,17 +77,16 @@ float fmaf(float x, float y, float z)
* If result is inexact, and exactly halfway between two float values,
* we need to adjust the low-order bit in the direction of the error.
*/
-#ifdef FE_TOWARDZERO
- fesetround(FE_TOWARDZERO);
-#endif
- volatile double vxy = xy; /* XXX work around gcc CSE bug */
- double adjusted_result = vxy + z;
- fesetround(FE_TONEAREST);
- if (result == adjusted_result) {
- u.f = adjusted_result;
+ double err;
+ int neg = u.i >> 63;
+ if (neg == (z > xy))
+ err = xy - result + z;
+ else
+ err = z - result + xy;
+ if (neg == (err < 0))
u.i++;
- adjusted_result = u.f;
- }
- z = adjusted_result;
+ else
+ u.i--;
+ z = u.f;
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/powerpc/fabs.c b/src/math/powerpc/fabs.c
index 0efc21ef..9453a3aa 100644
--- a/src/math/powerpc/fabs.c
+++ b/src/math/powerpc/fabs.c
@@ -1,6 +1,6 @@
#include <math.h>
-#if defined(_SOFT_FLOAT) || defined(BROKEN_PPC_D_ASM)
+#if defined(_SOFT_FLOAT) || defined(__NO_FPRS__) || defined(BROKEN_PPC_D_ASM)
#include "../fabs.c"
diff --git a/src/math/powerpc/fabsf.c b/src/math/powerpc/fabsf.c
index d88b5911..2e9da588 100644
--- a/src/math/powerpc/fabsf.c
+++ b/src/math/powerpc/fabsf.c
@@ -1,6 +1,6 @@
#include <math.h>
-#ifdef _SOFT_FLOAT
+#if defined(_SOFT_FLOAT) || defined(__NO_FPRS__)
#include "../fabsf.c"
diff --git a/src/math/powerpc/fma.c b/src/math/powerpc/fma.c
index 135c9903..0eb2ba1e 100644
--- a/src/math/powerpc/fma.c
+++ b/src/math/powerpc/fma.c
@@ -1,6 +1,6 @@
#include <math.h>
-#if defined(_SOFT_FLOAT) || defined(BROKEN_PPC_D_ASM)
+#if defined(_SOFT_FLOAT) || defined(__NO_FPRS__) || defined(BROKEN_PPC_D_ASM)
#include "../fma.c"
diff --git a/src/math/powerpc/fmaf.c b/src/math/powerpc/fmaf.c
index a99a2a3b..dc1a749d 100644
--- a/src/math/powerpc/fmaf.c
+++ b/src/math/powerpc/fmaf.c
@@ -1,6 +1,6 @@
#include <math.h>
-#ifdef _SOFT_FLOAT
+#if defined(_SOFT_FLOAT) || defined(__NO_FPRS__)
#include "../fmaf.c"
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;
generated by cgit v1.2.1 (git 2.18.0) at 2024-05-02 05:37:20 +0000