libm/s_nextafterf.c: forgot to "svn add" it

diff --git a/libm/s_nextafterf.c b/libm/s_nextafterf.c
new file mode 100644 (file)
index 0000000..8dee00f
--- /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