/***********************************************************************
Time the direct and iterative eps() 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 MAXTIMES 10

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

int	main ARGS((int argc__, char* argv__[]));
REAL	scale ARGS((REAL value, INTEGER n));
double	second ARGS((void));
void	time_eps ARGS((REAL value__));

#if STDC
int
main(int argc,char* argv[])
#else /* NOT STDC */
int
main(argc,argv)
int argc;
char* argv[];
#endif /* STDC */
{
    REAL_PARTS w;
    double t;

    t = second();
    time_eps(0.0);
    time_eps(0.5);
    time_eps(-0.5);
    time_eps(128.0);
    time_eps(-128.0);
    time_eps(1.0/128.0);
    time_eps(-1.0/128.0);
    time_eps((w.i = Inf_SP, w.r));
    time_eps((w.i = NegInf_SP, w.r));
    time_eps((w.i = NaN_SP, w.r));
    time_eps((w.i = MIN_NORMAL_SP, w.r));
    time_eps((w.i = MIN_NORMAL_SP, -w.r));
    time_eps((w.i = MIN_DENORMAL_SP, w.r));
    time_eps((w.i = MIN_DENORMAL_SP, -w.r));
    time_eps((w.i = MAX_NORMAL_SP, w.r));
    time_eps((w.i = MAX_NORMAL_SP, -w.r));
    t = second() - t;			/* 16 * MAXTIMES eps() calls */
    (void)printf("Average eps() 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
REAL
scale(REAL value, INTEGER n)		/* return value * 2**n */
#else /* NOT STDC */
REAL
scale(value,n)		/* return value * 2**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
void
store(float f)
#else
void
store(f)
float f;
#endif
{
}

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

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