/***********************************************************************
Test the single precision IEEE 754 floating-point support functions

Timing:
  Sun 386i (SunOS 4.0.3):	236.8u 0.1s 4:06 96% 0+32k 0+10io 0pf+0w
  IBM RS/6000-550 (AIX 3.1):	9.22u 0.7s 0:09 98% 8+30k 0+0io 0pf+0w

Notes:
	If gcc 1.36 is used for isnan(), it must be compiled with -g
	on Sun 386i; otherwise, it returns true when the arg is a NaN!

[20-Feb-1994]
***********************************************************************/

#include <stdio.h>
#include <math.h>
#include "ieeeftn.h"

/* 0.5 * 2**e is representable for EXPONENT_MIN <= e <= EXPONENT_MAX */
#define EXPONENT_MIN	(EXPONENT_DENORM_SP - T_SP + 1)
#define EXPONENT_MAX	(EXPONENT_INFNAN_SP - 1)

#define POW(x,y)	(REAL)pow((double)(x),(double)(y))
#define BASE_TO_THE(k)	POW(2.0,(k))

int		main ARGS((void));
REAL		eps ARGS((REAL *x));		/* Fortran version */
static REAL	eps2 ARGS((REAL *x));
static REAL	scale ARGS((REAL value, INTEGER n));
static char	*show ARGS((REAL x));
static void	show_eps ARGS((REAL x, REAL e_f, REAL e_c));
static void	test_adx ARGS((REAL value));
static void	test_all ARGS((const char *name, REAL value));
static void	test_eps ARGS((REAL value));
static void	test_intxp ARGS((REAL value));
static void	test_setxp ARGS((REAL value));

int
main(VOID_ARG)
{
    REAL_PARTS w;
    INTEGER k;
    char message[256];

    test_all("0.5", 0.5);
    test_all("128.0", 128.0);
    test_all("1.0/128.0", 1.0/128.0);

    test_all("zero", 0.0);

    test_all("min denormal", BASE_TO_THE(-149));
    test_all("neg min denormal", -BASE_TO_THE(-149));

    for (k = 1; k <= T_SP; ++k)
    {			/* test all bit positions in denormalized numbers */
	(void)sprintf(message,"min denormal * 2**%d", (int)k);
	test_all(message, BASE_TO_THE(-149 + k));
	(void)sprintf(message,"neg min denormal * 2**%d", (int)k);
	test_all(message, -BASE_TO_THE(-149 + k));
    }

    test_all("min normal", BASE_TO_THE(-126));
    test_all("neg min normal", -BASE_TO_THE(-126));

    for (k = EXPONENT_DENORM_SP + 1; k < (EXPONENT_INFNAN_SP -1); ++k)
    {			/* test all exponent values in normal numbers */
	(void)sprintf(message,"+2**%d", (int)k);
	test_all(message, BASE_TO_THE(k));
	(void)sprintf(message,"-2**%d", (int)k);
	test_all(message, -BASE_TO_THE(k));
    }

    test_all("max normal", (2.0 - BASE_TO_THE(-23)) * BASE_TO_THE(127));
    test_all("neg max normal", -(2.0 - BASE_TO_THE(-23)) * BASE_TO_THE(127));

    test_all("+Infinity", (w.i = Inf_SP, w.r));
    test_all("-Infinity", (w.i = NegInf_SP, w.r));
    test_all("NaN", (w.i = NaN_SP, w.r));

    exit(EXIT_SUCCESS);
    return (EXIT_SUCCESS);
}

#if STDC
static REAL
eps2(REAL *x)
#else /* NOT STDC */
static REAL
eps2(x)
REAL	*x;
#endif /* STDC */
{
    REAL	half = 0.5;
    INTEGER	n = intxp(x) - T_SP;
    INTEGER	zero = 0;
    REAL_PARTS	w;

    if ( (*x == 0.0) || isden(x) )
    {
	w.i = MIN_DENORMAL_SP;
	return (w.r);
    }
    else if ( isnan(x) || isinf(x) )
	return *x;
    else if ( (*x < 0.0) && (setxp(x,&zero) == -half) )
    {					/* special boundary case */
	w.r = setxp(&half,(n--, &n));
	if (w.r == 0.0)			/* then x = -1.17549e-38 80800000 */
	    w.i = MIN_DENORMAL_SP;	/* and setxp() -> 0, so fix up */
	return (w.r);
    }
    else
	return (setxp(&half,&n));
}

