/***********************************************************************
Test the double precision IEEE 754 floating-point support functions.
This test takes a long time:

   SGI Indigo (-O):		 707.5u 105.3s 14:19 94%
   Sun SPARCstation SLC (-g):	2180.2u 718.2s 51:07 94% 0+152k 45+50io 0pf+0w
   Sun 386i (-O -I.):		13272.3u 12.9s 3:51:05 95% 0+112k 0+29io 0pf+0w
[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_DP - T_DP + 1)
#define EXPONENT_MAX	(EXPONENT_INFNAN_DP - 1)

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

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

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

    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(-1076));
    test_all("neg min denormal", -BASE_TO_THE(-1076));

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

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

    for (k = EXPONENT_DENORM_DP + 1; k < (EXPONENT_INFNAN_DP -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(-T_DP)) * BASE_TO_THE(1023));
    test_all("neg max normal",-(2.0 - BASE_TO_THE(-T_DP)) * BASE_TO_THE(1023));

    test_all("+Infinity", (w.i[DP_HIGH] = Inf_DP, w.i[DP_LOW] = 0, w.r));
    test_all("-Infinity", (w.i[DP_HIGH] = NegInf_DP, w.i[DP_LOW] = 0, w.r));
    test_all("NaN", (w.i[DP_HIGH] = NaN_DP, w.i[DP_LOW] = 0, w.r));

    exit(EXIT_SUCCESS);
    return (EXIT_SUCCESS);
}

#if STDC
DOUBLEPRECISION
deps2(DOUBLEPRECISION *x)
#else /* NOT STDC */
DOUBLEPRECISION
deps2(x)
DOUBLEPRECISION	*x;
#endif /* STDC */
{
    DOUBLEPRECISION		half = 0.5;
    INTEGER			n = dintxp(x) - T_DP;
    INTEGER			zero = 0;
    DOUBLEPRECISION_PARTS	w;

    if ( (*x == 0.0) || disden(x) )
    {
	w.i[DP_HIGH] = MIN_DENORMAL_DP;
	w.i[DP_LOW] = MIN_DENORMAL_Low_DP;
	return (w.r);
    }
    else if ( disnan(x) || disinf(x) )
	return *x;
    else if ( (*x < 0.0) && (dsetxp(x,&zero) == -half) )
    {					/* special boundary case */
	w.r = dsetxp(&half,(n--, &n));
	if (w.r == 0.0)		/* then x = 2.22507e-308 0x00000000 00100000 */
	{
	    w.i[DP_HIGH] = MIN_DENORMAL_DP; /* and setxp() -> 0, so fix up */
	    w.i[DP_LOW] = MIN_DENORMAL_Low_DP;
	}
	return (w.r);
    }
    else
	return (dsetxp(&half,&n));
}

#if STDC
DOUBLEPRECISION
scale(DOUBLEPRECISION value, INTEGER n)	/* return value * 2**n */
#else /* NOT STDC */
DOUBLEPRECISION
scale(value,n)				/* return value * 2**n */
DOUBLEPRECISION value;
INTEGER n;
#endif /* STDC */
{
    INTEGER k;
    DOUBLEPRECISION 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(DOUBLEPRECISION x)
#else /* NOT STDC */
char*
show(x)
DOUBLEPRECISION x;
#endif /* STDC */
{
    static char buffer[128];
    DOUBLEPRECISION_PARTS w;
    w.r = x;
#if defined(__alpha)
    (void)sprintf(buffer, "%15g [0x%08x %08x]",
		  w.r,
		  (unsigned int)w.i[DP_HIGH],
		  (unsigned int)w.i[DP_LOW]);
#else
    (void)sprintf(buffer, "%15g [0x%08lx %08lx]",
		  w.r,
		  (unsigned long)w.i[DP_HIGH],
		  (unsigned long)w.i[DP_LOW]);
#endif
    return (buffer);
}

#if STDC
void
show_eps(DOUBLEPRECISION x, DOUBLEPRECISION e_f, DOUBLEPRECISION e_c)
#else /* NOT STDC */
void
show_eps(x, e_f, e_c)
DOUBLEPRECISION x;
DOUBLEPRECISION e_f;
DOUBLEPRECISION 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(DOUBLEPRECISION value)
#else /* NOT STDC */
void
test_adx(value)
DOUBLEPRECISION value;
#endif /* STDC */
{
    INTEGER n;
    DOUBLEPRECISION y;
    DOUBLEPRECISION z;

    for (n = EXPONENT_MIN; n <= EXPONENT_MAX; ++n)
    {
	y = scale(value,n);
	z = dadx(&value,&n);
	if (!(disnan(&y) && disnan(&z)) && (y != z))
	{
	    (void)printf("\tdadx(%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, DOUBLEPRECISION value)
#else /* NOT STDC */
void
test_all(name,value)
const char *name;
DOUBLEPRECISION 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
void
test_eps(DOUBLEPRECISION v)
#else /* NOT STDC */
void
test_eps(v)
DOUBLEPRECISION v;
#endif /* STDC */
{
    INTEGER n;
    DOUBLEPRECISION value = v;
    DOUBLEPRECISION x;
    DOUBLEPRECISION e_f;
    DOUBLEPRECISION e_c;

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

	e_f = deps(&x);
	e_c = deps2(&x);

	if (disnan(&x))
	{
	    if (!(disnan(&e_f) && disnan(&e_c)))
		show_eps(x, e_f, e_c);
	}
	else if (disinf(&x))
	{
	    if (!(disinf(&e_f) && disinf(&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(DOUBLEPRECISION value)
#else /* NOT STDC */
void
test_intxp(value)
DOUBLEPRECISION value;
#endif /* STDC */
{
    INTEGER e;
    INTEGER n;
    INTEGER m;
    DOUBLEPRECISION x;

    m = dintxp(&value);
    for (n = EXPONENT_MIN; n <= EXPONENT_MAX; ++n)
    {
	x = scale(value,n);
	e = dintxp(&x);
	if (x == 0.0)
	{
	    if (e != 0)
	        (void)printf("\tdintxp(%s) = %d != %d\n", show(x), (int)e, (int)n);
	}
	else if (disnan(&x) || disinf(&x))
	{
	    if (e != EXPONENT_INFNAN_DP)
	        (void)printf("\tdintxp(%s) = %d != %d\n", show(x), (int)e, (int)n);
	}
        else if (e != (n + m))
	    (void)printf("\tdintxp(%s) = %d != %d\n", show(x), (int)e, (int)n);
    }
}

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

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