def init_0 (N):
    "renvoie une liste composée de N zéros"
    return [0 for i in range(N)]

def addition (l_1, l_2): # On suppose qu'elles ont la même longueur
    return [l_1[i]+l_2[i] for i in range(len(l_1))]

def multiplication (liste, N):
    "On veut multiplier tous les termes par N et renvoyer une nouvelle liste contenant ces éléments"
    return [i*N for i in liste]

def matrice_lambda(S,A,q): 
    "A est la matrice du système mahlerien considéré, S la matrice remplie comme dans A-F"
    n = A.nrows()
    K = A[0,0].numerator().parent()
    x = K.gen()
    [M_grand,M_petit,h_grand,h_petit,borne_rumiz,borne_rumiz_petite,borne_AF,borne_AF_petite]=variables(A,q)
    M = borne_rumiz
    #M = 1000
    h = h_grand
    minimum = min(n*(h+1),M+1)
    if S.rank()==minimum:
        print ("toutes les fi sont linéairement indép")
        #return [i for i in [1..n]] #liste des indices des fonctions lin. indép
    Liste = S.left_kernel().basis()
    s = len(Liste)
    #print "s = ", s
    r = n-s
    #print len(Liste[0])
    liste_totale = []
    for element in Liste:
        liste_w = [] #liste contenant les w_i qui sont des vecteurs lignes de longueur n
        for i in range(0,len(element)-n+1,n):
            liste_n_termes = []
            for j in [0..n-1]:
                liste_n_termes.append(element[i+j])
            liste_w.append(liste_n_termes)
        liste_totale.append(liste_w)
    liste_polynomes_w_s = []
    #print liste_totale
    #print liste_totale[0][0]
    #print "len liste_totale =", len(liste_totale)
    #print "len liste_w = ", len (liste_w)
    #print "len len liste element =", len(liste_totale[0])
    #print "h = ", h
    #J = matrix(FractionField(K),1,len(Liste[0]))
    #J[0] = Liste[0]
    #print J*S
    for j in liste_totale:
        #print j[0]
        W = init_0(n)
        for k in [0..h]:
            temp = multiplication(j[k],x^k)
            W = addition(W,temp)
            #print W
        liste_polynomes_w_s.append(W)
    #print liste_polynomes_w_s
    #Construction de la matrice Lambda
    Lambda = matrix(FractionField(K),s,n)
    for i in range(s):
        Lambda[i] = liste_polynomes_w_s[i]
    return Lambda.echelon_form()

def liste_fonctions_indep(Lambda):
    "Idée: tester cette fonction sur une matrice A quelconque puis changer l'ordre des entrées de A "
    "et relancer les algos précédents pour que ça corresponde aux fonctions qu'on a envie d'isoler."
    n = Lambda.ncols()
    s = Lambda.rank()
    Lambda = Lambda[:s]
    r = n-s
    # On sait qu'il y a r = n-s fonctions linéairement indépendantes
    mineurs_s = Lambda.minors(s)
    liste_colonnes_retenues = Combinations(Lambda.ncols(),r).list()
    liste_fcts_indep = []
    for i in range(len(mineurs_s)):
        if mineurs_s[i] != 0:
            liste_fcts_indep.append(liste_colonnes_retenues[i])
    return liste_fcts_indep

def liste_fonctions_independantes(S,A,q):
    "Idée: tester cette fonction sur une matrice A quelconque puis changer l'ordre des entrées de A "
    "et relancer les algos précédents pour que ça corresponde aux fonctions qu'on a envie d'isoler."
    Lambda = matrice_lambda(S,A,q)
    n = Lambda.ncols()
    s = Lambda.rank()
    Lambda = Lambda[:s]
    r = n-s
    # On sait qu'il y a r = n-s fonctions linéairement indépendantes
    mineurs_s = Lambda.minors(s)
    liste_colonnes_retenues = Combinations(Lambda.ncols(),r).list()
    liste_fcts_indep = []
    for i in range(len(mineurs_s)):
        if mineurs_s[i] != 0:
            liste_fcts_indep.append(liste_colonnes_retenues[i])
    return liste_fcts_indep