#if STDC
static REAL
scale(REAL value, INTEGER n)		/* return value * 2**n */
#else /* NOT STDC */
static REAL
scale(value,n)
REAL value;
INTEGER n;
#endif /* STDC */
{
    INTEGER k;
    REAL y;

    if (n > (EXPONENT_MAX - 2))		/* then 2**n overflows */
    {
	k = n - (EXPONENT_MAX - 2);	/* k > 0 */
	y = (value * BASE_TO_THE(n - k)) * BASE_TO_THE(k);
    }
    else if (n < EXPONENT_MIN)		/* then 2**n underflows to zero */
    {
	k = n - EXPONENT_MIN;		/* k < 0 */
	y = (value * BASE_TO_THE(k)) * BASE_TO_THE(EXPONENT_MIN);
    }
    else				/* 2**n is representable */
	y = value * BASE_TO_THE(n);
    return (y);
}

#if STDC
static char*
show(REAL x)
#else /* NOT STDC */
static char*
show(x)
REAL x;
#endif /* STDC */
{
    static char buffer[128];
    REAL_PARTS w;
    w.r = x;
    (void)sprintf(buffer, "%lg [0x%08lx]", (double)w.r, (unsigned long)w.i);
    return (buffer);
}

#if STDC
static void
show_eps(REAL x, REAL e_f, REAL e_c)
#else /* NOT STDC */
static void
show_eps(x,e_f,e_c)
REAL x;
REAL e_f;
REAL e_c;
#endif /* STDC */
{
    (void)printf("\tFortran eps(%s) = ",show(x));
    (void)printf("%s != ",show(e_f));
    (void)printf("C eps() = %s\n",show(e_c));
}

#if STDC
static void
test_adx(REAL value)
#else /* NOT STDC */
static void
test_adx(value)
REAL value;
#endif /* STDC */
{
    INTEGER n;
    REAL y;
    REAL z;

    for (n = EXPONENT_MIN; n <= EXPONENT_MAX; ++n)
    {
	y = scale(value,n);
	z = adx(&value,&n);
	if (!(isnan(&y) && isnan(&z)) && (y != z))
	{
	    (void)printf("\tadx(%s,%d) = ", show(value), (int)n);
	    (void)printf("%s != ", show(z));
	    (void)printf("%s\n", show(y));
	}
    }
}

#if STDC
static void
test_all(const char *name, REAL value)
#else /* NOT STDC */
static void
test_all(name,value)
const char *name;
REAL value;
#endif /* STDC */
{
    (void)printf("Test value: %s\t%s\n", show(value), name);
    test_adx(value);
    test_intxp(value);
    test_setxp(value);
    test_eps(value);
}

#if STDC
static void
test_eps(REAL value)
#else /* NOT STDC */
static void
test_eps(value)
REAL value;
#endif /* STDC */
{
    INTEGER n;
    REAL x;
    REAL e_f;
    REAL e_c;

    for (n = EXPONENT_MIN; n <= EXPONENT_MAX; ++n)
    {
	x = scale(value,n);

	e_f = eps(&x);
	e_c = eps2(&x);

	if (isnan(&x))
	{
	    if (!(isnan(&e_f) && isnan(&e_c)))
		show_eps(x, e_f, e_c);
	}
	else if (isinf(&x))
	{
	    if (!(isinf(&e_f) && isinf(&e_c)))
		show_eps(x, e_f, e_c);
	}
        else if (e_f != e_c)
	    show_eps(x, e_f, e_c);
    }
}

#if STDC
static void
test_intxp(REAL value)
#else /* NOT STDC */
static void
test_intxp(value)
REAL value;
#endif /* STDC */
{
    INTEGER e;
    INTEGER n;
    INTEGER m;
    REAL x;

    m = intxp(&value);
    for (n = EXPONENT_MIN; n <= EXPONENT_MAX; ++n)
    {
	x = scale(value,n);
	e = intxp(&x);
	if (x == 0.0)
	{
	    if (e != 0)
	        (void)printf("\tintxp(%s) = %d != %d\n", show(x), (int)e, (int)n);
	}
	else if (isnan(&x) || isinf(&x))
	{
	    if (e != EXPONENT_INFNAN_SP)
	        (void)printf("\tintxp(%s) = %d != %d\n", show(x), (int)e, (int)n);
	}
        else if (e != (n + m))
	    (void)printf("\tintxp(%s) = %d != %d\n", show(x), (int)e, (int)n);
    }
}

#if STDC
static void
test_setxp(REAL value)
#else /* NOT STDC */
static void
test_setxp(value)
REAL value;
#endif /* STDC */
{
    INTEGER k;
    INTEGER n;
    REAL x;
    REAL y;
    REAL z;

    k = intxp(&value);
    x = scale(value,-k);		/* x = fraction(value) * 2**0 */
    for (n = EXPONENT_MIN; n <= EXPONENT_MAX; ++n)
    {
	y = scale(x,n);
	z = setxp(&value,&n);
	if (!(isnan(&value) && isnan(&z)) && (y != z))
	{
	    (void)printf("\tsetxp(%s,%d) = ", show(value), (int)n);
	    (void)printf("%s != ", show(z));
	    (void)printf("%s\n", show(y));
	}
    }
}