//////////////////////////////////////////////////////////////////////////////
// ratmin.m
//
// Nils Bruin (April 20, 2012)
//
// Magma functions to compute minimal models of rational functions
// under the conjugation action of PGL_2(Q).
//////////////////////////////////////////////////////////////////////////////

intrinsic Resultant(f::FldFunRatUElt[RngInt])->RngIntElt
{Returns the resultant of numerator and denominator of the rational function f.
Note that the resultant is taken by considering numerator and denominator to
be of the same degree. Also note that numerator and denominator are normalized
in the sense that the coefficients are integer and without common factor.}
  d:=Degree(f);
  N:=Numerator(f);
  D:=Denominator(f);
  return LeadingCoefficient(N)^(d-Degree(D))*
         ((-1)^d*LeadingCoefficient(D))^(d-Degree(N))*Resultant(N,D);
end intrinsic;

intrinsic kbounds(f::FldFunRatUElt[RngInt],p::RngIntElt)->RngIntElt,RngIntElt
{Returns lower and upper bounds on k involved in minimizing the p-adic valuation
of the resultant of f.
See Theorems 3.3.2 and 3.3.3 in [Molnar, M.Sc. thesis]}
  F:=Numerator(f);
  x:=Parent(F).1;
  G:=Denominator(f);
  d:=Degree(f);
  dG:=Degree(G);
  
  //Lemma 3.3
  if dG gt (d+1)/2 then
    // k > -2/(2*dG+1)*ord_p(lc(G))
    // to get an inclusive lower bound we take floor and add one.
    lowerk:=1+Floor(-2/(2*dG-d+1)*Valuation(LeadingCoefficient(G),p));
  else
    // k > -2/(d+1) * ord_p(lc(F))
    // to get an inclusive lower bound we take floor and add one.
    lowerk:=1+Floor(-2/(d+1)*Valuation(LeadingCoefficient(F),p));
  end if;

  //Lemma 3.3
  // k < 2*ord_p(Res(F-x*G,G))/(d-1)
  // to get an inclusive upper bound we take ceiling and subtract one.
  upperk:=Ceiling(2*Valuation(Resultant(F-x*G,G),p)/(d-1))-1;
  
  return lowerk,upperk;
end intrinsic;

intrinsic interpolate(c::SeqEnum)->FldFunRatUElt
{Given a sequence of distinct numbers, interpolates a rational function with
that sequence as part of an orbit. The interpolation method is naive, so it may
fail in some degenerate cases.}
  m:=#c;
  d:=Ceiling((m-2)/2);
  M:=Matrix([[c[i]^(j): j in [0..d]] cat
             [-c[i+1]*c[i]^(j): j in [0..d]]: i in [1..m-1]]);
  b:=Basis(Kernel(Transpose(M)));
  b:=Eltseq(b[1]);
  return Polynomial(b[1..d+1])/Polynomial(b[d+2..2*d+2]);
end intrinsic;

