#!/usr/bin/env python3

"""
    auteur: Dylan Voisin - M1 Informatique Luminy
    MNI: TP 4
"""

import numpy as np
from numpy.polynomial import polynomial as P
from copy import deepcopy as dp
from math import isclose

def vect_propre(A, µ=1, nb_iter=15):
    I = np.identity(A.shape[0])
    v = np.ones(A.shape[0])
    inv = np.linalg.inv(A - µ*I)
    for i in range(nb_iter):
        v = inv @ v / np.linalg.norm(inv@v)
    return v

def vect_propre2(A, eps = 0.1):
    v = np.random.random(A.shape[0])
    v /= np.linalg.norm(v)
    lambda_ = v.transpose() @ A @ v
    eps = 1e-13
    écart = float('inf')
    while écart > eps:
        w = lambda_ * A @ v
        v = w / np.linalg.norm(w)
        old = lambda_
        lambda_ = (v.transpose() @ A @ v) / (v.transpose() @ v)
        écart = abs(old - lambda_)
    return v

def vect_propre_généralisé(A, eps = 1-9):
    A_ = A @ A.transpose()
    v = vect_propre2(A_, eps)
    v = (v.transpose() @ A @ v) / (v.transpose() * v)
    return v

def eval_poly(p, x):
    """un polynome sera représenté par une liste pleine de ses coefficients
    f = a_0 + a_1 * x ** 1 + … + a_n * x ** b_n
    [a_0, a_1, a_2, …, a_n]
    """
    return sum(p[i] * x ** i for i in range(len(p)))

def dichotomie(size, step):
    l = [0, size]
    lc = dp(l)
    for i in range(step):
        for j in range(len(l)-1):
            m = (l[j+1] + l[j]) / 2
            lc.insert(2*(j+1)-1, m)
            yield m
        l = dp(lc)

def newton(p, epsilon_f = 1e-13, epsilon_c = 1e-9):
    """peut faire une boucle infinie à cause de l'imprécision
    à régler avec les paramètres epsilon, mais l'imprécision va
    se répercuter sur les valeurs des racines"""
    x = np.random.random() * 100
    pp = np.polyder(p[::-1])[::-1]
    roots = []
    while True:
        # toujours de l'imprécision à gérer
        if abs(eval_poly(p, x)) < epsilon_f:
            roots.append(x)
            p = P.polydiv(p, (-x, 1))[0]
            # on purge l'imprécision
            p = list(map(lambda x: x if abs(x) > epsilon_c else 0, p))
            pp = np.polyder(p[::-1])[::-1]
            if len(p) == 1:
                return roots # f = c
        fpx = eval_poly(pp, x)
        if fpx == 0:
            print(x)
            # extremum (local ou non)
            old = x
            step = 5
            maxrange = 1000
            eps = 1*10**(-8)
            sg = np.sign(eval_poly(pp, x-eps))
            sd = np.sign(eval_poly(pp, x+eps))
            for diff in dichotomie(maxrange, step):
                if np.sign(eval_poly(pp, x-diff)) != sg:
                    x = x - diff
                    break
                elif np.sign(eval_poly(pp, x+diff)) != sd:
                    x = x + diff
                    break
            if x == old:
                return roots # on considère qu'on a toutes les racines
        x -= eval_poly(p, x) / eval_poly(pp, x)
    return roots

if __name__ == '__main__':
    m = np.array(
        [
            [2, 3],
            [3, 2]
        ]
    )
    v = vect_propre(m)
    v /= np.linalg.norm(v)
    print(v)
    v = vect_propre2(m)
    print(v)
    v = vect_propre_généralisé(m)
    v /= np.linalg.norm(v)
    print(v)
    # p = (0, 1, -2) # -2x²+x
    p = (30, -19, -15, 3, 1) # (x+5)(x-3)(x-1)(x+2)
    print(newton(p))