import numpy as np
import math
import random
import numpy.polynomial as nppols

#coefficient de stathme minimal
def min_stathm(A):
    n,p=A.shape
    m=99999999
    res=(0,0)
    for i in range(n):
        for j in range(p):
            a=A[i,j]
            if a!=0 and abs(a)<m:
                m=abs(a)
                res=(i,j)
    return res[0],res[1],m

#pgcd des coefficients d'une matrice
def matrix_gcd(A):
    n,p=A.shape
    res=0
    for i in range(n):
        for j in range(p):
            res=math.gcd(res,A[i,j])
    return res

#matrice aléatoire
def random_matrix(n,p):
    res=np.zeros((n,p),dtype=int)
    for i in range(n):
        for j in range(p):
            res[i,j]=random.randint(-1,1)
    return res

#quotient division euclidienne
def quotient(x,q):
    if q>0:
        return x//q
    else:
        return -(x//(-q))

#échange les lignes
def swap_rows(A,i,j):
    A[i],A[j]=A[j].copy(),A[i].copy()

#ppcm
def lcm(a,b):
    return (a*b)//math.gcd(a,b)

#retourne une équation cartésienne minimale de Im A sous la forme d'une matrice B de rang maximal telle que Im A = Ker B
def cart_of_im(A):
    n,p=A.shape
    C=A.copy()
    I=np.eye(n,n,dtype=int)
    zero_rows = hermite(C,I)
    return I[n-zero_rows:n]

#retourne une base de Ker A sous la forme d'une matrice B de rang maximal telle que Im B = Ker A
def base_of_ker(A):
    return cart_of_im(A.transpose()).transpose()

#retourne une matrice C telle que A=BC lorsque c'est possible (A et B ont autant de lignes)
def factor(A,B):
    n,p=A.shape
    q=B.shape[1]
    AC=A.copy()
    BC=B.copy()
    rank=full_gauss(BC,AC)
    return AC[0:rank]

#calcule les facteurs invariants du quotient Ker Dn / Im Dn+1
def homology_factors(Dn,Dnp1):
    A=base_of_ker(Dn)
    C=factor(Dnp1,A)
    factors=smith(C.astype(int))
    return np.array(factors)

#calcule le quotient des images
def quotient_images(A,B):
    C=factor(B,A)
    factors=smith(C.astype(int))
    factors += [0]*(A.shape[1]-len(factors))
    return np.array(factors)

#effectue des transformations sur les lignes jusqu'à avoir une matrice Jr lorsque c'est possible (peut sortir de Z)
def full_gauss(A,I=None):
    n,p=A.shape
    modify_I= not (I is None)
    rank=gauss(A,I)

    #pas besoin d'aller jusqu'à la première ligne
    for i1 in range(rank-1,0,-1):
        col_pivot=-1
        for j in range(p-1,-1,-1):
            if A[i1,j]!=0:
                col_pivot=j
                #normalisation
                val=A[i1,j]
                A[i1]=A[i1]/val
                if modify_I:
                    I[i1]=I[i1]/val

        for i in range(i1):
            c=A[i,col_pivot]
            if c!=0:
                A[i]=A[i]-c*A[i1]
                if modify_I:
                    I[i]=I[i]-c*I[i1]

    return rank


#effectue des transformations sur les lignes jusqu'à avoir une matrice échelonnée
def gauss(A, I=None):
    n,p=A.shape
    min_row=0
    modify_I= not (I is None)
    for i in range(p):
        found=False
        for j in range(min_row,n):
            c=A[j,i]
            if c!=0:
                if found:
                    #pivotage
                    A[j]=A[j]-c*A[min_row]
                    if modify_I:
                        I[j]=I[j]-c*I[min_row]
                else:
                    found=True
                    A[j]=A[j]/c      #normalisation
                    if modify_I:
                        I[j]=I[j]/c

                    #échange des lignes
                    if min_row!=j:
                        swap_rows(A,min_row,j)
                        if modify_I:
                            swap_rows(I,min_row,j)

        if found:
            min_row+=1

    #retourne le rang
    return min_row

#algorithme de hermite
#modifie A et applique les memes opérations à I, si I est fournie
def hermite(A, I=None):
    n,p=A.shape
    min_row=0
    modify_I= not (I is None)
    for i in range(p):
        done = False
        zero_col = True
        while not done:
            #recherche d'un élément non nul de stathme minimal
            min_index=-1
            min_stat = -1

            for j in range(min_row,n):
                c=A[j,i]
                if c!=0:
                    s=abs(c)
                    if min_index==-1 or s<min_stat:
                        min_index=j
                        min_stat=s

            #colonne à partir de min_row non nulle
            if min_index>=0:
                zero_col=False
                #échange de la ligne pivot avec la première ligne du bloc en cours
                if min_index!=min_row:
                    swap_rows(A,min_index,min_row)
                    if modify_I:
                        swap_rows(I,min_index,min_row)

                pivot=A[min_row,i]
                done=True
                #pivotage
                for j in range(min_row+1,n):
                    c=A[j,i]
                    if c!=0:
                        done=False
                        q=c//pivot
                        A[j]=A[j]-q*A[min_row]
                        if modify_I:
                            I[j]=I[j]-q*I[min_row]
            else:
                done=True

        if not zero_col:
            min_row+=1

    #nombre de lignes nulles finales
    return n-min_row



#algorithme de smith : calcule les facteurs invariants du conoyau de A
def smith(A):
    factors = rec_smith(A)
    factors += [0]*(A.shape[0]-len(factors)) #voir théorie
    return factors

def rec_smith(A):
    n,p=A.shape

    if np.all(A==0):
        return []

    diag=min(n,p)

    if diag == 0:
        return []

    i,j,m=min_stathm(A)

    if m==0:
        return [0]*diag

    d=matrix_gcd(A)

    if d==m:
        #on force la positivité du coefficient de valeur absolue d
        if A[i,j]<0:
            A[i]=-A[i]

        #on fait en sorte que a00=d avec des transvections
        if i>0:
            A[0]=A[0]-(A[0,j]//d-1)*A[i]
            i=0
        if j>0:
            A[:,0]=A[:,0]-(A[i,0]//d-1)*A[:,j]
            j=0

        #on pivote
        for k in range(1,n):
            A[k]=A[k]-(A[k,0]//d)*A[0]
        for k in range(1,p):
            A[:,k]=A[:,k]-(A[0,k]//d)*A[:,0]

        return [d]+rec_smith(A[1:,1:])

    else:
        found_k=-1
        for k in range(n):
            a=A[k,j]
            if a%m!=0:
                found_k=k
                break
        if found_k>=0:
            A[found_k]=A[found_k]- (A[found_k,j]//m)*A[i]
            #on se ramène à une matrice de stathme minimal strictement plus petit
            return rec_smith(A)

        #on fait en sorte d'avoir que des m sur la colonne j
        for k in range(n):
            if k!=i:
                A[k]=A[k]-(A[k,j]//m-1)*A[i]

        #on cherche quelqu'un qui n'est pas divisible par m
        found_i, found_j = -1,-1
        for ii in range(n):
            for jj in range(p):
                if A[ii,jj] % m != 0:
                    found_i,found_j=ii,jj
                    break
        if found_i>=0:
            A[:,found_j]=A[:,found_j]-(A[found_i,found_j]//m)*A[:,j]
            #on se ramène à une matrice de stathme minimal strictement plus petit
            return rec_smith(A)

        print("probleme rec_smith", A)

#affiche les groupes d'homologie
def disp_group(factors):
    facs=list(factors)
    if facs == []:
        return "0"
    res=""
    cur_fac=1
    cur_count=0
    first=True

    #on ajoute -1 à la fin pour que le dernier facteur soit quand même traité
    for f in (facs+[-1]):
        if int(f)==cur_fac:
            cur_count+=1
        else:
            if cur_fac != 1:
                if not first:
                    res+=" + "

                res+="Z"

                if cur_fac>1:
                    res+=str(cur_fac)

                if cur_count>1:
                    res+="^"+str(cur_count)

                first=False

            cur_count=1
            cur_fac=int(f)

    if res=="":
        return "0"

    return res

#chaine qui représente un polynome
def str_pol(p):
    res=""
    coefs = p.coef

    for i, a in enumerate(coefs):
        if int(a) == a:              #enlever la virgule
            a = int(a)
        if i == 0:                   #coefficient constant
            if a > 0:
                res += "{a} + ".format(a=a)
            elif a < 0:              #parenthèses autout des coefficients négatifs
                res += "({a}) + ".format(a=a)

        elif i == 1:                 #coefficient en X
            if a == 1:               #le 1 ne doit pas être affiché
                res += "X + "
            elif a > 0:
                res += "{a} X + ".format(a=a)
            elif a < 0:
                res += "({a}) X + ".format(a=a)
        else:
            if a == 1:
                res += "X**{i} + ".format(i=i)
            elif a > 0:
                res += "{a} X**{i} + ".format(a=a, i=i)
            elif a < 0:
                res += "({a}) X**{i} + ".format(a=a, i=i)
    res = res[:-3]
    if res=="":
        return "0"
    return res

#affiche un polynome
def print_pol(p):
    print(str_pol(p))

#matrice de gram
def gram_matrix(vectors):
    n=len(vectors)
    res = np.zeros((n,n))
    for i in range(n):
        for j in range(i+1):
            res[i,j] = res[j,i] = np.dot(vectors[i],vectors[j])
    return res

#déterminant de la matrice de gram
def det_gram(vectors):
    return np.linalg.det(gram_matrix(vectors))

#procédé d'orthonormalisation de gram
def gram_process(vectors):
    res=[]
    for i in range(len(vectors)):
        u=vectors[i].copy()
        for v in res:
            u = u - np.dot(u,v) * v #pas grave si u change pendant le calcul
        nrm = np.linalg.norm(u)
        if abs(nrm)>=0.00000001:  #permet de gérer aussi les familles liées
                                  #(on oublie complétement les vecteurs combinaison linéaires des autres)
            u=u/nrm
            res.append(u)

    return res










