##################################################################
# AffineMinimal.sage                                             #
#                                                                #
# Alex Molnar (May 22, 2012)                                     #
#                                                                #
# Sage functions to compute minimal models of rational functions #
# under the conjugation action of PGL_2(Q).                      #
##################################################################


R.<z> = PolynomialRing(QQ)
z = R.0


def Resultant(F,G,d):
#Compute the rational function resultant of the polynomials F and G. Namely, we treat both numerator and denominator as being degree max(deg(F),deg(G)).
    dF = R(F).degree()
    dG = R(G).degree()
    if dF < dG:
        return abs(G[d]^(dG-dF)*F.resultant(G))
    else:
        return abs(F[d]^(dF-dG)*F.resultant(G))


def normalize(F,G):
#'Normalizes' a rational function to have coprime integer coefficients.

    H = (F/G).factor()
    H = H.expand()
    F = H.numerator()
    G = H.denominator()
    N = lcm(F.denominator(),G.denominator())
    nF = N*F
    nG = N*G
    return nF,nG


def rational_degenerate(H,d):
#Test to check that F/G has no common factor, and either F or G has degree greater than 1 to remove degenrate cases. See Theorem 2.2.6 in [Molnar, M.Sc. thesis].

    nd = R(H.numerator()).degree()
    dd = R(H.denominator()).degree()
    if (nd == d or dd == d) and d > 1:
        return false
    else:
        return true


def bCheck(c,v,p,b):
#Compute a lower bound on the value of b needed, for a transformation A(z) = z*p^k + b to satisfy ord_p(Res(phi^A)) < ord_p(Res(phi)) for a rational map phi. See Theorem 3.3.5 in [Molnar, M.Sc. thesis].

    #Input:
    # c - a list of polynomials in b. See v for their use.
    # v - a list of rational numbers, where we are considering the inequalities ord_p(c[i]) > v[i].
    # p - a prime.
    # b - local variable.

    #Output:
    # bval - lower bound in Theorem 3.3.5

    val = floor(v+1)
    S.<b> = PolynomialRing(QQ)
    deg = S(c).degree()
    coeffs = S(c).coeffs()
    lcoeff = coeffs[deg]; coeffs.remove(lcoeff)
    check1 = [(coeffs[i].valuation(p) - lcoeff.valuation(p))/(deg - i) for i in range(0,len(coeffs)) if coeffs[i] <> 0]
    check2 = (val - lcoeff.valuation(p))/deg
    check1.append(check2)
    bval = min(check1)
    return ceil(bval)


def scale(c,v,p):
#Given an integral polynomial c, we can write c = p^i*c', where p does not divide c'. We return c' and v - i.

    scaleval = min([coeff.valuation(p) for coeff in c.coefficients()])
    if scaleval > 0:
        c = c/(p^scaleval)
        v = v - scaleval
    if v <= 0:
        flag = false
    else:
        flag = true
    return [flag,c,v]


def blift(LF,Li,p):
#Search for a solution to the given list of inequalities. If found, lift the solution to an appropriate valuation. See Lemma 3.3.6 in [Molnar, M.Sc. thesis]

    P.<b> = PolynomialRing(QQ)
    #Determine which inequalities are trivial, and scale the rest, so that we only lift as many times as needed.
    keepScaledIneqs = [scale(P(coeff),Li,p) for coeff in LF if coeff <> 0]
    keptVals = [i[2] for i in keepScaledIneqs if i[0]]
    if keptVals <> []:
        #Determine the valuation to lift until.
        liftval = max(keptVals)
    else:
        #All inequalities are satisfied.
        return true,1
    S.<b> = PolynomialRing(Integers(p))
    keptScaledIneqs = [S(i[1]) for i in keepScaledIneqs if i[0]]
    #We need a solution for each polynomial on the left hand side of the inequalities, so we need only find a solution for their gcd.
    g = gcd(keptScaledIneqs)
    rts = g.roots(multiplicities=False)
    for r in rts:
        #Recursively try to lift each root
        newInput = P(r) + P(p)*P(b)
        LG = [P(F)(newInput) for F in LF]
        lift,lifted = blift(LG,Li,p)
        if lift:
            #Lift successful.
            return true,P(r) + P(p)*P(lifted)
    #Lift non successful.
    return False,0


