/***********************************************************************
[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 MAXTIMES 10

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

int		main(VOID_ARG);
DOUBLEPRECISION	scale ARGS((DOUBLEPRECISION value, INTEGER n));
double		second ARGS((void));
void		time_deps ARGS((DOUBLEPRECISION value__));

int
main(VOID_ARG)
{
    DOUBLEPRECISION_PARTS w;
    double t;

    t = second();
    time_deps(0.0);
    time_deps(0.5);
    time_deps(-0.5);
    time_deps(128.0);
    time_deps(-128.0);
    time_deps(1.0/128.0);
    time_deps(-1.0/128.0);
    time_deps((w.i[DP_HIGH] = Inf_DP,
	       w.i[DP_LOW] = Inf_Low_DP, w.r));
    time_deps((w.i[DP_HIGH] = NegInf_DP,
	       w.i[DP_LOW] = NegInf_Low_DP, w.r));
    time_deps((w.i[DP_HIGH] = NaN_DP,
	       w.i[DP_LOW] = NaN_Low_DP, w.r));
    time_deps((w.i[DP_HIGH] = MIN_NORMAL_DP,
	       w.i[DP_LOW] = MIN_NORMAL_Low_DP, w.r));
    time_deps((w.i[DP_HIGH] = MIN_NORMAL_DP,
	       w.i[DP_LOW] = MIN_NORMAL_Low_DP, -w.r));
    time_deps((w.i[DP_HIGH] = MIN_DENORMAL_DP,
	       w.i[DP_LOW] = MIN_DENORMAL_Low_DP, w.r));
    time_deps((w.i[DP_HIGH] = MIN_DENORMAL_DP,
	       w.i[DP_LOW] = MIN_DENORMAL_Low_DP, -w.r));
    time_deps((w.i[DP_HIGH] = MAX_NORMAL_DP,
	       w.i[DP_LOW] = MAX_NORMAL_Low_DP, w.r));
    time_deps((w.i[DP_HIGH] = MAX_NORMAL_DP,
	       w.i[DP_LOW] = MAX_NORMAL_Low_DP, -w.r));
    t = second() - t;			/* 16 * MAXTIMES deps() calls */
    (void)printf("Average deps() time = %10.2lf microsec/call\n",
		 t * 1000000.0 / (16 * (double)MAXTIMES *
				  (double)(1 + EXPONENT_MAX - EXPONENT_MIN)));

    exit (EXIT_SUCCESS);
    return (EXIT_SUCCESS);
}

#if STDC
DOUBLEPRECISION
scale(DOUBLEPRECISION value,INTEGER 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
void
store(double d)
#else
void
store(d)
double d;
#endif
{
}

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

    for (n = EXPONENT_MIN; n <= EXPONENT_MAX; ++n)
    {
	x = scale(value,n);
	for (k = 0; k < MAXTIMES; ++k)
	    store(deps(&x));
    }
}