/***********************************************************************
Test the single precision IEEE 754 floating-point support functions.
[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))

/* NB: In order to work with pre-Standard C compilers, values of type
REAL are passed under the name `v', and then copied internally to a
REAL under the name `value'.  That way, promotion to double for
argument passing does not interfere with subsequent access to the bits
of the REAL. */

int	main ARGS((void));
REAL	scale ARGS((REAL v__, INTEGER n__));
char*	show ARGS((REAL x__));
void	show_eps ARGS((REAL x, REAL e_f, REAL e_c));
void	test_adx ARGS((REAL v__));
void	test_all ARGS((const char *name__, REAL v__));
void	test_eps ARGS((REAL v__));
void	test_intxp ARGS((REAL v__));
void	test_setxp ARGS((REAL v__));

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
REAL
eps2(REAL *x)
#else /* NOT STDC */
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 */
	    /* and setxp() -> 0, so fix up */
	    w.i = MIN_DENORMAL_SP;
	return (w.r);
    }
    else
	return (setxp(&half,&n));
}

#if STDC
REAL
scale(REAL v, INTEGER n)		/* return v * 2**n */
#else /* NOT STDC */
REAL
scale(v,n)
REAL v;
INTEGER n;				/* return v * 2**n */
#endif /* STDC */
{
    INTEGER k;
    REAL value = v;			/* force REAL precision */
    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
char*
show(REAL x)
#else /* NOT STDC */
char*
show(x)
REAL x;
#endif /* STDC */
{
    static char buffer[128];
    REAL_PARTS w;

    w.r = x;
#if defined(__alpha)
    (void)sprintf(buffer, "%15g\t0x%08x", w.r, (unsigned int)w.i);
#else
    (void)sprintf(buffer, "%15g\t0x%08lx", w.r, (unsigned long)w.i);
#endif
    return (buffer);
}

#if STDC
void
show_eps(REAL x, REAL e_f, REAL e_c)
#else /* NOT STDC */
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
void
test_adx(REAL v)
#else /* NOT STDC */
void
test_adx(v)
REAL v;
#endif /* STDC */
{
    INTEGER n;
    REAL value = v;			/* force REAL precision */
    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
void
test_all(const char *name, REAL v)
#else /* NOT STDC */
void
test_all(name,v)
const char *name;
REAL v;
#endif /* STDC */
{
    REAL value = v;

    (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
void
test_eps(REAL v)
#else /* NOT STDC */
void
test_eps(v)
REAL v;
#endif /* STDC */
{
    INTEGER n;
    REAL value = v;
    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
void
test_intxp(REAL v)
#else /* NOT STDC */
void
test_intxp(v)
REAL v;
#endif /* STDC */
{
    INTEGER e;
    INTEGER n;
    INTEGER m;
    REAL value = v;			/* force REAL precision */
    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
void
test_setxp(REAL v)
#else /* NOT STDC */
void
test_setxp(v)
REAL v;
#endif /* STDC */
{
    INTEGER k;
    INTEGER n;
    REAL value = v;			/* force REAL precision */
    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));
	}
    }
}