/* Laguerre methods */
/* laguerre.c */
#include <stdio.h>
#include <math.h>
#define EPS 2e-15
#define NMAX 20
#define N 4
double f(double x);
double df(double x);
double ddf(double x);

int main(void)
{
  int n;
  double x,fx,dfx,ddfx;
  double y,fy,dfy,ddfy;  

  x = 1000.0; /* 右側の初期値 */
  y = 1000.0; /* 左側の初期値 */
  printf("K(x)=%20.15f\n",1-(1.0*N/(N-1))*f(x)*ddf(x)/(df(x)*df(x)));
  n = 0;
  fx = f(x);
  dfx = df(x);
  ddfx = ddf(x);
  fy = f(y);
  dfy = df(y);
  ddfy = ddf(y);    
  printf("%3s %20s %20s\n","n","Laguerre+","Laguerre-");
  printf("%3d %20.15f %10.3e %20.15f %10.3e\n",n,x,fx,y,fy);
  while((fabs(fx) > EPS || fabs(fy) > EPS) && n < NMAX){
    x = x - N * fx /
      (dfx + sqrt((N-1) * (N-1) * dfx * dfx - N * (N-1) * fx * ddfx));
    y = y - N * fy /
      (dfy - sqrt((N-1) * (N-1) * dfy * dfy - N * (N-1) * fy * ddfy));   
    n++;
    fx = f(x);
    dfx = df(x);
    ddfx = ddf(x);
    fy = f(y);
    dfy = df(y);
    ddfy = ddf(y);        
    printf("%3d %20.15f %10.3e %20.15f %10.3e\n",n,x,fx,y,fy);
  }

  return 0;
}

double f(double x)
{
  double y;

  y = 64 + x * (-120 + x * (70 + x * (-15 + x)));
  
  return y;
}

double df(double x)
{
  double y;

  y = -120 + x * (140 + x * (-45 + 4 * x));
    
  return y;
}

double ddf(double x)
{
  double y;

  y = 140 + x * (-90 + x * 12);
    
  return y;
}