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

  コンパイル方法: cc -lm newton-method-1.c
*/

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

/* f(x) = 0 の x を求める問題となる関数 */
double f(double x)
{
  return sin(exp(x));
}

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

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

main()
{

  /* 各種初期値設定 */
  int    n     = 0;                            /* 繰り返し回数 */
  double ans   = log(3.14159265358979323846);  /* 真の x (本来は不明) */
  double delta = 3E-16;                        /* 許容誤差 */
                                               /* 3E-16 .. 3.0x10の-16乗 */
  double xk    = 1;                            /* 初期値 */

  /* 開始メッセージを表示 */
  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");
}