/* sqrt.c */
#include <stdio.h>
#include <math.h>

int main(int argc, char * argv[])
{
  double x, y, xp;
  int i;
  
  printf("sqrt.c running \n");
  x = 2.0;
  y = 1.0; /* initial guess */
  printf("find  y = sqrt(x=%g)  y_next = (x/y+y)/2\n", x);
  printf("note solution x = y*y, thus when solved x/y = y \n");
  printf("expect quadratic convergence \n");
  
  for(i=0; i<9; i++)
  {
    y = (x/y+y)/2.0; /* next approximation */
    printf("i=%d, y=%18.15f, err=%g \n", i, y, y-sqrt(x));
  }

  x = 1.0e6;
  y = 1.0;
  printf("find sqrt(x=%g) y_next = (x/y+y)/2\n", x);

  for(i=0; i<20; i++) 
  {
    y = (x/y+y)/2.0; /* next approximation */
    printf("i=%d, y=%18.15f, err=%g \n", i, y, y-sqrt(x));
  }

  /* better  sqrt(x) = sqrt(x/2^2n) * 2^n */
  printf("find sqrt(x=%g) y_next = (x/y+y)/2\n", x);
  printf("xp=x/2^2n, y=sqrt(xp)*2^n \n");
  xp = x/(1024.0*1024.0);
  y=1.0;
  printf("find sqrt(%g) \n", xp);
  for(i=0; i<9; i++) /* expect quadratic convergence */
  {
    y = (xp/y+y)/2.0; /* next approximation */
    printf("i=%d, y=%18.15f, err=%g \n", i, y, y-sqrt(xp));
  }
  y = y*1024.0; /* 2^10 */
  printf("y=y*2^10=%18.15f, err=%g \n", y, y-sqrt(x)); 
  
  return 0;
}