def systeme_fcts_indep (S, A, q):
    "Renvoie le système vérifié par les r premières fonctions indépendantes, "
    "qui correspondent aux r premières lignes de A" 
    "après avoir lancé la fonction liste_fonctions_independantes au préalable pour les séclectionner"
    n = A.nrows()
    K = A[0,0].numerator().parent()
    Lambda = matrice_lambda(S, A, q) # S dépend de A 
    s = Lambda.rank()
    r = n-s
    #print ("s = ", s)
    new_A = matrix(FractionField(K),n,n)
    liste_premiere_famille_indep = liste_fonctions_indep(Lambda)[0]
    #print "liste_premiere_famille_indep =",liste_premiere_famille_indep
    liste_termes_restants = list(range(n))
    for w in liste_premiere_famille_indep:
        del liste_termes_restants[liste_termes_restants.index(w)]
    #print "liste_termes_restants = ", liste_termes_restants
    for k in range(len(liste_premiere_famille_indep)):
        new_A[k] = A[liste_premiere_famille_indep[k]]
    for l in [0..n-1-r]:
        new_A[r+l] = A[liste_termes_restants[l]]
    #print(new_A)
    matrice_S = matrix(FractionField(K),n,n)
    for i in [0..r-1]:
        matrice_S[i,i] = 1
    for j in [0..s-1]:
        matrice_S[r+j] = Lambda[j]

    #print(matrice_S)
    B = ((matrice_S(x^q))^(-1))
    B = new_A(x)*B
    B = matrice_S(x)*B
    #print((matrice_S(x) * new_A(x) * (matrice_S(x^q)^(-1)))[:r, :r])
    return B[:r,:r]

def equa_inhom_min(B,q): 
#### Implicitement, on se place ici dans le cas où la fct n'est pas rationnelle 
    s = B.nrows()-1
    K = B[0,0].numerator().parent()
    x = K.gen()
    vecteur_0 = vecteur(s+1,s)
    liste_vecteurs = [vecteur_0]
    matrice = matrix.identity(FractionField(K),s+1)
    matrice_inv = matrix.identity(FractionField(K),s+1)
    for l in [0..s-1]:
        temp = B(x^(q^l))
        matrice =  matrice*temp
        matrice_inv = temp^(-1)*matrice_inv
        liste_vecteurs.append(matrice_inv[0])
    C = matrix(FractionField(K),s+1,s+1)
    for i in [0..s]:
        C[i] = liste_vecteurs[i]
    #print C
    r = C.rank()
    #print "le rang de C est ", r
    mbre_droit = matrix(K,1,s+1)
    mbre_droit[0,0] = K(1)
    #print mbre_droit
    u = C[:r,:r].solve_left(mbre_droit) # ceci a bien un sens car le noyau à gauche est de dimension >= 1
    vecteur_u = padding_0(u[0].list(),s+1)
    print("vecteur_u =", vecteur_u)
    liste_u_i = [vecteur_u[i] for i in [0..r-1]]
    p_0 = LCM(liste_u_i[k].denominator() for k in range(r))
    liste_p_i = [K(-p_0*liste_u_i[0]),K(p_0)] + [K(-p_0*liste_u_i[j]) for j in [1..r-1]]
    return liste_p_i

def composition(liste_pol_1,liste_pol_2,b):
    "Le but est de renvoyer L_1oL_2 où les L_i sont des polynômes en M avec M l'opérateur de Mahleur d'indice b"
    K = liste_pol_1[0].parent()
    A.<M> = PolynomialRing(K)
    liste_res = []
    n_1 = len(liste_pol_1)
    n_2 = len(liste_pol_2)
    res = A(0)
    for i in range(n_1):
        for j in range(n_2):
            res += liste_pol_1[i]*M^i*liste_pol_2[j](x^(b^i))*M^j
    return res.list()

def elimine_polynome(P,b):
    K = P.parent()
    x = K.gen()
    A.<M> = PolynomialRing(K)
    return P*M - P(x^b)

