///////////////////////////////////////////////////////////////////////////////
version="$Id: intBasis.lib,v 1.0 2010/05/19 Exp$";
category="Commutative Algebra";
info="
LIBRARY: integralBasis.lib   Integral basis in algebraic function field
AUTHORS: Santiago Laplagne,  slaplagn@dm.uba.ar

MAIN PROCEDURES:
 integralBasis(f, vari);     Integral basis of an algebraic function field
";

LIB "normal.lib";

proc integralBasis(poly f, int vari)
"USAGE:  integralBasis(f, vari); f polynomial in two variables, vari integer indicating
         that the vari-th variable of the ring is the integral element
ASSUME:  The basering must be a ring in two variables.
         The polynomial f must be irreducible and monic as polynomial in the
         variable indicated by vari.@*
NOTE:    The procedure might fail or give a wrong output if the assumptions
         do not hold.
RETURN:  a list, say L, of size 2.
@format  L[1] is an ideal I and L[2] is a polynomial D such that the integral basis is
         b_0 = I[1] / D, b_1 = I[2] / D, ..., b_{n-1} = I[n] / D.@*
         That is, the integral closure of k[x] in the algebraic function
         field L(x,y) is @*
         k[x] b_0 + k[x] b_1 + ... + k[x] b_{n-1},@*
         where we assume that x is the trascendental variable, y is the integral
         element (indicated by vari), f gives the integral equation and n is
         the degree of f as a polynomial in y.@*
@end format
THEORY:  We compute the integral basis of the integral closure of k[x] in k(x,y)
         by computing the normalization of the affine ring k[x,y]/<f> and
         converting the k[x,y]-module generators into a k[x]-basis.@*
KEYWORDS: integral basis; normalization.
SEE ALSO: normal.
EXAMPLE: example integralBasis; shows an example
"
{
  int i, j;
  def R = basering;

  // The degree of f with respect to the variable vari
  int n = size(coeffs(f, var(vari))) - 1;

  // If the integral variable is the first, then the universal denominator
  // must be a polynomial in the second variable (and viceversa).
  string conduStr;
  if(vari == 1){
    conduStr = "var2";
  } else {
    conduStr = "var1";
  }

  // We compute the normalization of the affine ring k[x,y]/f(y)
  ideal I = f;
  list nor = normal(I, conduStr);
  ideal normalGen = nor[2][1];
  poly D = normalGen[size(normalGen)];  // The universal denominator

  //Debug information
  // "The denominator is: ", D;
  // "It must be a polynomial in the ", var(3-vari), " variable.";

  // We define a new ring where the integral variable is the first variable
  // (needed for reduction) and has the appropiate ordering.
  list rl = ringlist(R);
  rl[2] = list(var(vari), var(3-vari));
  rl[3] = list(list("C", 0), list("lp", intvec(1,1)));
  def S = ring(rl);
  setring S;

  // We map the elements in the previous ring to the new one
  poly f = imap(R, f);
  ideal normalGen = imap(R, normalGen);

  // We create the system of generatos y^i*f_j.
  list l;
  ideal red = groebner(f);
  for(j = 1; j <= size(normalGen); j++){
    l[j] = reduce(normalGen[j], red);
  }
  for(i = 1; i <= n-1; i++){
    for(j = 1; j <= size(normalGen); j++){
      l[size(l)+1] = reduce(var(1)^i*normalGen[j], red);
    }
  }

  // To eliminate the redundant elements, we look at the polynomials as
  // elements of a free module where the coordinates are the coefficients
  // of the polynomials regarded as polynomials in y.
  // The groebner basis of the module generated by these elements
  // gives the desired basis.
  matrix vecs[n + 1][size(l)];
  matrix coeffi[n + 1][2];

  for(i = 1; i<= size(l); i++){
    coeffi = coeffs(l[i], var(1));
    vecs[1..nrows(coeffi), i] = coeffi[1..nrows(coeffi), 1];
  }
  module M = vecs;
  M = std(M);

  // We go back to the original ring.
  setring R;
  module M = imap(S, M);

  // We go back from the module to the ring in two variables
  ideal G;
  poly g;
  for(i = 1; i <= size(M); i++){
    g = 0;
    for(j = 0; j <= n; j++){
      g = g + M[i][j+1] * var(vari)^j;
    }
    G[i] = g;
  }

  // The first element in the output is the ideal of numerators.
  // The second element is the denominator.
  list outp = G, D;

  return(outp);
}

example
{ "EXAMPLE:";
  printlevel = printlevel+1;
  echo = 2;
  ring s = 0,(x,y),dp;
  poly f = y5-y4x+4y2x2-x4;
  list l = integralBasis(f, 2);
  l;
// The integral basis of the integral closure of Q[x] in Q(x,y) consists
// of the elements of l[1] divided by the polynomial l[2].
  echo = 0;
  printlevel = printlevel-1;
}