/*
 * complex-newton.C -- Newton 法で方程式 f(x)=0 を解く
 *  コンパイル: g++ -o complex-newton complex-newton.C
 *  実行: ./complex-newton
 */

#include <iostream.h>
#include <complex.h>

int main()
{
    int i, maxitr = 100;
    complex<double> f(complex<double>), dfdz(complex<double>), x, dx;
    double eps;

    cout << " 初期値 x0, 許容精度ε=";
    cin >> x >> eps;

    cout.precision(16);
    for (i = 0; i < maxitr; i++) {
        dx = - f(x) / dfdz(x);
        x += dx;
        cout << "f(" << x << ")=" << f(x) << endl;
        if (abs(dx) <= eps) break;
    }
    return 0;
}

/*
 * この関数の例は杉原・室田『数値計算法の数理』岩波書店, p.67 による
 *   f(z) = z^3-2z+2
 *   1実根(≒-1.769292354238631),2虚根(≒0.8846461771193157±0.5897428050222054)
 *   を持つが、原点の近くに初期値を取ると、0 と 1 の間を振動する。
 */

complex<double> f(complex<double> z)
{
  /* return z * z * z - 2 * z + 2; */
  return z * (z * z - 2) + 2;
}

/* 関数 f の導関数 (df/dx のつもりで名前をつけた) */
complex<double> dfdz(complex<double> z)
{
  return 3 * z * z - 2;
}