def inhom_to_hom(liste_polynomes,b):
    "On part d'une équation inhomogène et on veut en obtenir une équation homogène "
    P = elimine_polynome(-liste_polynomes[0],b)
    K = liste_polynomes[0].parent()
    A = P.parent()
    M = A.gen()
    Q = A(0)
    for i in [1..len(liste_polynomes)-1]:
        Q += liste_polynomes[i]*M^(i-1)
    liste_P = P.list()
    liste_Q = Q.list()
    equation_hom = composition(liste_P,liste_Q,b)
    pgcd = GCD(t for t in equation_hom)
    #print(pgcd)
    K = equation_hom[0].parent()
    equation_hom = [K(t/pgcd) for t in equation_hom]
    return equation_hom

def rajout_fct_1_fin(A): #rajoute la fonction 1 à la fin
    K = A[0,0].parent()
    n = A.ncols()
    new_A = matrix(K, n + 1, n + 1)
    for i in range(n):
        new_A[i] = list(A[i]) + [0]
    new_A[n] = vecteur(n + 1, n)
    return new_A

def rajout_fct_1_debut(A): #rajoute la fonction 1 au début
    K = A[0,0].parent()
    n = A.ncols()
    new_A = matrix(K, n + 1, n + 1)
    new_A[0] = vecteur(n + 1, 0)
    for i in [1..n]:
        new_A[i] = [0] +list(A[i - 1])
    return new_A

def liste_alea(M, nbre_col): 
### Renvoie une liste d'entiers de nbre_col entiers compris entre 0 et M inclus. 
    l = []
    while(len(l)) != nbre_col:
        a = randint(0,M)
        if not a in l:
            l.append(a)
    return l

def liste_alea_log(N, K, b):
    "Renvoie une liste de N entiers compris entre 0 et K en privilégiant les entiers plus petits "
    liste = []
    while(len(liste) < N):
        a = floor(b^(log(uniform(0, K), b)))
        if not a in liste:
            liste.append(a)
    return sorted(liste)



# Remplissage de la nouvelle matrice S, qu'on nommera S_A:

def calcul_rang_alea (b, function, A, nbre_max, last_column, p):
    t_0 = cputime()
    [M_relat,M_lindep,h_relat,h_lindep,borne_rumiz_relat,
    borne_rumiz_lindep,borne_AF_relat,borne_AF_lindep]=variables(A,b)
    print("h_lindep = ",h_lindep)
    print("h_relat = ", h_relat)
    print("M_relat = ", "%e"%M_relat)
    print("M_lindep = ", "%e"%M_lindep)
    print("borne_rumiz_relat = ", "%e"%borne_rumiz_relat)
    print("borne_rumiz_lindep = ", "%e"%borne_rumiz_lindep)
    h = h_relat
    M = M_relat
    n = A.nrows()
    nombre_lignes = n*(h+1)
    print("nombre de lignes = ", nombre_lignes)
    nombre_cols = nombre_lignes # Pour au moins voir le cas de rang plein
    rho = 3
    #print("calcul des variables (en s):", cputime() - t_0)
    
    #### On choisit aléatoirement Borne_max colonnes à remplir :

    liste_indices = liste_alea_log(nbre_max, last_column, b)

    #### initialisation :
    
    rang = 0
    temp = matrix()
    RES = matrix() 
    t_tot = cputime()

    while(nombre_cols <= nbre_max):
        t_ini = cputime(subprocesses = True)
        S_A = matrix(nombre_cols, nombre_lignes)
        S_A[:temp.nrows(), :nombre_lignes] = temp 
        for i in range(temp.nrows(), nombre_cols):
            J = liste_indices[i]
            #print("on remplit la colonne numéro:", "%e"%J)
            ti = cputime(subprocesses = True)
            m = max(J - h -1, -1)
            ligne_i = flatten([function(j, p)  for j in range(J, m, -1)])
            ligne_i = padding_0(ligne_i, n*(h+1))
            S_A[i] = ligne_i
            #print("temps pour remplir cette colonne = ", cputime(subprocesses = True) - ti)
            #print("nombre de colonnes remplies = ", i + 1)
        print("   Temps de construction de la sous-matrice (en s):", cputime(subprocesses = True) - t_ini)

        temp = S_A
        
        #### Calcul du rang :
        
        t_rang = cputime(subprocesses = True)
        rang_temp = S_A.rank()
        print("Temps de calcul du rang (en s):", cputime(subprocesses = True) - t_rang)
        print("Le rang de la sous-matrice de S_A est ", S_A.rank())
        
        # On veut garder la matrice de plus petite taille avec le rang maximal.
        # Rq:  avec notre construction, le rang ne peut que grandir. 
        if rang_temp > rang:
            rang = rang_temp
            RES = S_A.transpose()
        
        #### Vérifs
        
        print("nombre de colonnes = ", S_A.nrows())
        #print(S_A.nrows() == nombre_cols)
        print("Le rang max qu'on a est ", rang_temp)
        
        #### On incrémente :
        
        if nombre_cols*rho > nbre_max:
            print("rang max obtenu = ", rang)
            print("obtenu pour un nombre de colonnes égal à", RES.ncols())
            print("temps total de calcul =", cputime() - t_0)

        nombre_cols *= rho
        
    return RES


