///////////////////////////////////////////////////////////////////////////////
version = "$Id$";
category="Commutative Algebra";
info="
LIBRARY: cimonom.lib   Determines if the toric ideal of an affine monomial curve is a complete intersection

AUTHORS: I.Bermejo,               ibermejo@ull.es
@*       I.Garcia-Marco,          iggarcia@ull.es
@*       J.-J.Salazar-Gonzalez,   jjsalaza@ull.es

OVERVIEW:
  A library for determining if the toric ideal of an affine monomial curve is a complete intersection with NO
  NEED of computing explicitly a system of generators of such ideal. It also contains procedures to obtain the
  minimum positive multiple of an integer which is in a semigroup of positive integers.
  The procedures are based on a paper by Isabel Bermejo, Ignacio Garcia and Juan Jose Salazar-Gonzalez: 'An
  algorithm to check whether the toric ideal of an affine monomial curve is a complete intersection', Preprint.

SEE ALSO: Integer programming

PROCEDURES:
 BelongSemig(n,v[,sup]);    checks whether n is in the semigroup generated by v;
 MinMult(a,b);              computes k, the minimum positive integer such that k*a is in the semigroup of
                            positive integers generated by the elements in b.
 CompInt(d);                checks wether I(d) is a complete intersection or not.
";

LIB "general.lib";

