/*
** eps2.c
** Mth 351 Spring 2001 - Bent Petersen
**
** Microsoft C/C++ cl 32-bit 12.00.8168      cl /Od eps2.c
** GNU C/C++ gcc egcs-2.91.66 (egcs 1.1.2)   gcc -o eps2 -i eps2.c
**
** This code looks pretty silly but computing the unit
** round (epsilon) is a bit tricky. The obvious code may
** not work, or may work but yield incorrect results. It 
** is necessary to force stores of intermediate results to 
** prevent the compiler from doing the calculations at the 
** full precision of the FPU or, at any rate, at too high
** precision.
*/

#include <stdio.h>
#include <float.h>

int main( void )  {

/* gcc wants main() to return an int - void main(void) produces */
/* a warning message. MS cl doesn't care. */

  float x1=1.0, z1=2.0;
  double x2=1.0, z2=2.0;
  long double x3=(long double)1.0, z3=(long double)2.0;
  float x4=(float)1.0;

  while ( z1 > 1.0 )  {
    x1 /= 2.0;
    z1 = 1.0 + x1;
}
  printf("\nUnit round - float: %e\t(float.h claims %e)", 
        2*x1, FLT_EPSILON );

  while ( z2 > 1.0 )  {
    x2 /= 2.0;
    z2 = 1.0 + x2;
}
  printf("\nUnit round - double: %e\t(float.h claims %e)",
        2*x2, DBL_EPSILON );

  while ( z3 > 1.0 )  {
    x3 /= 2.0;
    z3 = 1.0 + x3;
}
  printf("\nUnit round - long double: %e\t(float.h claims %e)",
        2.0*(double)x3, (double)LDBL_EPSILON );

/* In float.h for gcc LDBL_EPSILON is declared as a long double */
/* and therefore the cast to a double is required for printf(). */
/* Likewise x3 must be cast to a double or printf() will print  */
/* nonsense. It is a bit curious since LDBL_EPSILON easily fits */
/* in a double. */

  while ( 1.0 + x4 > 1.0 )
    x4 /= 2.0 ;

  printf("\nUnit round - %e\t\t(bogus, float without stores)\n",
        2.0*x4 );
  
  printf("\n");
  return 0;
}