def calcul_rang (b, function, A, nbre_max, p):
    t_0 = cputime()
    [M_relat,M_lindep,h_relat,h_lindep,borne_rumiz_relat,
    borne_rumiz_lindep,borne_AF_relat,borne_AF_lindep]=variables(A,b)
    print("h_lindep = ",h_lindep)
    print("h_relat = ", h_relat)
    print("M_relat = ", "%e"%M_relat)
    print("M_lindep = ", "%e"%M_lindep)
    print("borne_rumiz_relat = ", "%e"%borne_rumiz_relat)
    print("borne_rumiz_lindep = ", "%e"%borne_rumiz_lindep)
    h = h_relat
    M = M_relat
    n = A.nrows()
    nombre_lignes = n*(h+1)
    print("nombre de lignes = ", nombre_lignes)
    nombre_cols = nombre_lignes # Pour au moins voir le cas de rang plein
    rho = 3
    #print("calcul des variables (en s):", cputime() - t_0)
    
    # On choisit aléatoirement Borne_max colonnes à remplir :

    liste_indices = range(nbre_max)


    #### initialisation :
    
    rang = 0
    temp = matrix()
    RES = matrix() 
    t_tot = cputime()

    while(nombre_cols < nbre_max):
        t_ini = cputime(subprocesses = True)
        S_A = matrix(nombre_cols, nombre_lignes)
        S_A[:temp.nrows(), :nombre_lignes] = temp 
        for i in range(temp.nrows(), nombre_cols):
            #print("on remplit la colonne numéro:", "%e"%i)
            ti = cputime(subprocesses = True)
            m = max(i - h -1, -1)
            ligne_i = flatten([function(j, p)  for j in range(i, m, -1)])
            ligne_i = padding_0(ligne_i, n*(h+1))
            S_A[i] = ligne_i
            #print("temps pour remplir cette colonne = ", cputime(subprocesses = True) - ti)
            #print("nombre de colonnes remplies = ", i + 1)
        print("   Temps de construction de la sous-matrice (en s):", cputime(subprocesses = True) - t_ini)

        temp = S_A
        
        #### Calcul du rang :
        
        t_rang = cputime(subprocesses = True)
        rang_temp = S_A.rank()
        print("Temps de calcul du rang (en s):", cputime(subprocesses = True) - t_rang)
        print("Le rang de la sous-matrice de S_A est ", S_A.rank())
        
        # On veut garder la matrice de plus petite taille avec le rang maximal.
        # Rq:  avec notre construction, le rang ne peut que grandir. 
        if rang_temp > rang:
            rang = rang_temp
            RES = S_A.transpose()
        
        #### Vérifs
        
        print("nombre de colonnes = ", S_A.nrows())
        #print(S_A.nrows() == nombre_cols)
        print("Le rang max qu'on a est ", rang_temp)
        
        #### On incrémente :
        
        if nombre_cols*rho >= nbre_max:
            print("rang max obtenu = ", rang)
            print("obtenu pour un nombre de colonnes égal à", RES.ncols())
            print("temps total de calcul =", cputime() - t_0)

        nombre_cols *= rho
        
    return RES