summaryrefslogtreecommitdiff
path: root/libm
diff options
context:
space:
mode:
Diffstat (limited to 'libm')
-rw-r--r--libm/s_nextafterf.c103
1 files changed, 103 insertions, 0 deletions
diff --git a/libm/s_nextafterf.c b/libm/s_nextafterf.c
new file mode 100644
index 000000000..8dee00ff7
--- /dev/null
+++ b/libm/s_nextafterf.c
@@ -0,0 +1,103 @@
+/* s_nextafterf.c -- float version of s_nextafter.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifndef math_opt_barrier
+# define math_opt_barrier(x) ({ __typeof (x) __x = x; __asm ("" : "+m" (__x)); __x; })
+# define math_force_eval(x) __asm __volatile ("" : : "m" (x))
+#endif
+
+float nextafterf(float x, float y)
+{
+ int32_t hx, hy, ix, iy;
+
+ GET_FLOAT_WORD(hx, x);
+ GET_FLOAT_WORD(hy, y);
+ ix = hx & 0x7fffffff; /* |x| */
+ iy = hy & 0x7fffffff; /* |y| */
+
+ /* x is nan or y is nan? */
+ if ((ix > 0x7f800000) || (iy > 0x7f800000))
+ return x + y;
+
+ if (x == y)
+ return y;
+
+ if (ix == 0) { /* x == 0? */
+ float u;
+ /* return +-minsubnormal */
+ SET_FLOAT_WORD(x, (hy & 0x80000000) | 1);
+ u = math_opt_barrier(x);
+ u = u * u;
+ math_force_eval(u); /* raise underflow flag */
+ return x;
+ }
+
+ if (hx >= 0) { /* x > 0 */
+ if (hx > hy) { /* x > y: x -= ulp */
+ hx -= 1;
+ } else { /* x < y: x += ulp */
+ hx += 1;
+ }
+ } else { /* x < 0 */
+ if (hy >= 0 || hx > hy) { /* x < y: x -= ulp */
+ hx -= 1;
+ } else { /* x > y: x += ulp */
+ hx += 1;
+ }
+ }
+ hy = hx & 0x7f800000;
+ if (hy >= 0x7f800000) {
+ x = x + x; /* overflow */
+//?? if (FLT_EVAL_METHOD != 0)
+// asm ("" : "+m"(x));
+ return x; /* overflow */
+ }
+ if (hy < 0x00800000) {
+ float u = x * x; /* underflow */
+ math_force_eval(u); /* raise underflow flag */
+ }
+ SET_FLOAT_WORD(x, hx);
+ return x;
+}
+
+#if 0
+/* "testprog N a b"
+ * calculates a = nextafterf(a, b) and prints a as float
+ * and as raw bytes; repeats it N times.
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+int main(int argc, char **argv)
+{
+ int cnt, i;
+ float a, b;
+ cnt = atoi(argv[1]);
+ a = strtod(argv[2], NULL);
+ b = strtod(argv[3], NULL);
+ while (cnt-- > 0) {
+ for (i = 0; i < sizeof(a); i++) {
+ unsigned char c = ((char*)(&a))[i];
+ printf("%x%x", (c >> 4), (c & 0xf));
+ }
+ printf(" %f\n", a);
+ a = nextafterf(a, b);
+ }
+ return 0;
+}
+#endif