/*
*   Modified NEWTON ALGORITHM 2.3
*
*   To find a solution to f(x) = 0 given an
*   initial approximation p0.
*   
*   f(x) = exp(x) - x - 1 has multiplicity 2 at root p = 0.
*
*   INPUT:   initial approximation p0; tolerance TOL;
*            maximum number of iterations NO.
*
*   OUTPUT:  approximate solution p or a message of failure
*/

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

#define true    1
#define false   0
#define PI      3.1415926535897932384626433832795

static   double F(double);
static   double FP(double);
static   double F_doubleP(double);
static   void OUTPUT(int *);

int main()
{
      double TOL,P0,D,F0,FP0, F_doubleP0;
      int OK,I,NO,FLAG;

      P0 = 1.0;
      TOL = 1.0e-13;
      NO = 20;

      OUTPUT(&FLAG);
      F0 = F(P0);

      /* STEP 1 */ 
      I = 1;      
      OK = true;        
      /* STEP 2 */ 
      while ((I<=NO) && true == OK) 
      {  
	 /* STEP 3 */ 
	 /* compute P(I) */ 
	 FP0 = FP(P0);
         F_doubleP0 = F_doubleP(P0);
	 D = F0*FP0/(FP0*FP0 - F0*F_doubleP0);
	 /* STEP 6 */ 
	 P0 = P0 - D;  
	 F0 = F(P0);
         F_doubleP0 = F_doubleP(P0);
	 if (FLAG == 2) printf("%3d   %14.8e   %14.7e\n",I,P0,F0); 
	 /* STEP 4 */ 
	 if (fabs(D) < TOL) { 
	    /* procedure completed successfully */
	    printf("\nNewton's method successed\n");
	    printf("\nApproximate solution = %.10e\n",P0);
	    printf("with F(P) = %.10e\n",F0);
	    printf("Number of iterations = %d\n",I);
	    printf("Tolerance = %.10e\n",TOL);
	    OK = false;
	 }
	 /* STEP 5 */ 
	 else I++;      
      }
      if (true == OK) {
	 /* STEP 7 */ 
	 /* procedure completed unsuccessfully */
	 printf("\nNewton's method failed\n");
	 printf("\nIteration number %d",NO);
	 printf(" gave approximation %.10e\n",P0);
	 printf("with F(P) = %.10e not within tolerance  %.10e\n",F0,TOL);
      }
      return 0;
}

/* Change functions F and FP for a new problem */
static double F(double X)
{
   double f; 

   f = exp(X) - X - 1.0; 
   return f;
}
/*  FP is the derivative of F */
static double FP(double X)
{
   double fp; 

   fp = exp(X) - 1;
   return fp;
}
/* F_doubleP is the second derivative of F*/
static double F_doubleP(double X)
{
   double fp; 

   fp = exp(X);
   return fp;
}

static void OUTPUT(int *FLAG)
{
   char NAME[30];

   printf("Select amount of output\n");
   printf("1. Answer only\n");
   printf("2. All intermeditate approximations\n");
   printf("Enter 1 or 2: ");
   scanf("%d", FLAG);
   printf("Modified Newtons Method\n");
   if (*FLAG == 2) printf("  I                P             F(P)\n");
}