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

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

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

/* f(x) = 0 の x を求める問題となる関数 typ: 関数の種類  */
double f(int typ, double x)
{
  if (typ == 1) return sin(exp(x));
  if (typ == 2) return x * x * x - 3 * x - 2;  /* (x+1)(x+1)(x-2) */

  return 0; /* xxx 本来は error */
}

/* f の一次導関数 typ: 関数の種類 */
double f1(int typ, double x)
{
  if (typ == 1) return exp(x) * cos(exp(x));
  if (typ == 2) return 3 * x * x - 3;

  return 0; /* xxx 本来は error */
}

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

main()
{
  /* 各種初期値設定 */
  int    k     = 0;           /* 繰り返し回数 */
  double delta = 3E-16;       /* 許容誤差 */
  double xk    = 1.3;         /* 初期値 */
  int    typ   = 0;

  printf("function type: ");
  scanf("%d", &typ);

  /* 開始メッセージを表示 */
  printf("Newton method program start function type = %d.\n", typ);

  while (fabs(f(typ, xk)) > delta){ /* 収束条件 : f(x) が delta 以下になる */
    k = k + 1;
    xk = newton(typ, xk);
    printf("k:%d xk:%f f(xk):%f\n", k, xk, f(typ, xk));

    if (fabs(f(typ, xk)) < delta) break;
  }

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