"""Метод Ньютона с центральным Якобианом и уменьшением длины шага.

Неизвестные и невязки предварительно приводятся к согласованным масштабам.
step — шаг разностной производной; tolerance — допуск максимальной невязки.
Пример аналитический: полевой решатель FEMM не вызывается.
"""

import math
import numpy as np


def real_newton(residual, start, step=1e-4, tolerance=1e-7, maximum=20):
    """Найти корень квадратной вещественной системы с конечными входами."""
    if not math.isfinite(step) or step <= 0:
        raise ValueError("Шаг производной должен быть конечным и положительным")
    if not math.isfinite(tolerance) or tolerance <= 0:
        raise ValueError("Допуск должен быть конечным и положительным")
    if type(maximum) is not int or maximum <= 0:
        raise ValueError("Число итераций должно быть положительным целым")

    # Явно создаём вещественный массив: целый start не должен округлять ±step.
    initial = [float(value) for value in start]
    if not initial or any(not math.isfinite(value) for value in initial):
        raise ValueError("Начальное приближение должно содержать конечные числа")
    x = np.array(initial)
    dimension = len(initial)

    def evaluate(state):
        """Проверить размер и конечность невязки перед численным действием."""
        values = [float(value) for value in residual(state)]
        if len(values) != dimension or any(not math.isfinite(v) for v in values):
            raise ValueError("Невязка должна быть конечной и совпадать с размером x")
        return np.array(values)

    for iteration in range(maximum):
        r = evaluate(x)
        norm = max(abs(float(value)) for value in r)
        # Выводим номер итерации и максимальную нормированную невязку.
        print(iteration, norm)
        if norm <= tolerance:
            return x

        # Столбец k содержит производные всех уравнений по неизвестной x[k].
        columns = []
        for k in range(dimension):
            plus = np.array(x)
            minus = np.array(x)
            plus[k] += step
            minus[k] -= step
            columns.append((evaluate(plus) - evaluate(minus)) / (2 * step))
        jacobian = np.array(
            [[columns[k][j] for k in range(dimension)] for j in range(dimension)]
        )
        if any(not math.isfinite(float(v)) for column in columns for v in column):
            raise ValueError("Разностная матрица Якоби содержит неконечное значение")

        # Решаем J·Δx = −r; обратную матрицу явно не вычисляем.
        change = np.linalg.solve(jacobian, -r)
        if any(not math.isfinite(float(v)) for v in change):
            raise ValueError("Поправка Ньютона содержит неконечное значение")
        damping = 1.0
        accepted = False
        for trial in range(12):
            candidate = x + damping * change
            nc = max(abs(float(v)) for v in evaluate(candidate))
            if nc < norm:
                x = candidate
                accepted = True
                break
            damping *= 0.5
        # Без уменьшения невязки метод завершается явным сообщением об отказе.
        if not accepted:
            raise ValueError("Демпфирование не уменьшило невязку")

    # Последняя разрешённая поправка тоже может дать требуемую точность.
    if max(abs(float(v)) for v in evaluate(x)) <= tolerance:
        return x
    raise ValueError("Число итераций исчерпано")


def control(x):
    """Два независимых кубических уравнения с монотонными левыми частями."""
    return [x[0] + 0.2 * x[0] ** 3 - 2.0, x[1] + 0.1 * x[1] ** 3 - 1.0]


# Обе составляющие невязки должны стремиться к нулю; это не расчёт поля.
print(real_newton(control, [1.0, 0.5]))