clear();
// f(z) = 0
function ret = func(z)
    //(1) sin(z) - 1/2 = 0
   // ret = sin(z) - 1/2;
    //(2) z^3 - 1 = 0
    ret = z^3 - 1;
endfunction

// f'(z)
function ret = dfunc(z)
    //(1)
   // ret = cos(z);
    //(2)
    ret = 3 * z^2;
endfunction

// Newton法
function [ret_z, times] = newton(init_z, maxtimes, reps, aeps)
    old_z = init_z;
    for i = 1:maxtimes
        new_z = old_z - func(old_z) / dfunc(old_z);
        if abs(new_z - old_z) <= abs(old_z) * reps + aeps then
            times = i; // 反復回数
            break;
        end
        old_z = new_z;
//        mprintf("%d, %g + %g * I\n", i, real(old_z), imag(old_z));
//        times = i;
    end
    if i >= maxtimes then
        disp("収束しませんでした。(反復回数: )", maxtimes)
    end
    ret_z = new_z;
endfunction

// Julia集合描画
div = 1000;// Re(=x), Im(=y)軸の分割数
max_abs_xy = 5; // Re, Im範囲の絶対値
julia = zeros(div, div);
xmin = -max_abs_xy, xmax = max_abs_xy;
ymin = -max_abs_xy, ymax = max_abs_xy;
x = linspace(xmin, xmax, div);
y = linspace(ymin, ymax, div);

// 収束チェック
for i = 1:div
    for j = 1:div
        init_z = x(i) + y(j) * %i;
        [z, julia(i, j)] = newton(init_z, 100, 1.0e-5, 1.0e-15)
    end
end

// Julia集合描画
grayplot(x, y, julia)
