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

/**************************************************
*     (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,
*     deps(&1.0) is known as the largest relative
*     spacing, beta**(1-t), and deps(&-1.0) is the
*     smallest relative spacing, beta**(-t).
*
*     (14-Feb-1992)
**************************************************/

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

    if (disinf(&x))
	epsilon = x;
    else if (disnan(&x))
	epsilon = x;
    else
    {
	epsilon = (x == 0.0) ? 1.0 : ABS(x);
	while ((temp = epsilon*betain, temp += x, temp != x))
	    epsilon *= betain;

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