# Учебная редакция: назначение, единицы и этапы алгоритма пояснены по-русски.
# Исходная программа остаётся самостоятельным свидетельством прежней проверки.
# Основа — II03-3.py: центральные разности Якоби и метод Ньютона с уменьшением шага.
# Аналитические невязки для действующих токов: два комплексных тока, четыре неизвестных.
import math
import numpy as np


def real_newton(residual, start, step=1e-4, tolerance=1e-7, maximum=20, max_calls=300):
    """Решить четыре действительных уравнения методом Ньютона с уменьшением шага.

    Начальные неизвестные переводятся в float, поэтому целочисленный вход
    не обнуляет малые возмущения при вычислении Якоби. step и tolerance —
    положительные конечные числа. maximum ограничивает число поправок,
    max_calls — число вычислений residual, включая пробные состояния.
    """
    if type(maximum) is not int or maximum <= 0:
        raise ValueError("Нужно положительное целое число итераций")
    if type(max_calls) is not int or max_calls <= 0:
        raise ValueError("Нужен положительный целый бюджет вызовов")
    step, tolerance = float(step), float(tolerance)
    if not math.isfinite(step) or step <= 0:
        raise ValueError("Шаг производной должен быть конечным и положительным")
    if not math.isfinite(tolerance) or tolerance <= 0:
        raise ValueError("Допуск должен быть конечным и положительным")
    start_values = [float(value) for value in start]
    if len(start_values) != 4 or any(not math.isfinite(value) for value in start_values):
        raise ValueError("Нужны четыре конечных действительных начальных значения")
    x = np.array(start_values)
    calls = [0]

    def evaluate(state):
        """Учесть вызов и вернуть четыре конечных действительных невязки."""
        if calls[0] >= max_calls:
            raise ValueError("Исчерпан бюджет вычислений невязки или поля")
        calls[0] += 1
        values = [float(value) for value in residual(state)]
        if len(values) != 4 or any(not math.isfinite(value) for value in values):
            raise ValueError("Нужны четыре конечных действительных невязки")
        return np.array(values)

    # Сохраняем невязку принятого состояния: дорогой полевой расчёт не повторяется.
    r = evaluate(x)
    history = []
    for iteration in range(maximum):
        norm = max(abs(float(value)) for value in r)
        history.append(norm)
        if norm <= tolerance:
            return [float(value) for value in x], norm, calls[0], history
        # Четыре столбца Якоби требуют восьми независимых оценок невязки.
        columns = []
        for k in range(4):
            plus, minus = np.array(x), np.array(x)
            plus[k] += step
            minus[k] -= step
            columns.append((evaluate(plus)-evaluate(minus))/(2*step))
        jacobian = np.array([[columns[k][row] for k in range(4)] for row in range(4)])
        # Решаем J·Δx = −r без явного обращения матрицы.
        change = np.linalg.solve(jacobian, -r)
        if any(not math.isfinite(float(value)) for value in change):
            raise ValueError("Линейная система дала неконечную поправку")
        damping, accepted = 1.0, False
        for trial in range(12):
            candidate = x + damping*change
            candidate_r = evaluate(candidate)
            if max(abs(float(value)) for value in candidate_r) < norm:
                x, r, accepted = candidate, candidate_r, True
                break
            damping *= 0.5
        if not accepted:
            raise ValueError("Уменьшение шага не снизило невязку")
    # Последняя разрешённая поправка тоже может достичь заданного допуска.
    final_norm = max(abs(float(value)) for value in r)
    history.append(final_norm)
    if final_norm <= tolerance:
        return [float(value) for value in x], final_norm, calls[0], history
    raise ValueError("Исчерпано число итераций")


def port_residual(x):
    # Как и в учебнике, оба тока направлены внутрь соответствующих обмоток.
    """Вернуть четыре невязки напряжений двух связанных обмоток, В.
    
    x = [Re I1, Im I1, Re I2, Im I2], токи — действующие значения, А.
    Оба тока положительны при входе в обмотку. Источник 40 В относится
    к первой действительной составляющей; нагрузка второй цепи включена явно."""
    a, b, c, d = x
    omega, L1, L2, M = 2*math.pi*50, .12, .03, .05
    # Зависимость активного сопротивления от квадрата действующего тока — учебная модель.
    nonlinear_R1 = .4 + .2*(a*a+b*b)
    Rload, Xload, R2 = 20., 8., .1
    return [nonlinear_R1*a-omega*(L1*b+M*d)-40,
            nonlinear_R1*b+omega*(L1*a+M*c),
            (R2+Rload)*c-Xload*d-omega*(M*b+L2*d),
            (R2+Rload)*d+Xload*c+omega*(M*a+L2*c)]


def run_checks():
    """Проверить уменьшение невязки, достигнутый допуск и отказы при исчерпании бюджета."""
    x, norm, calls, history = real_newton(port_residual, [0., 0., 0., 0.])
    assert norm < 1e-7 and all(history[k+1] < history[k] for k in range(len(history)-1))
    rejected = 0
    for kwargs in [{"maximum": 1}, {"max_calls": 1}]:
        try:
            real_newton(port_residual, [0., 0., 0., 0.], **kwargs)
        except ValueError:
            rejected += 1
    assert rejected == 2
    return {"order": "Re I1,Im I1,Re I2,Im I2; A RMS", "solution": x, "max_port_residual_V": norm, "calls": calls, "budget_failures_rejected": rejected}