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

#include <stdio.h>			/* DEBUG */


/**************************************************
*     (Machine Epsilon)
*
*     Return the machine epsilon, the smallest
*     positive number such that (x + epsilon) .ne.
*     x.  This code works for x = 0.0, NaN, and
*     Inf as well.
*
*     Note that by supplying negative x values,
*     you can determine the smallest positive
*     number such that (x - epsilon) .ne. x.
*
*     For base beta and t-digit significands,
*     eps(&1.0) is known as the largest relative
*     spacing, beta**(1-t), and eps(&-1.0) is the
*     smallest relative spacing, beta**(-t).
*
*     (14-Feb-1992)
**************************************************/

void store ARGS((REAL *x));

#if STDC
REAL
eps(REAL *px)
#else /* NOT STDC */
REAL
eps(px)
REAL *px;
#endif /* STDC */
{
    REAL beta = 2.0;
    REAL betain = 1.0/beta;
    REAL epsilon;
    REAL temp;
    REAL x = *px;

    if (isinf(&x))
	epsilon = x;
    else if (isnan(&x))
	epsilon = x;
    else
    {
	epsilon = (x == 0.0) ? 1.0 : ABS(x);
#if 1
	while ((temp = epsilon*betain, temp += x, temp != x))
	    epsilon *= betain;
#else
	for (;;)
	{
	    {
		union {float r; unsigned int i; } wf;
		wf.r = epsilon;
	    }
	    if (epsilon == 0.0)
		break;
	    temp = epsilon*betain;
	    temp += x;
	    store(&temp);
	    if (temp == x)
		break;
#if 0
	    {
		union {float r; unsigned int i; } wf;
		wf.r = x;
		printf("x = %15.7g %08x\t", wf.r, wf.i);
		wf.r = temp;
		printf("temp = %15.7g %08x\t", wf.r, wf.i);
		wf.r = temp + x;
		printf("temp+x = %15.7g %08x\n",wf.r, wf.i);
	    }
#endif
	    epsilon *= betain;
	    store(&epsilon);
	}
#endif

	/* purify epsilon to eliminate all but high-order bit */
	temp = x + epsilon;
	if (isinf(&temp))
	{
	    epsilon = x - epsilon;
	    epsilon = x - epsilon;
	}
	else
	{
	    epsilon = x + epsilon;
	    epsilon = epsilon - x;
	}
    }
    return (epsilon);
}

#if STDC
void
store(REAL *x)
#else /* NOT STDC */
void
store(x)
REAL *x;
#endif /* STDC */
{
    *x = *x;
}