function pad(L,b)
  //pad a list with 0 to given length.
  return L cat [0: i in [1..b-#L]];
end function;

intrinsic bequations(phi::FldFunRatUElt[RngInt],p::RngIntElt,k::RngIntElt)->.
{Constructs the coefficient inequalities involved in determining
Resultant reducing affine transformations

Input:
 phi -- a rational function
 p   -- a prime
 k   -- an integer
 
Output:
 L   -- a sequence consisting of tuples <c(b'), e>
 t   -- an exponent
 
We consider affine linear transformations A: x :-> p^k*x + b such that
ord_p(Res(phi^A)) < ord_p(Res(phi)).

If such a transformation exists then b = b' * p^k for some integer b', such that
for each tuple <c(b'),e> we have ord_p(c(b')) > e
}
  d:=Degree(phi);
  x:=Parent(Numerator(phi)).1;
  f:=pad(Eltseq(Numerator(phi)),d+1);
  g:=pad(Eltseq(Denominator(phi)),d+1);
  
  //these are the coefficients of the polynomials G(X+b) and F(X+b)-b*G(X+b)
  //respectively (note that x=b here)
  //MAGMA list are indexed base 1, hence the (i+1)
  
  h:=[&+[Binomial(i,j)*g[i+1]*x^(i-j):i in [j..d]]:j in [0..d]];
  r:=[&+[Binomial(i,j)*(f[i+1]-x*g[i+1])*x^(i-j): i in [j..d]]:j in [0..d]]; 

  //we want that ord_p( F(p^k*X+b)-b*G(p^k*X+b)) and ord_p(p^k*G(p^k*X+b))
  //both are bigger than alpha=(d+1)*k/2
  //That means that we must find x such that
  // ord_p( h[i+1])+k+k*i  > alpha
  // ord_p( r[i+1]*p^(k*i) ) > alpha
  // we make a list of tuples <C(b),e> standing for the inequality ord_p(C(b)) > e
  
  L:=[];
  alpha:=(d+1)*k/2;
  for i in [0..d] do
    Append(~L, <h[i+1],alpha-k*(i+1)>);
    Append(~L, <r[i+1],alpha-k*i>);
  end for;
  
  //for each of these equations, we know that if the valuation of b is low enough, then
  //ord_p(C(b)) is determined by the leading term. We make these bounds explicit and
  //compute for each equation the minimal valuation for b that could possibly lead to
  //a solution (i.e., the minimum for whatever the leading term would lead to together
  //with the bound below which the leading term becomes dominant).
  //See Lemma 3.7
  
  valbbounds:=[];
  for l in L do
    C:=l[1];
    alphaC:=l[2];
    c:=Eltseq(C);
    di:=Degree(C);
    if di gt 0 then
      val_cdi:=Valuation(c[di+1],p);
      m:=Minimum([(alphaC-val_cdi)/di] cat
           [(Valuation(c[j+1],p)-val_cdi)/(di-j):j in [0..di-1]]);
      Append(~valbbounds,m);
    end if;
  end for;
  
  //since each bound must be satisfied, we can take the maximum of them as a lower
  //bound on ord_p(b).

  bminval:=Ceiling(Maximum(valbbounds));
  
  //We rewrite the equations in terms of b' = b/p^bminval
  //so that we know that integral b' should be sufficient.
  //this means substituting b = b' * p^bminval
  //which is straightforward if bminval >=0. Otherwise we have to clear denominators.
  
  if bminval ge 0 then
    Lnew:=[<Evaluate(l[1],p^bminval*x),l[2]> : l in L];
  else
    Lnew:=[];
    for l in Lnew do
      d:=Degree(l[1]);
      C:=Eltseq(l[1]);
      polnew:=Parent(C)![p^(-bminval*(d-i))*C[i+1]: i in [0..d]];
      Append(~Lnew,<polnew,l[2]-bminval*d>);
    end for;
  end if;
  
  //we report the inequalities, together with the change of coordinates applied to b.    
  return L,bminval;
end intrinsic;
  
intrinsic solve_ineq(L::SeqEnum,p::RngIntElt)->BoolElt,RngIntElt
{Determines whether a system of valuation inequalities can be satisfied and determine a
solution.

Input:
  L -- Sequence of tuples <C(x),e>
  p -- prime

Output:
  b -- Boolean signifying if the system has a solution
  X -- if b is true, an integer that provides a solution.

A solution X is an integer such that ord_p(C(X)) > e for each tuple in L.
}

  //first we normalize our equations. We loop through them and remove content
  //and short-circuit some trivial cases.
  
  Lnew:=[];
  for l in L do
    c:=Eltseq(l[1]);
    if #c eq 0 then //0 polynomials have no coefficients and impose no conditions
      continue;
    end if;

    // determine valuation of the content of the polynomial.
    v:=Minimum([Valuation(a,p): a in c]);
    
    // if l[2] < v then any integer solves the equation, so we can drop the constraint.
    // if l[2] >= v we need to take a closer look.
    
    if l[2] ge v then
      if #c eq 1 then //a constant polynomial here represents an unsatisfiable constraint.
        return false,_;
      else
        //otherwise, we remove content from both sides, so that the
        //polynomial has non-zero reduction mod p.
        Append(~Lnew,<Parent(l[1])![a/p^v: a in c],l[2]-v>);
      end if;
    end if;
  end for;
  
  //at this point Lnew consists of non-trivial constraints.
  //meaning all polynomials in Lnew are non-constant and have non-zero reduction mod p
  //and all valuation bounds are non-negative.
  //if there are no constraints, then any integer provides a solution.

  if #Lnew eq 0 then
    return true,0;
  end if;
  
  //otherwise, any solution must be a root mod p of each polynomial.
  //so we reduce all polynomials mod p, take the GCD of them and
  //iterate over the roots in GF(p).

  x:=Parent(L[1][1]).1;
  kX:=PolynomialRing(GF(p));
  for r in Roots(GCD([kX!Eltseq(l[1]):l in Lnew])) do
    
    //we investigate if x=r+p*y leads to any integer solutions y
    //we can do so by recursion: we simply substitute the new x into
    //each equation. Upon recursion, the normalization procedure above will remove
    //the content from each polynomial that the substitution introduces.
    
    newx:=(Integers()!r[1])+p*x;
    bl,sol:=solve_ineq([<Evaluate(l[1],newx),l[2]>:l in Lnew],p);

    //if a solution is found, we shortcut and return the solution for x
    //that results from this
    if bl then
      return true, (Integers()!r[1])+p*sol;
    end if;
  end for;

  //if we arrive here, we have exhausted all possibilities and there is no solution,
  //so we return false.
  return false,_;
end intrinsic;

intrinsic IsAffineMinimal(f::FldFunRatUElt[RngInt],p::RngIntElt)->.
{Determine if f is affine minimal at p and return an affine minimal version
together with the affine transformation relating the two.

Input:
  f -- rational function
  p -- prime
  
Output:
  B   -- boolean; true if f is affine minimal at p; false otherwise
  g   -- rational function affine minimal at p
  <a,b> -- transformation such that g=(f(a*x+b)-b)/a
  
We have ord_q(Res(g))=ord_q(Res(f)) for all primes q != p
}
  x:=Parent(f).1;
  f_minimal:=true;
  a_total:=1;b_total:=0;
  repeat
    klower,kupper:=kbounds(f,p);
    bl:=false;
    for k in [klower..kupper] do
      L,bval:=bequations(f,p,k);
      bl,bsol:=solve_ineq(L,p);
      if bl then
        f_minimal:=false;
        b:=bsol*p^bval;
        a:=p^k;
        b_total:=b_total+a_total*b;
        a_total:=a_total*a;
        f:=(Evaluate(f,a*x+b)-b)/a;
        break;  
      end if;
    end for;
  until not(bl);
  return f_minimal,f,<a_total,b_total>;
end intrinsic;
  
intrinsic IsAffineMinimal(phi::FldFunRatUElt[RngInt])->.
{Determine if f is affine minimal and return an affine minimal version
together with the affine transformation relating the two.

Input:
  f -- rational function
  
Output:
  B   -- boolean; true if f is affine minimal; false otherwise
  g   -- rational function affine minimal
  <a,b> -- transformation such that g=(f(a*x+b)-b)/a
}
  R:=Resultant(phi);
  P:=[p[1]: p in Factorization(R)];
  a:=1;b:=0;
  B:=true;
  for p in P do
    bl,phinew,tr:=IsAffineMinimal(phi,p);
    if not(bl) then
      B:=false;
      b:=b+a*tr[2];
      a:=a*tr[1];
      phi:=phinew;
    end if;
  end for;
  return B,phi,<a,b>;
end intrinsic;

intrinsic ProjectiveMap(f::FldFunRatUElt)->.
{Return the projective map on P^1 corresponding to the rational function f}
  d:=Degree(f);
  R:=BaseRing(Parent(f));
  P1<x,z>:=ProjectiveSpace(R,1);
  CP1:=CoordinateRing(P1);
  fhom:=Evaluate(f,x/z);
  return map<P1->P1|[Numerator(fhom),Denominator(fhom)]>;
end intrinsic;

intrinsic Reduction(f::FldFunRatUElt[RngInt],p::RngIntElt)->.
{Return reduction of f modulo p}
  k:=GF(p);
  P:=PolynomialRing(k);
  return (P!Numerator(f))/(P!Denominator(f));
end intrinsic;

intrinsic Cycles(phi::MapSch[Prj[FldFin], Prj[FldFin]])->.
{Return the cycles of rational points under a rational morphism on a projective space.}

  visited:={}; //keep track of points already looked at
  cycles:={};  //and cycles found
  
  P1:=Domain(phi);
  for Q0 in RationalPoints(P1) do
    Q1:=Q0;
    orbit:={@ @}; //build orbit of Q0
    while true do
      //if Q1 equals a point already in orbit, we have found a cycle.
      if Q1 in orbit then
        i:=Index(orbit,Q1); //find index of Q1 in orbit
        Include(~cycles,Setseq(orbit[i..#orbit])); //that gives us the cycle
        Include(~visited,Q1); //remember we have visited Q1 already
        break;
      //if we have already visited Q1 then we end up in a cycle that we
      //have already found, so we can stop.
      elif Q1 in visited then
        break;
      end if;
      
      //otherwise keep track of Q1
      Include(~visited,Q1);
      Include(~orbit,Q1);
      
      //and step to next point in orbit.
      Q1:=phi(Q1);
    end while;
  end for;
  return cycles;
end intrinsic;
 
intrinsic Multiplier(phi::MapSch[Prj, Prj], cycle::SeqEnum)->.
{Compute the multiplier of phi along a cycle

Input:
  phi --   rational map
  cycle -- periodic orbit

Output:
  multiplier
}

  //check that the cycle is valid
  require #Seqset(cycle) eq #cycle and
      forall{i: i in [1..#cycle-1] |cycle[i+1] eq phi(cycle[i])} and
      cycle[1] eq phi(cycle[#cycle]): "Parameter is not a cycle under phi";

  P1:=Domain(phi);
  R:=BaseRing(P1);
  x:=PolynomialRing(R).1;
  infty:=P1![1,0];

  //multipliers at infinity require a change of coordinates.
  //if infinity is a fixed point, we simply use 1/x
  if cycle eq [infty] then
      f:=L[2]/L[1] where L:=[Evaluate(P,[1,x]): P in DefiningPolynomials(phi)];
      fprime:=Derivative(f);
      return Evaluate(fprime,0);
  end if;
  
  f:=L[1]/L[2] where L:=[Evaluate(P,[x,1]): P in DefiningPolynomials(phi)];
  fprime:=Derivative(f);

  //check if infinity occurs in the cycle.
  i:=Index(cycle,infty);
  if i ne 0 then
  
    //if so, normalize such that cycle is [?,infty,...]
    if i eq 1 then
      cycle:= [cycle[#cycle],infty] cat cycle[2..#cycle-1];
    elif i gt 2 then
      cycle:= cycle[i-1..#cycle] cat cycle[1..i-2];
    end if;

    //we have now normalized the cycle to [Q1,infty,...]
    //we know that f(f(Q1)) is not infinity, so we simply get the
    //multiplier of that bit by taking the derivative of f^2 directly.
    fsq:=Evaluate(f,f);
    fsqprime:=Derivative(fsq);
    multiplier:=Evaluate(fsqprime,cycle[1][1]/cycle[1][2]);
    cycle:=cycle[3..#cycle];
  else
    multiplier:=1;
  end if;
  
  //at this point, cycle does not contain infty, so we can compute the rest
  //of the multiplier via the chain rule.
  for P in cycle do
    multiplier *:=Evaluate(fprime,P[1]/P[2]);
  end for;
  return multiplier;
end intrinsic;

intrinsic RationalPeriods(f::FldFunRatUElt[RngInt]: Primes:=50)->.
{Compute the possible periods of rational points under f based on p-adic properties at several
primes

Input:
  f -- rational map
  Primes -- (optional) A bound on the size of primes used or a set of primes to be used

Output:
  Set of possible periods
  
Only primes at least 3 and of good reduction are used. Furthermore, the smallest good
prime in Primes is also not used, to ensure that the remaining primes will not occur in periods
of rational points.
}
  if ISA(Type(Primes),RngIntElt) then
    Primes:=PrimesInInterval(5,Primes);
  end if;
  require 2 notin Primes: "Cannot use prime 2";
  Sort(~Primes);
  rs:=Resultant(f);
  PossiblePeriods:=false;
  lengths:=false;
  first:=true;
  for p in Primes do
    if rs mod p ne 0 then
      if first then //skip first good prime
        first:=false;
        continue;
      end if;
      phi:=ProjectiveMap(Reduction(f,p));
      C:=Cycles(phi);
      
      //for each cycle in reduction, #c or ord(multiplier)*#c is a possible
      //length of a cycle reducing to it. [this uses the good reduction at p
      //and to remove further powers of p from the consideration we use that p
      //is not the smallest good prime]
      
      lengths:={Integers()|#c: c in C} join
           {#c*Order(m) : c in C | m ne 0 where m:=Multiplier(phi,c)};
      if PossiblePeriods cmpeq false then
        PossiblePeriods:=lengths;
      else
        PossiblePeriods meet:=lengths;
        if #PossiblePeriods eq 0 then
          return PossiblePeriods;
        end if;
      end if;
    end if;
  end for;
  //we intersect the possible periods found for each prime.
  return PossiblePeriods;
end intrinsic;

function myratpoints(V)
  //helper function: Rational points for projective schemes over ZZ does not work
  //but its result should be the same as over QQ.
  if BaseRing(V) eq Integers() then
    Vt:=RationalPoints(BaseChange(V,Rationals()));
    //Magma automatically takes common factors out of coordinates.
    return {@V![Numerator(v[1])*Denominator(v[2]),Numerator(v[2])*Denominator(v[1])]: v in Vt@};
  end if;
  return RationalPoints(V);
end function;

intrinsic RationalPeriodicPoints(phi::MapSch[Prj, Prj], n::RngIntElt)->.
{Determine the rational periodic points of order dividing n}
  L:=DefiningPolynomials(phi^n);
  R:=Universe(L);
  P1:=Domain(phi);
  return myratpoints(Scheme(P1,R.1*L[2]-R.2*L[1]));
end intrinsic;

intrinsic RationalPreperiodicPoints(f::FldFunRatUElt[RngInt])->.
{Determines all rational preperiodic points of a rational map}

  //find possible periods
  P:=RationalPeriods(f);
  phi:=ProjectiveMap(f);
  
  //find rational points of those orders. If P is a large set, one could
  //likely be a bit smarter about this computation, but usually P is very small.
  //we need to investigate all these points for rational preimages.
  investigate:=Set(&join[RationalPeriodicPoints(phi,n): n in P]);

  visited:=investigate;
  while #investigate ne 0 do
    Q:=Rep(investigate); //get a point that needs its preimages determined.
    Exclude(~investigate,Q);
    new:=Set(myratpoints(Q@@phi)); //find its rational preimages
    investigate join:=new diff visited; //and add the ones not yet visited.
    visited join:=new; //make sure these points are not considered in the future.
  end while;
  //when we arrive here, no points with univestigated rational preimages exist.
  
  return visited; //return all points found in the process.
end intrinsic;
