#include #include "ieeeftn.h" #include /* 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; }