def Affine_minimal(F,G, D=[],prnt = false,quick = false):
#Given polynomials F and G, this procedure determines if the rational map F/G is minimal. In particular, it determines if the map is affine minimal, which is enough to decide if it is minimal or not. See Theorem 3.1.3 in [Molnar, M.Sc. thesis].

    #Input: 
    # F,G - polynomials in z.
    # D - a list of primes, in case one only wants to check minimality at those specific primes.
    # prnt - a boolean value, if true all transformations applied to minimized F/G are printed.
    # quick - a boolean value, if true the algorithm terminates once algorithm determines F/G is not minimal, otherwise algorithm only terminates once a minimal model has been found.

    #Output:
    # flag - boolean; true if F/G is affine minimal, and false otherwise.
    # minF,minG - polynomials such thta minF/minG is a minimal model of F/G
    
    #flag will be true until the algorithm finds an affine transformation showing F/G is not minimal.
    flag = true
    minF,minG = F,G
    d = max(R(minF).degree(),R(minG).degree())
    minF,minG = normalize(minF,minG)
    if rational_degenerate(minF/minG,d):
        if prnt:
            #Affine minimality is only considered for normalized maps of degree higher than 2, not of the form f or 1/f for a polynomial f.
            print "The rational function is degenerate."
        return (False,0,0)
    else:
        #If the valuation of a prime in the resultant is small enough, we can say the map is affine minimal at that prime without using the local minimality loop. See Theorem 3.2.2 in [Molnar, M.Sc. thesis]
        if d%2 == 0:
            g = d
        else:
            g = 2*d
        #Calculate the resultant of minF/minG
        Res = Resultant(R(minF),R(minG),d)

        #Some quantites needed for the local minimization loop, but we compute now since the value is constant, so we do not wish to compute in every local loop. See Theorem 3.3.3 in [Molnar, M.Sc thesis]
        H = R(minF) - z*R(minG)
        md = max(R(H).degree(),R(minG).degree())
        ubRes = Resultant(R(F) - R(z)*R(minG),R(minG),md)

        #Set the primes to check minimality at, if not already prescribed
        if D == []:
            D = prime_divisors(Res)

        #Check minimality at all primes in D. If D is all primes dividing Res(minF/minG), this is enough to show whether minF/minG is minimal or not. See Propositions 3.2.1 and 3.3.7 in [Molnar, M.Sc. thesis].
        for p in D:
            while true:
                if Res.valuation(p) < g:
                    #The model is minimal at p
                    min = true
                else:
                    #The model may not be minimal at p.
                    min,minF,minG = Min(R(minF),R(minG),p,d,ubRes,prnt)

                if min:
                    #The model is minimal at p
                    break
                elif F == minF and G == minG:
                    #The model is minimal at p
                    break
                else:
                    #The model is not minimal at p
                    flag = false
                    if quick:
                        break
            if quick and not flag:
                break
        return (flag,minF,minG)


