import numpy as np

def Kuengjoe_S12_Aufg5(input_matrix: np.ndarray,
                       initial_vector: np.ndarray,
                       tolerance: float = 1e-4,
                       max_iterations: int = 1000):
    current_normalized_vector = initial_vector.astype(float)
    current_normalized_vector = current_normalized_vector / np.linalg.norm(current_normalized_vector, ord=2)

    number_of_performed_iterations = 0
    last_difference_norm = None

    for iteration_index in range(max_iterations):
        multiplied_vector = input_matrix @ current_normalized_vector
        next_normalized_vector = multiplied_vector / np.linalg.norm(multiplied_vector, ord=2)

        difference_norm = np.linalg.norm(next_normalized_vector - current_normalized_vector, ord=2)
        number_of_performed_iterations = iteration_index + 1
        last_difference_norm = difference_norm

        if difference_norm < tolerance:
            current_normalized_vector = next_normalized_vector
            break

        current_normalized_vector = next_normalized_vector

    rayleigh_quotient_eigenvalue_estimate = float(
        current_normalized_vector.T @ (input_matrix @ current_normalized_vector)
    )

    return rayleigh_quotient_eigenvalue_estimate, current_normalized_vector, number_of_performed_iterations, last_difference_norm


if __name__ == "__main__":
    matrix_a = np.array([
        [1,  1, 0],
        [3, -1, 2],
        [2, -1, 3]
    ], dtype=float)

    initial_vector_x0 = np.array([1, 0, 0], dtype=float)

    dominant_eigenvalue_estimate, dominant_eigenvector_estimate, iteration_count, final_difference_norm = (
        Kuengjoe_S12_Aufg5(matrix_a, initial_vector_x0, tolerance=1e-4)
    )

    print("von-Mises result")
    print("iterations =", iteration_count)
    print("final ||x_{k+1}-x_k||2 =", final_difference_norm)
    print("dominant eigenvalue (Rayleigh) =", dominant_eigenvalue_estimate)
    print("dominant eigenvector =", dominant_eigenvector_estimate)

    eigenvalues_numpy, eigenvectors_numpy = np.linalg.eig(matrix_a)
    index_of_dominant_eigenvalue = int(np.argmax(np.abs(eigenvalues_numpy)))

    dominant_eigenvalue_numpy = float(eigenvalues_numpy[index_of_dominant_eigenvalue])
    dominant_eigenvector_numpy = eigenvectors_numpy[:, index_of_dominant_eigenvalue].astype(float)
    dominant_eigenvector_numpy = dominant_eigenvector_numpy / np.linalg.norm(dominant_eigenvector_numpy, ord=2)

    print("\nnp.linalg.eig check")
    print("dominant eigenvalue (eig) =", dominant_eigenvalue_numpy)
    print("dominant eigenvector (eig) =", dominant_eigenvector_numpy)

    # confronto robusto (segno +/-)
    difference_same_sign = np.linalg.norm(dominant_eigenvector_estimate - dominant_eigenvector_numpy, ord=2)
    difference_opposite_sign = np.linalg.norm(dominant_eigenvector_estimate + dominant_eigenvector_numpy, ord=2)
    print("\nvector agreement (min with +/- sign) =", min(difference_same_sign, difference_opposite_sign))