/* 
  ニュートン法による方程式の解法
  2005-7-5 Yoshinari Nomura.

  コンパイル方法: cc -lm newton-method-2b.c

  newton-method-2a.c から変更:
   (1) 収束条件を |f(xk)| < delta から |f(xk) - f(xk-1)| < delta に変更
   (2) k = 0 のときを表示
*/

#include <stdio.h>  /* printf を利用するのに必要 */
#include <math.h>   /* 各種算術関数のために必要  */

/* f(x) = 0 の x を求める問題となる関数 */
double f(double x)
{
  return x * x * x - 3 * x - 2;  /* (x+1)(x+1)(x-2) */
}

/* f の一次導関数*/
double f1(double x)
{
  return 3 * x * x - 3;
}

/* ニュートン法の繰り返し関数 */
double newton(double xk)
{
  return xk - (f(xk) / f1(xk));
}

main()
{

  /* 各種初期値設定 */
  int    k    = 0;           /* 繰り返し回数 */
  double ans   = 2;          /* 真の x (本来は不明) */
  double delta = 3E-16;      /* 許容誤差 3E-16 .. 3.0x10の-16乗 */
  double xk    = 10;         /* 初期値 */
  double pxk   = 0;          /* 1つ前の xk を保存する変数 */

  /* 開始メッセージを表示 */
  printf("Newton method program start.\n");

  /* k = 0 のときを表示しておく */
  printf("k:%d xk:%f f(xk):%f ans-xk:%f\n", k, xk, f(xk), ans - xk);

  /* 
     (2) f(xk-1) と f(xk) の差がなくなったら終わり
  */
  while (fabs(pxk - xk) > delta) {
    k = k + 1;
    pxk = xk;           /* pxk に xk を保存しておく */
    xk = newton(xk);    /* 新しい xk をニュートン方で求める */
    printf("k:%d xk:%f f(xk):%f ans-xk:%f\n", k, xk, f(xk), ans - xk);
  }

  /* 終了メッセージを表示 */
  printf("done.\n");
}