def Min(F,G,p,d,ubRes,prnt):
#Local loop for Affine_minimal, where we check minimality at the prime p.

    #Input:
    # F,G - polynomials.
    # p - a prime.
    # d - the degree of F/G
    # ubRes - The upper bound needed for Theorem 3.3.3 in [Molnar, M.Sc. thesis].
    # prnt - a boolean value, if true all transformations applied to minimize F/G are printed.

    #Output:
    # boolean, true if F/G is minimal at p, false otherwise
    # F,G - if F/G is minimal at p, or
    # minF,minG - polynomials such that minF/minG is local at p

    #First we bound the possible k in our transformations A = zp^k + b. See Theorems 3.3.2 and 3.3.3 in [Molnar, M.Sc. thesis].
    dG = R(G).degree()
    if dG > (d+1)/2:
        lowerBound = floor(-2*(G[dG]).valuation(p)/(2*dG - d + 1) + 1)
    else:
        lowerBound = floor(-2*(F[d]).valuation(p)/(d-1) + 1)
    upperBound = 2*(ubRes.valuation(p))

    if upperBound < lowerBound:
        #There are no possible transformations to reduce the resultant.
        return true,F,G
    else:
        #Looping over each possible k, we search for transformations to reduce the resultant of F/G
        k = lowerBound
        while k <= upperBound:
            Q.<z,b> = PolynomialRing(QQ)
            b = Q.1
            A = p^k*z + b
            Ft = Q(R(F)(A) - b*R(G)(A))
            Gt = Q(p^k*R(G)(A))
            Fcoeffs = [Ft.coefficient({z:i}) for i in range(0,Ft.degree(z)+1)]
            Gcoeffs = [Gt.coefficient({z:i}) for i in range(0,Gt.degree(z)+1)]
            coeffs = Fcoeffs + Gcoeffs
            RHS = (d + 1)*k/2
            #If there is some b such that Res(phi^A) < Res(phi), we must have ord_p(c) > RHS for each c in coeffs.

            #Make sure constant coefficients in coeffs satisfy the inequality.
            coeffsdeg = [Q(i).degree() for i in coeffs]
            satisfy = [QQ(coeffs[i]).valuation(p) > RHS for i in range(0,len(coeffs)) if coeffsdeg[i] == 0]

            if min(satisfy):
                #Constant coefficients in coeffs have large enough valuation, so check the rest.
                #We start by checking if simply picking b=0 works
                btrivial = [QQ(coeff.subs(b=0)) for coeff in coeffs]
                bbool = [coeff.valuation(p) > RHS for coeff in btrivial]
                if min(bbool):
                    #A = z*p^k satisfies the inequalities, and F/G is not minimal
                    Fmin = R(Q(Ft).subs(b=0))
                    Gmin = R(Q(Gt).subs(b=0))
                    Fmin,Gmin = normalize(Fmin,Gmin)  
                    if prnt:
                        print "Conjugating by", p,"^", k, "*z +", 0
                    return false,Fmin,Gmin
                
                #Otherwise we search if any value of b will work. We start by finding a minimum bound on the valuation of b that is necessary. See Theorem 3.3.5 in [Molnar, M.Sc. thesis].
                bval = max([bCheck(coeff,RHS,p,b) for coeff in coeffs if coeff.degree() > 0])
                
                #We scale the coefficients in coeffs, so that we may assume ord_p(b) is at least 0
                scaledCoeffs = [coeff.subs(b = b*p^(bval)) for coeff in coeffs]

                #We now scale the inequalities, ord_p(coeff) > RHS, so that coeff is in ZZ[b]
                scale = QQ(max([coeff.denominator() for coeff in scaledCoeffs]))
                normalizedCoeffs = [coeff*scale for coeff in scaledCoeffs]
                scaleRHS = RHS + scale.valuation(p)

                #We now search for integers that satisfy the inequality ord_p(coeff) > RHS. See Lemma 3.3.6 in [Molnar, M.Sc. thesis].
                bound = floor(scaleRHS+1)
                bool,sol = blift(normalizedCoeffs,bound,p)

                #If bool is true after lifting, we have a solution b, and F/G is not minimal.
                if bool:
                    #Rescale, conjugate and return new map
                    bsol = QQ(R(sol)*R(p^(bval)))
                    Ftn = Ft.subs(b=bsol)
                    Gtn = Gt.subs(b=bsol)
                    Ftn,Gtn = normalize(R(Ftn),R(Gtn))
                    if prnt:
                        print "Conjugating by", p,"^", k, "*z +", bsol
                    return false,Ftn,Gtn
            k = k + 1
    
        return true,F,G


def IntOrb(L,d):
#Given a list of integers, we construct a rational function of degree d, that maps the list of integers, under iteration, in an orbit. See Theorem 4.1.1 in [Molnar, M.Sc. thesis] and the discussion right after.
    F = 0
    G = 0
    M = []
    if d < floor((len(L) - 1)/2 - 1):
        #It is only possible to prescribe 2*d+2 points in an orbit of a rational map of degree d.
        return (F,G)
    IO = MatrixSpace(ZZ, 2*d + 2)
    for i in range(0,len(L)-1):
        M.append([L[i]^(d-k) for k in range(0,d+1)] + [L[i+1]*L[i]^(d-k) for k in range(0,d+1)])
    #pad M with rows of zeros if needed
    n = 2*d + 2 - len(M)
    for i in range(0,n):
        M.append([0 for k in range(0,2*d+2)])
    M = IO(M)
    K = M.transpose().kernel()
    #return the rational function as F and G
    if K <> []:
        F = R([K.gens()[0][d-i] for i in range(0,d+1)])
        G = -R([K.gens()[0][d-i] for i in range(d+1,2*d+2)])
    return (F,G)