Update.

1999-12-29  Ulrich Drepper  <drepper@cygnus.com>

	* soft-fp/*: Tons of new files to implement floating-point arithmetic
	in software.
	Contributed by  Richard Henderson, Jakub Jelinek and others.

1 files changed, 788 insertions, 0 deletions

diff --git a/soft-fp/testit.c b/soft-fp/testit.c
new file mode 100644
index 0000000..46a50c4
--- /dev/null
+++ b/soft-fp/testit.c
@@ -0,0 +1,788 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "soft-fp.h"
+#include "single.h"
+#include "double.h"
+
+#include <fpu_control.h>
+
+/*======================================================================*/
+/* declarations for the functions we are testing */
+
+double __adddf3(double, double);
+double __subdf3(double, double);
+double __muldf3(double, double);
+double __divdf3(double, double);
+double __negdf2(double);
+double __sqrtdf2(double);
+double __negdf3(double a, double dummy) { return __negdf2(a); }
+double __sqrtdf3(double a, double dummy) { return __sqrtdf2(a); }
+
+float __addsf3(float, float);
+float __subsf3(float, float);
+float __mulsf3(float, float);
+float __divsf3(float, float);
+float __negsf2(float);
+float __sqrtsf2(float);
+float __negsf3(float a, float dummy) { return __negsf2(a); }
+float __sqrtsf3(float a, float dummy) { return __sqrtsf2(a); }
+
+int __fixdfsi(double);
+int __fixsfsi(float);
+double __floatsidf(int);
+float __floatsisf(int);
+double __extendsfdf2(float);
+float __truncdfsf2(double);
+
+int __eqdf2(double, double);
+int __nedf2(double, double);
+int __gtdf2(double, double);
+int __gedf2(double, double);
+int __ltdf2(double, double);
+int __ledf2(double, double);
+
+int __eqsf2(float, float);
+int __nesf2(float, float);
+int __gtsf2(float, float);
+int __gesf2(float, float);
+int __ltsf2(float, float);
+int __lesf2(float, float);
+
+/*======================================================================*/
+/* definitions for functions we are checking against */
+
+double r_adddf3(double a, double b) { return a + b; }
+double r_subdf3(double a, double b) { return a - b; }
+double r_muldf3(double a, double b) { return a * b; }
+double r_divdf3(double a, double b) { return a / b; }
+double r_negdf3(double a, double b) { return -a; }
+double sqrt(double x);
+double r_sqrtdf3(double a, double b) { return sqrt(a); }
+
+float r_addsf3(float a, float b) { return a + b; }
+float r_subsf3(float a, float b) { return a - b; }
+float r_mulsf3(float a, float b) { return a * b; }
+float r_divsf3(float a, float b) { return a / b; }
+float r_negsf3(float a, float b) { return -a; }
+float sqrtf(float x);
+float r_sqrtsf3(float a, float b) { return sqrtf(a); }
+
+int r_fixdfsi(double a) { return (int)a; }
+int r_fixsfsi(float a) { return (int)a; }
+double r_floatsidf(int a) { return (double)a; }
+float r_floatsisf(int a) { return (float)a; }
+double r_extendsfdf2(float a) { return (double)a; }
+float r_truncdfsf2(double a) { return (float)a; }
+
+int r_eqdf2(double a, double b) { return !(a == b); }
+int r_nedf2(double a, double b) { return a != b; }
+int r_gtdf2(double a, double b) { return a > b; }
+int r_gedf2(double a, double b) { return (a >= b) - 1; }
+int r_ltdf2(double a, double b) { return -(a < b); }
+int r_ledf2(double a, double b) { return 1 - (a <= b); }
+
+int r_eqsf2(float a, float b) { return !(a == b); }
+int r_nesf2(float a, float b) { return a != b; }
+int r_gtsf2(float a, float b) { return a > b; }
+int r_gesf2(float a, float b) { return (a >= b) - 1; }
+int r_ltsf2(float a, float b) { return -(a < b); }
+int r_lesf2(float a, float b) { return 1 - (a <= b); }
+
+/*======================================================================*/
+
+void print_float(float x)
+{
+    union _FP_UNION_S ux;
+    ux.flt = x;
+    printf("%-20.8e %X %02X %06lX",
+	   x, ux.bits.sign, ux.bits.exp, (unsigned long)ux.bits.frac);
+}
+
+void print_double(double x)
+{
+    union _FP_UNION_D ux;
+    ux.flt = x;
+#if _FP_W_TYPE_SIZE < _FP_FRACBITS_D
+    printf("%-30.18e %X %04X %06lX%08lX",
+	   x, ux.bits.sign, ux.bits.exp,
+	   (unsigned long)ux.bits.frac1, (unsigned long)ux.bits.frac0);
+#else
+    printf("%-30.18e %X %04X %014lX",
+	   x, ux.bits.sign, ux.bits.exp,
+	   (unsigned long)ux.bits.frac);
+#endif
+}
+
+float rand_float(void)
+{
+    union {
+	union _FP_UNION_S u;
+	int i;
+    } u;
+
+    u.i = lrand48() << 1;
+
+    if (u.u.bits.exp == _FP_EXPMAX_S)
+	u.u.bits.exp--;
+    else if (u.u.bits.exp == 0 && u.u.bits.frac != 0)
+	u.u.bits.exp++;
+
+    return u.u.flt;
+}
+
+
+double rand_double(void)
+{
+    union {
+	union _FP_UNION_D u;
+	int i[2];
+    } u;
+
+    u.i[0] = lrand48() << 1;
+    u.i[1] = lrand48() << 1;
+
+    if (u.u.bits.exp == _FP_EXPMAX_D)
+	u.u.bits.exp--;
+#if _FP_W_TYPE_SIZE < _FP_FRACBITS_D
+    else if (u.u.bits.exp == 0 && !(u.u.bits.frac0 == 0 && u.u.bits.frac1 == 0))
+	u.u.bits.exp++;
+#else
+    else if (u.u.bits.exp == 0 && u.u.bits.frac != 0)
+	u.u.bits.exp++;
+#endif
+
+    return u.u.flt;
+}
+
+#define NSPECIALS  10
+
+float gen_special_float(int i)
+{
+    FP_DECL_EX;
+    FP_DECL_S(X);
+    float x;
+
+    switch (i & ~1)
+    {
+      case 0:
+	X_c = FP_CLS_NAN; X_f = 0x1234;
+	break;
+      case 2:
+	X_c = FP_CLS_NAN; X_f = 0x1;
+	break;
+      case 4:
+	X_c = FP_CLS_INF;
+	break;
+      case 6:
+	X_c = FP_CLS_ZERO;
+	break;
+      case 8:
+	X_c = FP_CLS_NORMAL; X_e = 0;
+	X_f = 0x4321;
+	break;
+    }
+    X_s = (i & 1);
+
+    FP_PACK_S(x, X);
+    return x;
+}
+
+double gen_special_double(int i)
+{
+    FP_DECL_EX;
+    FP_DECL_D(X);
+    double x;
+
+    switch (i & ~1)
+    {
+      case 0:
+	X_c = FP_CLS_NAN;
+#if _FP_W_TYPE_SIZE < _FP_FRACBITS_D
+	__FP_FRAC_SET_2(X, _FP_QNANBIT_D, 0x1234);
+#else
+	_FP_FRAC_SET_1(X, _FP_QNANBIT_D | 0x1234);
+#endif
+	break;
+      case 2:
+	X_c = FP_CLS_NAN;
+#if _FP_W_TYPE_SIZE < _FP_FRACBITS_D
+	__FP_FRAC_SET_2(X, _FP_QNANBIT_D, 0x1);
+#else
+	_FP_FRAC_SET_1(X, _FP_QNANBIT_D | 0x1);
+#endif
+	break;
+      case 4:
+	X_c = FP_CLS_INF;
+	break;
+      case 6:
+	X_c = FP_CLS_ZERO;
+	break;
+      case 8:
+	X_c = FP_CLS_NORMAL; X_e = 0;
+#if _FP_W_TYPE_SIZE < _FP_FRACBITS_D
+	__FP_FRAC_SET_2(X, 0, 0x87654321);
+#else
+	_FP_FRAC_SET_1(X, 0x87654321);
+#endif
+	break;
+    }
+    X_s = (i & 1);
+
+    FP_PACK_D(x, X);
+    return x;
+}
+
+float build_float(const char *s, const char *e, const char *f)
+{
+    union _FP_UNION_S u;
+
+    u.bits.sign = strtoul(s, 0, 16);
+    u.bits.exp = strtoul(e, 0, 16);
+    u.bits.frac = strtoul(f, 0, 16);
+
+    return u.flt;
+}
+      
+double build_double(const char *s, const char *e, const char *f)
+{
+    union _FP_UNION_D u;
+
+    u.bits.sign = strtoul(s, 0, 16);
+    u.bits.exp = strtoul(e, 0, 16);
+#if _FP_W_TYPE_SIZE < _FP_FRACBITS_D
+    {
+	size_t len = strlen(f)+1;
+	char *dup = memcpy(alloca(len), f, len);
+	char *low = dup + len - _FP_W_TYPE_SIZE/4 - 1;
+
+	u.bits.frac0 = strtoul(low, 0, 16);
+	*low = 0;
+	u.bits.frac1 = strtoul(dup, 0, 16);
+    }
+#else
+    u.bits.frac = strtoul(f, 0, 16);
+#endif
+	
+    return u.flt;
+}
+      
+/*======================================================================*/
+
+fpu_control_t fcw0, fcw1;
+
+void test_float_arith(float (*tf)(float, float),
+		      float (*rf)(float, float),
+		      float x, float y)
+{
+    float tr, rr;
+    rr = (*rf)(x, y);
+    tr = (*tf)(x, y);
+    if (memcmp(&tr, &rr, sizeof(float)) != 0)
+    {
+	fputs("error:\n\tx     = ", stdout); print_float(x);
+	fputs("\n\ty     = ", stdout); print_float(y);
+	fputs("\n\ttrue  = ", stdout); print_float(rr);
+	fputs("\n\tfalse = ", stdout); print_float(tr);
+	putchar('\n');
+    }
+}
+
+void test_double_arith(double (*tf)(double, double),
+		       double (*rf)(double, double),
+		       double x, double y)
+{
+    double tr, rr;
+#ifdef __i386__
+    /* Don't worry.  Even this does not make it error free
+       on ia32.  If the result is denormal,  it will not
+       honour the double precision and generate bad results
+       anyway.  On the other side,  who wants to use ia32
+       for IEEE math?  I don't.  */
+    _FPU_GETCW(fcw0);
+    fcw1 = ((fcw0 & ~_FPU_EXTENDED) | _FPU_DOUBLE);
+    _FPU_SETCW(fcw1);
+#endif    
+    rr = (*rf)(x, y);
+#ifdef __i386__
+    _FPU_SETCW(fcw0);
+#endif
+    tr = (*tf)(x, y);
+    if (memcmp(&tr, &rr, sizeof(double)) != 0)
+    {
+	fputs("error:\n\tx     = ", stdout); print_double(x);
+	fputs("\n\ty     = ", stdout); print_double(y);
+	fputs("\n\ttrue  = ", stdout); print_double(rr);
+	fputs("\n\tfalse = ", stdout); print_double(tr);
+	putchar('\n');
+    }
+}
+
+void test_float_double_conv(float x)
+{
+    double tr, rr;
+    rr = r_extendsfdf2(x);
+    tr = __extendsfdf2(x);
+    if (memcmp(&tr, &rr, sizeof(double)) != 0)
+    {
+	fputs("error:\n\tx     = ", stdout); print_float(x);
+	fputs("\n\ttrue  = ", stdout); print_double(rr);
+	fputs("\n\tfalse = ", stdout); print_double(tr);
+	putchar('\n');
+    }
+}
+
+void test_double_float_conv(double x)
+{
+    float tr, rr;
+    rr = r_truncdfsf2(x);
+    tr = __truncdfsf2(x);
+    if (memcmp(&tr, &rr, sizeof(float)) != 0)
+    {
+	fputs("error:\n\tx     = ", stdout); print_double(x);
+	fputs("\n\ttrue  = ", stdout); print_float(rr);
+	fputs("\n\tfalse = ", stdout); print_float(tr);
+	putchar('\n');
+    }
+}
+
+void test_int_float_conv(int x)
+{
+    float tr, rr;
+    rr = r_floatsisf(x);
+    tr = __floatsisf(x);
+    if (memcmp(&tr, &rr, sizeof(float)) != 0)
+    {
+	printf("error\n\tx     = %d", x);
+	fputs("\n\ttrue  = ", stdout); print_float(rr);
+	fputs("\n\tfalse = ", stdout); print_float(tr);
+	putchar('\n');
+    }
+}
+
+void test_int_double_conv(int x)
+{
+    double tr, rr;
+    rr = r_floatsidf(x);
+    tr = __floatsidf(x);
+    if (memcmp(&tr, &rr, sizeof(double)) != 0)
+    {
+	printf("error\n\tx     = %d", x);
+	fputs("\n\ttrue  = ", stdout); print_double(rr);
+	fputs("\n\tfalse = ", stdout); print_double(tr);
+	putchar('\n');
+    }
+}
+
+void test_float_int_conv(float x)
+{
+    int tr, rr;
+    rr = r_fixsfsi(x);
+    tr = __fixsfsi(x);
+    if (rr != tr)
+    {
+	fputs("error:\n\tx     = ", stdout); print_float(x);
+	printf("\n\ttrue  = %d\n\tfalse = %d\n", rr, tr);
+    }
+}
+	
+void test_double_int_conv(double x)
+{
+    int tr, rr;
+    rr = r_fixsfsi(x);
+    tr = __fixsfsi(x);
+    if (rr != tr)
+    {
+	fputs("error:\n\tx     = ", stdout); print_double(x);
+	printf("\n\ttrue  = %d\n\tfalse = %d\n", rr, tr);
+    }
+}
+
+int eq0(int x) { return x == 0; }
+int ne0(int x) { return x != 0; }
+int le0(int x) { return x <= 0; }
+int lt0(int x) { return x < 0; }
+int ge0(int x) { return x >= 0; }
+int gt0(int x) { return x > 0; }
+
+void test_float_cmp(int (*tf)(float, float),
+		    int (*rf)(float, float),
+		    int (*cmp0)(int),
+		    float x, float y)
+{
+    int tr, rr;
+    rr = (*rf)(x, y);
+    tr = (*tf)(x, y);
+    if (cmp0(rr) != cmp0(tr))
+    {
+	fputs("error:\n\tx     = ", stdout); print_float(x);
+	fputs("\n\ty     = ", stdout); print_float(y);
+	printf("\n\ttrue  = %d\n\tfalse = %d\n", rr, tr);
+    }
+}
+
+void test_double_cmp(int (*tf)(double, double),
+		     int (*rf)(double, double),
+		     int (*cmp0)(int),
+		     double x, double y)
+{
+    int tr, rr;
+    rr = (*rf)(x, y);
+    tr = (*tf)(x, y);
+    if (cmp0(rr) != cmp0(tr))
+    {
+	fputs("error:\n\tx     = ", stdout); print_double(x);
+	fputs("\n\ty     = ", stdout); print_double(y);
+	printf("\n\ttrue  = %d\n\tfalse = %d\n", rr, tr);
+    }
+}
+
+
+/*======================================================================*/
+
+
+int main(int ac, char **av)
+{
+#ifdef __alpha__
+    __ieee_set_fp_control(0);
+#endif
+    av++, ac--;
+    switch (*(*av)++)
+    {
+	{
+	    float (*r)(float, float);
+	    float (*t)(float, float);
+
+	    do {
+	      case 'a': r = r_addsf3; t = __addsf3; break;
+	      case 's': r = r_subsf3; t = __subsf3; break;
+	      case 'm': r = r_mulsf3; t = __mulsf3; break;
+	      case 'd': r = r_divsf3; t = __divsf3; break;
+	      case 'r': r = r_sqrtsf3; t = __sqrtsf3; break;
+	      case 'j': r = r_negsf3; t = __negsf3; break;
+	    } while (0);
+
+	    switch (*(*av)++)
+	    {
+	      case 'n':
+		{
+		    int count = (ac > 1 ? atoi(av[1]) : 100);
+		    while (count--)
+			test_float_arith(t, r, rand_float(), rand_float());
+		}
+		break;
+
+	      case 's':
+		{
+		    int i, j;
+		    for (i = 0; i < NSPECIALS; i++)
+			for (j = 0; j < NSPECIALS; j++)
+			    test_float_arith(t, r, gen_special_float(i),
+					      gen_special_float(j));
+		}
+		break;
+
+	      case 0:
+		if (ac < 7) abort();
+		test_float_arith(t, r, build_float(av[1], av[2], av[3]),
+				 build_float(av[4], av[5], av[6]));
+		break;
+	    }
+	}
+	break;
+
+	{
+	    double (*r)(double, double);
+	    double (*t)(double, double);
+
+	    do {
+	      case 'A': r = r_adddf3; t = __adddf3; break;
+	      case 'S': r = r_subdf3; t = __subdf3; break;
+	      case 'M': r = r_muldf3; t = __muldf3; break;
+	      case 'D': r = r_divdf3; t = __divdf3; break;
+	      case 'R': r = r_sqrtdf3; t = __sqrtdf3; break;
+	      case 'J': r = r_negdf3; t = __negdf3; break;
+	    } while (0);
+
+	    switch (*(*av)++)
+	    {
+	      case 'n':
+		{
+		    int count = (ac > 1 ? atoi(av[1]) : 100);
+		    while (count--)
+			test_double_arith(t, r, rand_double(), rand_double());
+		}
+		break;
+
+	      case 's':
+		{
+		    int i, j;
+		    for (i = 0; i < NSPECIALS; i++)
+			for (j = 0; j < NSPECIALS; j++)
+			    test_double_arith(t, r, gen_special_double(i),
+					      gen_special_double(j));
+		}
+		break;
+
+	      case 0:
+		if (ac < 7) abort();
+		test_double_arith(t, r, build_double(av[1], av[2], av[3]),
+				  build_double(av[4], av[5], av[6]));
+		break;
+	    }
+	}
+	break;
+
+      case 'c':
+	switch (*(*av)++)
+	{
+	  case 'n':
+	    {
+		int count = (ac > 1 ? atoi(av[1]) : 100);
+		while (count--)
+		    test_float_double_conv(rand_float());
+	    }
+	    break;
+
+	  case 's':
+	    {
+		int i;
+		for (i = 0; i < NSPECIALS; i++)
+		    test_float_double_conv(gen_special_float(i));
+	    }
+	    break;
+
+	  case 0:
+	    if (ac < 4) abort();
+	    test_float_double_conv(build_float(av[1], av[2], av[3]));
+	    break;
+	}
+	break;
+
+      case 'C':
+	switch (*(*av)++)
+	{
+	  case 'n':
+	    {
+		int count = (ac > 1 ? atoi(av[1]) : 100);
+		while (count--)
+		    test_double_float_conv(rand_double());
+	    }
+	    break;
+
+	  case 's':
+	    {
+		int i;
+		for (i = 0; i < NSPECIALS; i++)
+		    test_double_float_conv(gen_special_double(i));
+	    }
+	    break;
+
+	  case 0:
+	    if (ac < 4) abort();
+	    test_double_float_conv(build_double(av[1], av[2], av[3]));
+	    break;
+	}
+	break;
+
+      case 'i':
+	switch (*(*av)++)
+	{
+	  case 'n':
+	    {
+		int count = (ac > 1 ? atoi(av[1]) : 100);
+		while (count--)
+		    test_int_float_conv(lrand48() << 1);
+	    }
+	    break;
+
+	  case 0:
+	    if (ac < 2) abort();
+	    test_int_float_conv(strtol(av[1], 0, 0));
+	    break;
+	}
+	break;
+
+      case 'I':
+	switch (*(*av)++)
+	{
+	  case 'n':
+	    {
+		int count = (ac > 1 ? atoi(av[1]) : 100);
+		while (count--)
+		    test_int_double_conv(lrand48() << 1);
+	    }
+	    break;
+
+	  case 0:
+	    if (ac < 2) abort();
+	    test_int_double_conv(strtol(av[1], 0, 0));
+	    break;
+	}
+	break;
+
+      case 'f':
+	switch (*(*av)++)
+	{
+	  case 'n':
+	    {
+		int count = (ac > 1 ? atoi(av[1]) : 100);
+		while (count--)
+		    test_float_int_conv(rand_float());
+	    }
+	    break;
+
+	  case 's':
+	    {
+		int i;
+		for (i = 0; i < NSPECIALS; i++)
+		    test_float_int_conv(gen_special_float(i));
+	    }
+	    break;
+
+	  case 0:
+	    if (ac < 4) abort();
+	    test_float_int_conv(build_float(av[1], av[2], av[3]));
+	    break;
+	}
+	break;
+
+      case 'F':
+	switch (*(*av)++)
+	{
+	  case 'n':
+	    {
+		int count = (ac > 1 ? atoi(av[1]) : 100);
+		while (count--)
+		    test_double_int_conv(rand_double());
+	    }
+	    break;
+
+	  case 's':
+	    {
+		int i;
+		for (i = 0; i < NSPECIALS; i++)
+		    test_double_int_conv(gen_special_double(i));
+	    }
+	    break;
+
+	  case 0:
+	    if (ac < 4) abort();
+	    test_double_int_conv(build_double(av[1], av[2], av[3]));
+	    break;
+	}
+	break;
+
+	{
+	    int (*r)(float, float);
+	    int (*t)(float, float);
+	    int (*c)(int);
+
+	    do {
+	      case 'e': r = r_eqsf2; t = __eqsf2; c = eq0; break;
+	      case 'n': r = r_nesf2; t = __nesf2; c = ne0; break;
+	      case 'l':
+		switch (*(*av)++)
+		{
+		  case 'e': r = r_lesf2; t = __lesf2; c = le0; break;
+		  case 't': r = r_ltsf2; t = __ltsf2; c = lt0; break;
+		}
+		break;
+	      case 'g':
+		switch (*(*av)++)
+		{
+		  case 'e': r = r_gesf2; t = __gesf2; c = ge0; break;
+		  case 't': r = r_gtsf2; t = __gtsf2; c = gt0; break;
+		}
+		break;
+	    } while (0);
+
+	    switch (*(*av)++)
+	    {
+	      case 'n':
+		{
+		    int count = (ac > 1 ? atoi(av[1]) : 100);
+		    while (count--)
+			test_float_cmp(t, r, c, rand_float(), rand_float());
+		}
+		break;
+
+	      case 's':
+		{
+		    int i, j;
+		    for (i = 0; i < NSPECIALS; i++)
+			for (j = 0; j < NSPECIALS; j++)
+			    test_float_cmp(t, r, c, gen_special_float(i),
+					   gen_special_float(j));
+		}
+		break;
+
+	      case 0:
+		if (ac < 7) abort();
+		test_float_cmp(t, r, c, build_float(av[1], av[2], av[3]),
+				build_float(av[4], av[5], av[6]));
+		break;
+	    }
+	}
+	break;
+
+	{
+	    int (*r)(double, double);
+	    int (*t)(double, double);
+	    int (*c)(int);
+
+	    do {
+	      case 'E': r = r_eqdf2; t = __eqdf2; c = eq0; break;
+	      case 'N': r = r_nedf2; t = __nedf2; c = ne0; break;
+	      case 'L':
+		switch (*(*av)++)
+		{
+		  case 'E': r = r_ledf2; t = __ledf2; c = le0; break;
+		  case 'T': r = r_ltdf2; t = __ltdf2; c = lt0; break;
+		}
+		break;
+	      case 'G':
+		switch (*(*av)++)
+		{
+		  case 'E': r = r_gedf2; t = __gedf2; c = ge0; break;
+		  case 'T': r = r_gtdf2; t = __gtdf2; c = gt0; break;
+		}
+		break;
+	    } while (0);
+
+	    switch (*(*av)++)
+	    {
+	      case 'n':
+		{
+		    int count = (ac > 1 ? atoi(av[1]) : 100);
+		    while (count--)
+			test_double_cmp(t, r, c, rand_double(), rand_double());
+		}
+		break;
+
+	      case 's':
+		{
+		    int i, j;
+		    for (i = 0; i < NSPECIALS; i++)
+			for (j = 0; j < NSPECIALS; j++)
+			    test_double_cmp(t, r, c, gen_special_double(i),
+					    gen_special_double(j));
+		}
+		break;
+
+	      case 0:
+		if (ac < 7) abort();
+		test_double_cmp(t, r, c, build_double(av[1], av[2], av[3]),
+				build_double(av[4], av[5], av[6]));
+		break;
+	    }
+	}
+	break;
+
+      default:
+	abort();
+    }
+
+    return 0;
+}