// -----------------------------
//  iterate1.cpp
// -----------------------------
 #include "w.h"

 // Newton for x^2 = 2   =>  f(x) = x/2 + 1/x
 double f(double x) { return (x*0.5 + 1.0/x); }

 // Newton for exp(x) = 2  =>  g(x) = x -1 + 2/exp(x)
 double g(double x) { return ( x - 1 + 2/exp(x)); }

 // automate function iteration
 double fiterate(fptype1 fi, double x0, int n = 3)
 {
	double x = fi(x0);
	for (int i = 0; i < n; i++) x = fi(x);
	return x;
 }

void main()
{
  nl(0);
  banner("iterate1.cpp");

  // send message to output stream
  cout.precision(20);

  // double is 64 bit (about 15 sig fig)

  // check: our pl calls cout
  pl("1/3 = ", 1.0/3);

  pl("---------------------------------");

  cout.precision(15); // 15 is enough

  pl("Newton to find root(2)");
  double x = 2; // start
  for (int i = 0; i < 5; i++)
  {
	 x = f(x);
	 p("x = ",x); pl("  x^2 = ",x*x);
  }

  pl("---------------------------------");

  pl("Newton to find ln(2)");
  x = 2; //start
  for (i = 0; i < 7; i++)
  {
	 x = g(x);
	 p("x = ",x); pl("  exp(x) = ",exp(x));
  }

  pl("---------------------------------");
  pl("fiterate(g,1)   = ", fiterate(g,1)); // default n = 3
  pl("fiterate(g,1,7) = ", fiterate(g,1,7));


  }

  /* output ----------------------------------

	--------------
    iterate1.cpp
   --------------
   1/3 = 0.333333333333333315
   ---------------------------------
   Newton to find root(2)
   x = 1.5  x^2 = 2.25
   x = 1.41666666666667  x^2 = 2.00694444444444
   x = 1.41421568627451  x^2 = 2.00000600730488
   x = 1.41421356237469  x^2 = 2.00000000000451
   x = 1.4142135623731  x^2 = 2
   ---------------------------------
   Newton to find ln(2)
   x = 1.27067056647323  exp(x) = 3.56324115146432
   x = 0.831957303739969  exp(x) = 2.29781185737446
   x = 0.702350584017167  exp(x) = 2.01849176999946
   x = 0.693189402250512  exp(x) = 2.00008444516383
   x = 0.693147181451268  exp(x) = 2.00000000178265
   x = 0.693147180559945  exp(x) = 2
   x = 0.693147180559945  exp(x) = 2
   ---------------------------------
   fiterate(g,1)   = 0.693147180560026
   fiterate(g,1,7) = 0.693147180559945


 */