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

/* Demonstration of bisection and Newton's method for optimization
   (maximization) of a smooth function in 1-dimension. */


/* Objective function: x / (1 + log(x^2)).  On (-infinity,-1) it has
   one maximum, at x = sqrt(e). */
double fun(const double x)
{
  return x / (1 + log(x*x));
}


/* The first derivative of fun. */
double der1(const double x)
{
  double u = log(x*x);
  double v = 1+u;
  return (u-1) / v*v;
}


/* The second derivative of fun. */
double der2(const double x)
{
  double u = log(x*x);
  return 2*(3 - u) / ((1+u)*x);
}


int main(int argc, 
	 char** argv)
{
  double x1, x2, x3, f1, f2, f3, E, ee=0, xx, ff, ww, g1, g2, e1;
  int it;

  /* The exact value. */
  E = -sqrt(exp(1));

  /*           */
  /* Bisection */
  /*           */
  
  /* This is a bracket. */
  x1 = -5, x2 = -2, x3 = -1;
  f1 = fun(x1), f2 = fun(x2), f3 = fun(x3);
  
  while (x3-x1 > 1e-10)
    {
      if (x2-x1 > x3-x2)
	{
	  xx = (x1+x2)/2;
	  ff = fun(xx);

	  if (ff > f2) x3 = x2, f3 = f2, x2 = xx, f2 = ff;
	  else x1 = xx, f1 = ff;
	}

      else
	{
	  xx = (x2+x3)/2;
	  ff = fun(xx);

	  if (ff > f2) x1 = x2, f1 = f2, x2 = xx, f2 = ff;
	  else x3 = xx, f3 = ff;
	}

      /* The current log_10 error. */
      ee = log(fabs(E-x2))/log(10);

      printf("%20.18f  %18.15f  %18.15f\n", fabs(E-x2), ee, ww-log(x3-x1));

      ww = log(x3-x1);
    }

  /*                 */
  /* Newton's method */
  /*                 */

  printf("\n\n");

  xx = -2;
  for (it=0; it<10; ++it)
    {
      /* Take one Newton step. */
      g1 = der1(xx);
      g2 = der2(xx);
      xx -= g1/g2;

      /* The current log_10 error. */
      ee = log(fabs(E-xx))/log(10);

      printf("%20.18f  %18.15f  %18.5f\n", fabs(E-xx), ee, ee-2*e1);

      e1 = ee;
    }

  return(0);
}