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

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

  newton-method.c の解く方程式と初期値を変更しただけ
*/

#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    n     = 0;            /* 繰り返し回数 */
  double ans   = 2;            /* 真の x (本来は不明) */
  double delta = 3E-16;        /* 許容誤差 3E-16 .. 3.0x10の-16乗 */
  double xk    = 10;           /* 初期値 */

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

  /* ニュートン法の繰り返し部分 */
  /* xxx: 本来は n =0 の初期状態から表示すべき
     xxx: 収束の条件判定には工夫が必要
       (1) f(x) が delta 以下になったら終了
       (2) f(xk-1) と f(xk) の差がなくなったら終
       (3) n が ある回数を越えたら終了
       等々
     xxx: 初期値の決め方とその理由
  */

  /* fabs(x) … x の絶対値を求める関数*/
  while (fabs(f(xk)) > delta){ /* 収束条件 : f(x) が delta 以下になる */
    n = n + 1;
    xk = newton(xk);
    printf("n:%d xk:%f f(xk):%f ans-xk:%f\n", n, xk, f(xk), ans - xk);
  }

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