// test frexp() and Newton square-root
//
#include <stdio.h>
#include <stdlib.h> // atof()
#include <math.h>

int main( int argc, char *argv[])
{
  double a = 2; if( argc > 1) a = atof( argv[1]);

  int e; double f = frexp( a, &e);

  printf( "a = %g, f = %g, e = %i\n", a, f, e);

  double x = a/2; // like a3
  double y1 = sqrt(a);
  double y2 = (x + a/x)/2; // Newton

  printf( " x = %.18g\n", x); // %.18f format better only if numbers are near 1.0
  printf( "y1 = %.18g\n", y1);
  printf( "y2 = %.18g\n", y2);

 // using f,e from frexp
 //
  if( e % 2 == 1) { f *= 2; --e; } // e was odd, now it is even

  x = f/2;
  double g = (x + f/x)/2; // Newton, g ~= sqrt(f)
  y2 = ldexp( g, e/2);
  printf( "y2 = %.18g\n", y2);
}

/* sample runs:

$ ./a.out
a = 2, f = 0.5, e = 2
 x = 1
y1 = 1.41421356237309515
y2 = 1.5
y2 = 2.25

$ ./a.out 4
a = 4, f = 0.5, e = 3
 x = 2
y1 = 2
y2 = 2
y2 = 2.5

$ ./a.out 6
a = 6, f = 0.75, e = 3
 x = 3
y1 = 2.44948974278317788
y2 = 2.5
y2 = 2.75

$ ./a.out 10000
a = 10000, f = 0.610352, e = 14
 x = 5000
y1 = 100
y2 = 2501
y2 = 147.53125

$ ./a.out 1000000
a = 1e+06, f = 0.953674, e = 20
 x = 500000
y1 = 1000
y2 = 250001
y2 = 1268.140625

*/