Chebysav method Implementation in C

chebysav.c

/******* Chebysav method is like approimating the given Transcedental Equation into a quadratic equation f(x) = 0, f(x) ~ a0 + a1x + a2x^2

Let xk be an approximate root
f'(x) = a1 + a2x
f''(x) = 2a2  by subsituting the value xk in all the equations we get the values of f(xk), f'(xk), f''(xk)

we get,
      f(x) ~ fk + (x-xk)f'_k + (x-xk)^2*f''_k/2 ==> 0
      x_(k+1) = xk - fk/f'k - (fk^2 f''_k)/2((f'_k)^3) *********/

/*************** PROGRAM STARTS HERE ***************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/********* FUNCTION DECLARATION *********/
float function(float val);
float function_d(float val);
float function_d_d(float val);
void check_bound(float a);
void chebyshev(float a);

/********** MAIN STARTS HERE *********/
int main(int argc, char **argv)
{
   float         a;

   if (argc != 2)
   {
      fprintf(stderr, "Usage: %s <guess>\n", argv[0]);
      exit(1);
   }

   a = atof(argv[1]);

   check_bound(a);

   printf("By using chebyshev method: \n");
   printf("The equation is: \n");
   printf("\t");
   printf("f(x) = x^4 - 26x^2 + 49x - 25\n");

   printf("---------------------------\n");
   printf("   f(x)    f'(x)    f''(x) \n");
   printf("---------------------------\n");

   chebyshev(a);

   exit(0);
}

/********* FUNCTION DECLARATION *********/
void chebyshev(float a)
{
   float        fa, fa_d, fa_d_d, root, i = 0;

   while (1)  //Infinte Loop
   {
      fa = function(a);
      fa_d = function_d(a);
      fa_d_d = function_d_d(a);
      printf("%f  %f  %f\n", fa, fa_d, fa_d_d);

      root = a - (fa/fa_d) - (fa*fa*fa_d_d)/(2*fa_d*fa_d*fa_d);
      
      if (floor(a*10000) == floor(root*10000)) //Comparing the roots
      {
         printf("The given root is %f after %.1f iterations\n", root, i);
         break;  //Getting out of the loop
      }
      a = root;
      i++;
   }

   return ;
}

float function(float val)
{
   float        fx, x = val;

   fx = (x * x * x * x) - (26 * x * x) + (49 * x) - 25;

   return fx;
}

float function_d(float val)
{
   float        fx_d, x = val;

   fx_d = (4 * x * x * x) - (52 * x) + 49;
   return fx_d;
}

float function_d_d(float val)
{
   float        fx_d_d, x = val;

   fx_d_d = (12 * x * x) - 52;
   return fx_d_d;
}

void check_bound(float a)
{
   float        fa;

   fa = function(a);

   if (fa == 0)
   {
      printf("The root is one of the boundaries.\n");
      exit(0);
   }

   return ;
}