#include<stdio.h>
#include<stdlib.h>
 
int main(void)
{
  double x0,x1,a;
  int i;
 
  x1 = 1;
  x0  = 0;
  printf("Find square root of: ");
  fflush(stdout);
  scanf("%lf",&a);
  for(i = 0; i < 100 && x0 != x1; ++i)
  {
    x0 = x1;
    x1 = .5 * (x1 + a / x1);
#ifdef INCLUDE_PRINTF
    printf("Iteration %d - converging to %f - prev value %f\n",i,x1,x0);
#endif
  }
 
  printf("Finished in %d iterations\n",i);
  printf("Results: %.15f, %.15f\n",x0,x1);
  {
    union { double a; unsigned long b[2]; } v;
    v.a = x0;
    printf( "X0 in binary = %lu %lu\n", v.b[0], v.b[1] );
    v.a = x1;
    printf( "X1 in binary = %lu %lu\n", v.b[0], v.b[1] );
  }
  return 0;
}

$ gcc -W -Wall -ansi -pedantic -O2 -DINCLUDE_PRINTF foo.c
$ ./a.exe
Find square root of: 5
Iteration 0 - converging to 3.000000 - prev value 1.000000
Iteration 1 - converging to 2.333333 - prev value 3.000000
Iteration 2 - converging to 2.238095 - prev value 2.333333
Iteration 3 - converging to 2.236069 - prev value 2.238095
Iteration 4 - converging to 2.236068 - prev value 2.236069
Iteration 5 - converging to 2.236068 - prev value 2.236068
Iteration 6 - converging to 2.236068 - prev value 2.236068
Finished in 7 iterations
Results: 2.236067977499790, 2.236067977499790
X0 in binary = 2610427048 1073865591
X1 in binary = 2610427048 1073865591
 
$ gcc -W -Wall -ansi -pedantic -O2 foo.c
$ ./a.exe
Find square root of: 5
Finished in 100 iterations
Results: 2.236067977499790, 2.236067977499790
X0 in binary = 2610427048 1073865591
X1 in binary = 2610427048 1073865591
 
$ gcc -W -Wall -ansi -pedantic -O2 -ffloat-store foo.c
$ ./a.exe
Find square root of: 5
Finished in 7 iterations
Results: 2.236067977499790, 2.236067977499790
X0 in binary = 2610427048 1073865591
X1 in binary = 2610427048 1073865591
$

`-ffloat-store'
     Do not store floating point variables in registers, and inhibit
     other options that might change whether a floating point value is
     taken from a register or memory.
 
     This option prevents undesirable excess precision on machines such
     as the 68000 where the floating registers (of the 68881) keep more
     precision than a `double' is supposed to have.  Similarly for the
     x86 architecture.  For most programs, the excess precision does
     only good, but a few programs rely on the precise definition of
     IEEE floating point.  Use `-ffloat-store' for such programs, after
     modifying them to store all pertinent intermediate computations
     into variables.
 

$ gcc -W -Wall -ansi -pedantic -O2 -ffloat-store -Wfloat-equal foo.c
foo.c: In function `main':
foo.c:14: warning: comparing floating point with == or != is unsafe
 

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

int main(void)
{
   double x0 = 0, x1 = 1, a;
   int i;
   printf("Find square root of: ");
   fflush(stdout);
   scanf("%lf",&a);
   for ( i = 0; i < 100 && x0 != x1; ++i )
   {
      x0 = x1;
      x1 = .5 * (x1 + a / x1);
#ifdef INCLUDE_PRINTF
      printf("Iteration %d - converging to %f - prev value %f\n",i,x1,x0);
#endif
      {
         double epsilon = DBL_EPSILON * fabs(x1);
         if ( fabs(x1 - x0) <= epsilon )
         {
            printf("reached epsilon of %g\n", epsilon);
            break;
         }
      }
   }
   printf("Finished in %d iterations\n",i);
   printf("Results: %.15f, %.15f\n",x0,x1);
   return 0;
}