// epsilon.c  float 32 bit binary, 6 to 7 decimal digits
//           double 64 bit binary, 15 to 16 decimal digits
//           epsilon is smallest number added to 1.0 that is greater than 1.0
//           floating point registers are more that 64 bits, 80 or more
//           thus floating point arithmetic in registers may have 19 digits

#include <stdio.h>

int main()
{
  double aa, bb, cc, dd;
  float  a,  b,  c,  d;
  printf("epsilon.c running \n");
  aa = 1.0;      // note 1.0 is 64 bit constant
  bb = 1.0e-15;
  cc = aa + bb;  // add done in long registers, stored as 64 bit
  dd = cc - 1.0;

  printf("first double 1.0e-15 greater than double epsilon \n");
  printf("double cc = aa + bb;  dd = cc - 1.0 \n");
  printf("aa=%g  bb=%g  cc=%g \n",aa,bb,cc);
  printf("aa=%e  bb=%e  cc=%e \n",aa,bb,cc);
  printf("aa=%18.16e  bb=%18.16e  cc=%18.16e \n",aa,bb,cc);
  printf("dd=%g  dd=%e  dd=%17.15e \n\n", dd, dd, dd);

  a = 1.0f;    // note 1.0f is 32 bit constant
  b = 1.0e-7f;
  c = a + b;   // add done in long registers, stored as 32 bit
  d = c - 1.0f;

  printf("next float 1.0e-7 greater than float epsilon \n");
  printf("float c = a + b;  d = c - 1.0 \n");
  printf("a=%g  b=%g  c=%g \n",a,b,c);
  printf("a=%e  b=%e  c=%e \n",a,b,c);
  printf("a=%17.15e  b=%17.15e  c=%17.15e \n",a,b,c);
  printf("d=%g  d=%e  d=%17.15e \n\n", d, d, d);
  
  bb = 1.0e-16;
  cc = aa + bb;  // add done in long registers, stored as 64 bit
  dd = cc - 1.0;

  printf("then double 1.0e-16 smaller than double epsilon \n");
  printf("double cc = aa + bb;  dd = cc - 1.0 \n");
  printf("aa=%g  bb=%g  cc=%g \n",aa,bb,cc);
  printf("aa=%e  bb=%e  cc=%e \n",aa,bb,cc);
  printf("aa=%18.16e  bb=%18.16e  cc=%18.16e \n",aa,bb,cc);
  printf("dd=%g  dd=%e  dd=%17.15e \n", dd, dd, dd);
  printf("note dd = 0.0 \n\n");

  b = 1.0e-8f;
  c = a + b;   // add done in long registers, stored as 32 bit
  d = c - 1.0f;

  printf("then float 1.0e-8 smaller than float epsilon \n");
  printf("float c = a + b;  d = c - 1.0 \n");
  printf("a=%g  b=%g  c=%g \n",a,b,c);
  printf("a=%e  b=%e  c=%e \n",a,b,c);
  printf("a=%17.15e  b=%17.15e  c=%17.15e \n",a,b,c);
  printf("d=%g  d=%e  d=%17.15e \n", d, d, d);
  printf("note d = 0.0 \n\n");
  
  printf("epsilon.c finished \n");
  return 0;
}