///////////////////////////////////////////////////////////////////////////////////////////////////////////
//
proc BelongSemig(bigint n, intvec v, list #)
"
USAGE:   BelongSemig (n,v[,sup]);  n bigint, v and sup intvec
RETURN:  In the default form, it returns 1 if n is in the semigroup generated by
         the elements of v or 0 otherwise. If the argument sup is added and in case
         n belongs to the semigroup generated by the elements of v, it returns
         a monomial in the variables {x(i) | i in sup} of degree n if we set
         deg(x(sup[j])) = v[j].
ASSUME:  v and sup positive integer vectors of same size, sup has no
         repeated entries, x(i) has to be an indeterminate in the current ring for
         all i in sup.
EXAMPLE: example BelongSemig; shows some examples
"
{
//--------------------------- initialisation ---------------------------------
  int i, j, num;
  bigint PartialSum;
  num = size(v);
  int e = size(#);

  if  (e > 0)
  {
    intvec sup = #[1];
    poly mon;
  }

  for (i = 1; i <= nrows(v); i++)
  {
    if ((n % v[i]) == 0)
    {
      // ---- n is multiple of v[i]
      if (e)
      {
        mon = x(sup[i])^(int(n/v[i]));
        return(mon);
      }
      else
      {
        return (1);
      }
    }
  }

  if (num == 1)
  {
    // ---- num = 1 and n is not multiple of v[1] --> FALSE
    return(0);
  }

  intvec counter;
  counter[num] = 0;
  PartialSum = 0;

  intvec w = sort(v)[1];
  intvec cambio = sort(v)[2];

  // ---- Iterative procedure to determine if n is in the semigroup generated by v
  while (1)
  {
    if (n  >= PartialSum)
    {
      if (((n - PartialSum) % w[1]) == 0)
      {
        // ---- n belongs to the semigroup generated by v,
        if (e)
        {
          // ---- obtain the monomial.
          mon = x(sup[cambio[1]])^(int((n - PartialSum) / w[1]));
          for (j = 2; j <= num; j++)
          {
            mon = mon * x(sup[cambio[j]])^(counter[j]);
          }
          return(mon);
        }
        else
        {
          // ---- returns true.
          return (1);
        }
      }
    }
    i = num;
    while (!defined(end))
    {
      if (i == 1)
      {
        // ---- Stop, n is not in the semigroup
        return(0);
      }
      if (i > 1)
      {
        // counters control
        if (counter[i] >= ((n - PartialSum) / w[i]))
        {
          PartialSum = PartialSum - (counter[i]*w[i]);
          counter[i] = 0;
          i--;
        }
        else
        {
          counter[i] = counter[i] + 1;
          PartialSum = PartialSum + w[i];
          int end;
        }
      }
    }
    kill end;
  }
}
example
{ "EXAMPLE:";
   ring r=0,x(1..5),dp;
   int a = 125;
   intvec v = 13,17,51;
   intvec sup = 2,4,1;
   BelongSemig(a,v,sup);
   BelongSemig(a,v);
}

///////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////

proc MinMult(int a, intvec b)
"
USAGE:   MinMult (a, b); a integer, b integer vector.
RETURN:  an integer k, the minimum positive integer such that ka belongs to the
         semigroup generated by the integers in b.
ASSUME:  a is a positive integer, b is a positive integers vector.
EXAMPLE: example MinMult; shows some examples.
"
{
//--------------------------- initialisation ---------------------------------
  int i, j, min, max;
  int n = nrows(b);

  if (n == 1)
  {
    // ---- trivial case
    return(b[1]/gcd(a,b[1]));
  }

  max = b[1];
  for (i = 2; i <= n; i++)
  {
    if (b[i] > max)
    {
      max = b[i];
    }
  }
  int NumNodes = a + max; //----Number of nodes in the graph

  int dist = 1;
  // ---- Auxiliary structures to obtain the shortest path between the nodes 1 and a+1 of this graph
  intvec queue = 1;
  intvec queue2;

  // ---- Control vector:
  //      control[i] = 0 -> node not reached yet
  //      control[i] = 1 -> node in queue1
  //      control[i] = 2 -> node in queue2
  //      control[i] = 3 -> node already processed
  intvec control;
  control[1] = 3;         // Starting node
  control[a + max] = 0;   // Ending node
  int current = 1;        // Current node
  int next;               // Node connected to corrent by arc (current, next)

  int ElemQueue, ElemQueue2;
  int PosQueue = 1;

  // Algoritmo de Dijkstra
  while (1)
  {
    if (current <= a)
    {
      // ---- current <= a, arcs are (current, current + b[i])
      for (i = 1; i <= n; i++)
      {
        next = current + b[i];
        if (next == a+1)
        {
          kill control;
          kill queue;
          kill queue2;
          return (dist);
        }
        if ((control[next] == 0)||(control[next] == 2))
        {
          control[next] = 1;
          queue = queue, next;
        }
      }
    }
    if (current > a)
    {
      // ---- current > a, the only possible ars is (current, current - a)
      next = current - a;
      if (control[next] == 0)
      {
        control[next] = 2;
        queue2[nrows(queue2) + 1] = next;
      }
    }
    PosQueue++;
    if (PosQueue <= nrows(queue))
    {
      current = queue[PosQueue];
    }
    else
    {
      dist++;
      if (control[a+1] == 2)
      {
        return(dist);
      }
      queue = queue2[2..nrows(queue2)];
      current = queue[1];
      PosQueue = 1;
      queue2 = 0;
    }
    control[current] = 3;
  }
}
example
{ "EXAMPLE:";
   "int a = 46;";
   "intvec b = 13,17,59;";
   "MinMult(a,b);";
   int a = 46;
   intvec b = 13,17,59;
   MinMult(a,b);
   "// 3*a = 8*b[1] + 2*b[2]"
}

///////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////

proc CompInt(intvec d)
"
USAGE:   CompInt(d); d intvec.
RETURN:  1 if the toric ideal I(d) is a complete intersection or 0 otherwise.
ASSUME:  d is a vector of positive integers.
NOTE:    If printlevel > 0, additional info is displayed in case
         I(d) is a complete intersection:
         if printlevel >= 1, it displays a minimal set of generators of the toric
         ideal formed by quasihomogeneous binomials. Moreover, if printlevel >= 2
         and gcd(d) = 1, it also shows the Frobenius number of the semigroup
         generated by the elements in d.
EXAMPLE: example CompInt; shows some examples
"
{
//--------------------------- initialisation ---------------------------------

  int i,j,k,l,divide,equal,possible;

  int n = nrows(d);
  int max = 2*n - 1;
  ring r = 0, x(1..n), dp;

  int level = printlevel - voice + 2;
  // ---- To decide how much extra information calculate and display
  if (level > 1)
  {
    int e = d[1];
    for (i = 2; i <= n; i++)
    {
      e = gcd(e,d[i]);
    }
    if (e <> 1)
    {
      print ("// Semigroup generated by d is not numerical!");
    }
  }
  if (level > 0)
  {
    ideal id;
    vector mon;
    mon[max] = 0;
    if ((level > 1)&&(e == 1))
    {
      bigint frob = 0;
    }
  }

  // ---- Trivial cases: n = 1,2 (it is a complete intersection).
  if (n == 1)
  {
    print("// Ideal is (0)");
    return (1);
  }

  if (n == 2)
  {
    if (level > 0)
    {
      intvec d1 = d[1];
      intvec d2 = d[2];
      int f1 = MinMult(d[1],d2);
      int f2 = MinMult(d[2],d1);
      id = x(1)^(f1) - x(2)^(f2);
      print ("// Toric ideal:");
      id;
      if ((level > 1)&&(e == 1))
      {
        frob = d[1]*f1 - d[1] - d[2];
        print ("// Frobenius number of the numerical semigroup:");
        frob;
      }
    }
    return (1);
  }

  // ---- For n >= 3 (non-trivial cases)
  matrix mat[max][n];
  intvec using, bound, multiple;
  multiple[max] = 0;
  bound[max] = 0;
  using[max] = 0;

  for (i = 1; i <= n; i++)
  {
    using[i] = 1;
    multiple[i] = 0;
    mat[i,i] = 1;
  }
  if (level > 1)
  {
    if (e == 1)
    {
      for (i = 1; i <= n; i++)
      {
        frob = frob - d[i];
      }
    }
  }

  int new, new1, new2;
  for (i = 1; i <= n; i++)
  {
    for (j = 1; j < i; j++)
    {
      if (i <> j)
      {
        new = gcd(d[i],d[j]);
        new1 = d[j]/new;
        new2 = d[i]/new;
        if (!bound[i] ||(new1 < bound[i]))
        {
          bound[i] = new1;
        }
        if (!bound[j] ||(new2 < bound[j]))
        {
          bound[j] = new2;
        }
      }
    }
  }

  // ---- Begins the inductive part
  for (i = 1; i < n; i++)
  {
    // ---- n-1 stages
    for (j = 1; j < n + i; j++)
    {
      if ((using[j])&&(multiple[j] == 0))
      {
        possible = 0;
        for (k = 1; (k < n + i)&&(!possible); k++)
        {
          if ((using[k])&&(k != j)&&(bigint(bound[k])*d[k] == bigint(bound[j])*d[j]))
          {
            possible = 1;
          }
        }
        if (possible)
        {
          // ---- If possible == 1, then c_j has to be computed
          intvec aux;
          // ---- auxiliary vector containing all d[l] in use except d[j]
          k = 1;
          for (l = 1; l < n + i; l++)
          {
            if (using[l] && (l != j))
            {
              aux[k] = d[l];
              k++;
            }
          }

          multiple[j] = MinMult(d[j], aux);
          kill aux;

          if (j <= n)
          {
            if (level > 0)
            {
              mon = mon + (x(j)^multiple[j])*gen(j);
            }
          }
          else
          {
            // ---- if j > n, it has to be checked if c_j belongs to a certain semigroup
            intvec aux, sup;
            k = 1;
            for (l = 1; l <= n; l++)
            {
              if (mat[j, l] <> 0)
              {
                sup[k] = l;
                aux[k] = d[l];
                k++;
              }
            }
            if (level > 0)
            {
              mon = mon + (BelongSemig(bigint(multiple[j])*d[j], aux, sup))*gen(j);
              if (mon[j] == 0)
              {
                // ---- multiple[j]*d[j] does not belong to the semigroup generated by aux,
                // ---- then it is NOT a complete intersection
                return (0);
              }
            }
            else
            {
              if (!BelongSemig(bigint(multiple[j])*d[j], aux))
              {
                // ---- multiple[j]*d[j] does not belong to the semigroup generated by aux,
                // ---- then it is NOT a complete intersection
                return (0);
              }
            }
            kill sup;
            kill aux;
          }

          // ---- Searching if there exist k such that multiple[k]*d[k]= multiple[j]*d[j]
          equal = 0;
          for (k = 1; k < n+i; k++)
          {
            if ((k <> j) && multiple[k] && using[k])
            {
              if (d[j]*bigint(multiple[j]) == d[k]*bigint(multiple[k]))
              {
                // found
                equal = k;
                break;
              }
            }
          }
          // ---- if equal = 0 no coincidence
          if (!equal)
          {
            if (j == n + i - 1)
            {
              // ---- All multiple[k]*d[k] in use are different -> NOT complete intersection
              return (0);
            }
          }
          else
          {
            // ---- Next stage is prepared
            if (level > 0)
            {
              //---- New generator of the toric ideal
              id[i] = mon[j] - mon[equal];
              if ((level > 1)&&(e == 1))
              {
                frob = frob + bigint(multiple[j])*d[j];
              }
            }
            //---- Two exponents are removed and one is added
            using[j] = 0;
            using[equal] = 0;
            using[n + i] = 1;
            d[n + i] = gcd(d[j], d[equal]);  //---- new exponent
            for (l = 1; l <= n; l++)
            {
              mat[n + i, l] = mat[j, l] + mat[equal, l];
            }

            // Bounds are reestablished
            for (l = 1; l < n+i; l++)
            {
              if (using[l])
              {
                divide = gcd(d[l],d[n+i]);
                new = d[n+i] / divide;
                if ((multiple[l])&&(multiple[l] > new))
                {
                  return (0);
                }
                if (new < bound[l])
                {
                  bound[l] = new;
                }
                new = d[l] / divide;
                if  ( !bound[n+i] || (new < bound[n+i]))
                {
                  bound[n+i] = new;
                }
              }
            }
            break;
          }
        }
      }
      if (j == n + i - 1)
      {
        // ---- All multiple[k]*d[k] in use are different -> NOT complete intersection
        return (0);
      }
    }
  }
  if (level > 0)
  {
    "// Toric ideal: ";
    id;
    if ((level > 1)&&(e == 1))
    {
      "// Frobenius number of the numerical semigroup: ";
      frob;
    }
  }
  return(1);
}
example
{ "EXAMPLE:";
   printlevel = 0;
   intvec d = 14,15,10,21;
   CompInt(d);
   printlevel = 3;
   d = 36,54,125,150,225;
   CompInt(d);
   d = 45,70,75,98,147;
   CompInt(d);
};
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////