// This program corresponds to section 7.2 of the mvtensor paper with
// the following correspondence for variables:
//
// paper    program
// x_2-x_1  x
// a_1      a(1)
// ...      ...
// a_{10}   a(10)
// b_1      a(11)
// ...      ...
// b_8      a(18)
// b_9      -a(19)
// b_{10}   -a(20)
// a'_1     -b(1)
// a'_2     b(2)
// a'_3     b(3)
// a'_4     -b(4)
// a'_5     -b(5)
// a'_6     b(6)
// a'_7     b(7)
// a'_8     -b(8)
// a'_9     b(9)
// a'_{10}  b(10)
// b'_1     -b(11)
// b'_2     b(12)
// b'_3     b(13)
// b'_4     -b(14)
// b'_5     -b(15)
// b'_6     b(16)
// b'_7     b(17)
// b'_8     -b(18)
// b'_9     -b(19)
// b'_{10}  -b(20)
//
// The program computes two ideals eqnX and eqnY. The ideal \mathfrak p
// in the paper is eqnY + (x_2-x_1). The ideal \mathfrak q is eqnX.
//
// Singular adds redundant generators to eqnX. The last part in the program
// manually removes them: one repeatedly calls
// eqnX = try(i, eqnX);
// for various values of i: if the generator eqnX[i] is superfluous, then
// this commands removes it, otherwise this command leaves eqnX unchanged.
// At the end of this process, when no more simplification can be done,
// eqnX is just the list of 18 generators given in the paper.
//
// To execute the code, launch an interactive Singular session from a
// directory that contains this program and type
// <"codeD4-paper";
// 
// This program runs on Singular 4.0.2.
//
// Adaptations are needed to let it run on subsequent versions of the
// software:
// - The reverse lexicographic order rp (used on line 587 below) has been
// removed in Singular 4.1.0. It can now be simulated with a construct
// M(intmat_expression), but the simplest workaround is to just use the
// usual lexicographic order lp instead.
// - The format of the output of sat() has changed (and the change is not
// documented). The [1] on line 601 should be removed.
//

// INPUTS 
// We are in type D_n.
// Specifically the group is SO(2n) (equal to its own Laglands dual) and
// we focus on tensor products of the natural representation and the adjoint
// representations.
int n = 4;
list ten_type = 2, 2;
int nb_fact = size(ten_type);

// We compute the fusion of two blocks, consisting of the first (split_pos)
// factors and the last (nb_fact)-(split_pos) ones.
int split_pos = 1;

// Description of the target cycle (Y).
// For the adjoint representation, we use Kashiwara and Nakashima's tableau
// (here just column) notation, from top to bottom.
list keyY =   1,  2,
             -2, -1;

// Description of the source cycle (X).
list keyX =  2, -3,
             3, -2;

// CHECKS WHETHER KEYS HAVE THE SAME WEIGHT
int i, u, v;
int s = 0;
for (i=1; i<=nb_fact; i++)
{
  s = s + ten_type[i];
}
intvec wtY, wtX;
for (i=1; i<=n; i++)
{
  wtY[i] = 0;
  wtX[i] = 0;
}
for (i=1; i<=s; i++)
{
  u = keyY[i];
  if (u>0)
  {
    wtY[u] = wtY[u]+1;
  }
  else
  {
    wtY[-u] = wtY[-u]-1;
  }
  v = keyX[i];
  if (v>0)
  {
    wtX[v] = wtX[v]+1;
  }
  else
  {
    wtX[-v] = wtX[-v]-1;
  }
}
for (i=1; i<=n; i++)
{
  if (wtY[i]!=wtX[i])
  {
    print("keyY and keyX have different weights.");
    break;
  }
}

// INITIALIZATION
int flag, j, k, t, u1, u2, v1, v2;
list pardim;
pardim[1] = 0;
for (k=1; k<=nb_fact; k++)
{
  if (ten_type[k]==1)
  {
    pardim[k+1] = pardim[k]+(2*n-2);
  }
  if (ten_type[k]==2)
  {
    pardim[k+1] = pardim[k]+(4*n-6);
  }
}
int totdim = pardim[nb_fact+1];
string str;
ring R = (0,a(1..totdim),x), z, lp;
number b(1..totdim), denom_b(1..totdim), c, d, q;
list vareqY, vareqX;
matrix dia[2*n][2*n];
matrix mat[2*n][2*n];
matrix MAT[2*n][2*n];
matrix MUT[2*n][2*n];
matrix pri[2][2];
matrix prj[2][2];
matrix deu[2][1];
matrix cob[2][1];
matrix coc[1][1];
MAT = 0;
MAT = MAT+1;

// COMPUTATION OF THE TRANSITION MAP BETWEEN THE CHARTS
for (k=1; k<=nb_fact; k++)
{
  s = 0;
  for (i=1; i<k; i++)
  {
    s = s+ten_type[i];
  }
// We first compute x(a(...)) t^{keyY} and multiply MAT on the right by
// this matrix. Here x(a(...)) can be written either as a product of
// one parameter additive subgroups, or as the exponential of the nilpotent
// element \sum a(j) e_{\alpha_j}. We multiply everything by z because
// Singular prefer polynomials.
  MUT = 0;
  MUT = MUT+1;
  if (ten_type[k]==1)
  {
    u = keyY[s+1];
    if (u<0)
    {
      u = (2*n+1)+u;
    }
    dia = 0;
    dia = dia+z;
    dia[u,u] = z2;
    dia[(2*n+1)-u,(2*n+1)-u] = 1;
    mat = 0;
    mat = mat+1;
    t = pardim[k]+1;
    for (i=1; i<=2*n; i++)
    {
      if (i!=u && i!=(2*n+1)-u)
      {
        mat[u,i] = a(t);
        mat[(2*n+1)-i,(2*n+1)-u] = -a(t);
        if (i>u)
        {
          vareqY[size(vareqY)+1] = t;
        }
        t = t+1;
      }
    }
    t = pardim[k];
    q = 0;
    for (i=1; i<=n-1; i++)
    {
      q = q+a(t+i)*a(t+(2*n-1)-i);
    }
    mat[u,(2*n+1)-u] = -q;
    MAT = MAT*mat*dia;
    MUT = MUT*mat*dia;
  }
  if (ten_type[k]==2)
  {
    u1 = keyY[s+1];
    u2 = keyY[s+2];
    if (u1+u2==0)
    {
      flag = 1;
      if (u1==1 || (u1<0 && u1!=-n))
      {
        ERROR("Convention error.");
      }
      if (u1>0)
      {
        u1 = u1-1;
        u2 = (2*n+1)+u2;
      }
      else
      {
        u1=-u1-1;
      }
    }
    else
    {
      flag = 0;
      if (u1<0)
      {
        u1 = (2*n+1)+u1;
      }
      if (u2<0)
      {
        u2 = (2*n+1)+u2;
      }
    }
    mat = 0;
    mat = mat+1;
    t = pardim[k]+1;
    for (i=1; i<=2*n; i++)
    {
      if (i!=u1 && i!=(2*n+1)-u1 && i!=u2 && i!=(2*n+1)-u2)
      {
        mat[u1,i] = a(t);
        mat[(2*n+1)-i,(2*n+1)-u1] = -a(t);
        if (i>u1)
        {
          vareqY[size(vareqY)+1] = t;
        }
        t = t+1;
      }
    }
    t = pardim[k];
    q = 0;
    for (i=1; i<=n-2; i++)
    {
      q = q+a(t+i)*a(t+(2*n-3)-i);
    }
    mat[u1,(2*n+1)-u1] = -q;
    MAT = MAT*mat;
    MUT = MUT*mat;
    mat = 0;
    mat = mat+1;
    t = pardim[k]+(2*n-3);
    for (i=1; i<=2*n; i++)
    {
      if (i!=u1 && i!=(2*n+1)-u1 && i!=u2 && i!=(2*n+1)-u2)
      {
        mat[u2,i] = a(t);
        mat[(2*n+1)-i,(2*n+1)-u2] = -a(t);
        if (i>u2)
        {
          vareqY[size(vareqY)+1] = t;
        }
        t = t+1;
      }
    }
    t = pardim[k]+(2*n-4);
    q = 0;
    for (i=1; i<=n-2; i++)
    {
      q = q+a(t+i)*a(t+(2*n-3)-i);
    }
    mat[u2,(2*n+1)-u2] = -q;
    MAT = MAT*mat;
    MUT = MUT*mat;
    mat = 0;
    mat = mat+1;
    t = pardim[k]+(4*n-7);
    mat[u1,(2*n+1)-u2] = a(t);
    mat[u2,(2*n+1)-u1] = -a(t);
    if (u1+u2<(2*n+1))
    {
      vareqY[size(vareqY)+1] = t;
    }
    MAT = MAT*mat;
    MUT = MUT*mat;
    mat = 0;
    mat = mat+z;
    t = pardim[k]+(4*n-6);
    if (flag)
    {
      mat[(2*n+1)-u2,u1] = a(t);
      mat[(2*n+1)-u1,u2] = -a(t);
      u1 = (2*n+1)-u1;
      u2 = (2*n+1)-u2;
    }
    else
    {
      mat[u1,u1] = z2;
      mat[(2*n+1)-u1,(2*n+1)-u1] = 1;
      mat[u2,u2] = z2;
      mat[(2*n+1)-u2,(2*n+1)-u2] = 1;
      mat[u1,(2*n+1)-u2] = z*a(t);
      mat[u2,(2*n+1)-u1] = -z*a(t);
    }
    if (u1+u2<(2*n+1))
    {
      vareqY[size(vareqY)+1] = t;
    }
    MAT = MAT*mat;
    MUT = MUT*mat;
  }
// We now find the parameters b(pardim[k]+1..pardim[k+1]) such that
// ( x(b(...)) t^{keyX} )^{-1} MAT belong to G(O) and compute this
// product. We save the denominators because they define the open set
// where the transition map between the two charts is actually defined.
  if (ten_type[k]==1)
  {
    v = keyX[k];
    if (v<0)
    {
      v = (2*n+1)+v;
    }
    dia = 0;
    dia = dia+z;
    dia[v,v] = 1;
    dia[(2*n+1)-v,(2*n+1)-v] = z2;
    mat = 0;
    mat = mat+1;
    t = pardim[k];
    for (i=1; i<=2*n; i++)
    {
      if (i!=v && i!=(2*n+1)-v)
      {
        b(t) = number(-subst(MAT[(2*n+1)-i,(2*n+1)-u], z, 0)/subst(MAT[(2*n+1)-v,(2*n+1)-u], z, 0));
        denom_b(t) = denominator(b(t));
        mat[v,i] = -b(t);
        mat[(2*n+1)-i,(2*n+1)-v] = b(t);
        if (i>v)
        {
          vareqX[size(vareqX)+1] = t;
        }
        t = t+1;
      }
    }
    t = pardim[k];
    q = 0;
    for (i=1; i<=n+1; i++)
    {
      q = q+b(t+i)*b(t+(2*n-1)-i);
    }
    mat[v,(2*n+1)-v] = -q;
    if (k==nb_fact)
    {
      break;
    }
    MAT = dia*mat*MAT;
    MAT = MAT/z2;
  }
  if (ten_type[k]==2)
  {
    v1 = keyX[s+1];
    v2 = keyX[s+2];
    if (v1+v2==0)
    {
      flag = 1;
      if (v1==1 || (v1<0 && v1!=-n))
      {
        ERROR("Convention error.");
      }
      if (v1>0)
      {
        v1 = v1-1;
        v2 = (2*n+1)+v2;
      }
      else
      {
        v1=-v1-1;
      }
    }
    else
    {
      flag = 0;
      if (v1<0)
      {
        v1 = (2*n+1)+v1;
      }
      if (v2<0)
      {
        v2 = (2*n+1)+v2;
      }
    }
    pri[1,1] = subst(MAT[(2*n+1)-v1,(2*n+1)-u1], z, 0);
    pri[1,2] = subst(MAT[(2*n+1)-v2,(2*n+1)-u1], z, 0);
    pri[2,1] = subst(MAT[(2*n+1)-v1,(2*n+1)-u2], z, 0);
    pri[2,2] = subst(MAT[(2*n+1)-v2,(2*n+1)-u2], z, 0);
    t = pardim[k]+1;
    for (i=1; i<=2*n; i++)
    {
      if (i!=v1 && i!=(2*n+1)-v1 && i!=v2 && i!=(2*n+1)-v2)
      {
        deu[1,1] = subst(MAT[(2*n+1)-i,(2*n+1)-u1], z, 0);
        deu[2,1] = subst(MAT[(2*n+1)-i,(2*n+1)-u2], z, 0);
        cob = lift(pri, deu);
        b(t) = number(-cob[1,1]);
        b(t+(2*n-4)) = number(-cob[2,1]);
        denom_b(t) = denominator(b(t));
        denom_b(t+(2*n-4)) = denominator(b(t+(2*n-4)));
        t = t+1;
      }
    }
    mat = 0;
    mat = mat+1;
    t = pardim[k]+1;
    for (i=1; i<=2*n; i++)
    {
      if (i!=v1 && i!=(2*n+1)-v1 && i!=v2 && i!=(2*n+1)-v2)
      {
        mat[v1,i] = -b(t);
        mat[(2*n+1)-i,(2*n+1)-v1] = b(t);
        if (i>v1)
        {
          vareqX[size(vareqX)+1] = t;
        }
        t = t+1;
      }
    }
    t = pardim[k];
    q = 0;
    for (i=1; i<=n-2; i++)
    {
      q = q+b(t+i)*b(t+(2*n-3)-i);
    }
    mat[v1,(2*n+1)-v1] = -q;
    MAT = mat*MAT;
    mat = 0;
    mat = mat+1;
    t = pardim[k]+(2*n-3);
    for (i=1; i<=2*n; i++)
    {
      if (i!=v1 && i!=(2*n+1)-v1 && i!=v2 && i!=(2*n+1)-v2)
      {
        mat[v2,i] = -b(t);
        mat[(2*n+1)-i,(2*n+1)-v2] = b(t);
        if (i>v2)
        {
          vareqX[size(vareqX)+1] = t;
        }
        t = t+1;
      }
    }
    t = pardim[k]+(2*n-4);
    q = 0;
    for (i=1; i<=n-2; i++)
    {
      q = q+b(t+i)*b(t+(2*n-3)-i);
    }
    mat[v2,(2*n+1)-v2] = -q;
    MAT = mat*MAT;
    t = pardim[k]+(4*n-7);
    prj[1,1] = subst(MAT[(2*n+1)-v1,(2*n+1)-u1], z, 0);
    prj[2,1] = subst(MAT[(2*n+1)-v1,(2*n+1)-u2], z, 0);
    deu[1,1] = subst(MAT[v2,(2*n+1)-u1], z, 0);
    deu[2,1] = subst(MAT[v2,(2*n+1)-u2], z, 0);
    coc = lift(prj, deu);
    b(t) = number(-coc[1,1]);
    denom_b(t) = denominator(b(t));
    mat = 0;
    mat = mat+1;
    mat[v1,(2*n+1)-v2] = -b(t);
    mat[v2,(2*n+1)-v1] = b(t);
    if (v1+v2<(2*n+1))
    {
      vareqX[size(vareqX)+1] = t;
    }
    MAT = mat*MAT;
    t = pardim[k]+(4*n-6);
    if (flag)
    {
      prj[1,1] = subst(MAT[v2,(2*n+1)-u1]/z, z, 0);
      prj[2,1] = subst(MAT[v2,(2*n+1)-u2]/z, z, 0);
      deu[1,1] = subst(MAT[(2*n+1)-v1,(2*n+1)-u1], z, 0);
      deu[2,1] = subst(MAT[(2*n+1)-v1,(2*n+1)-u2], z, 0);
    }
    else
    {
      prj[1,1] = subst(MAT[(2*n+1)-v1,(2*n+1)-u1], z, 0);
      prj[2,1] = subst(MAT[(2*n+1)-v1,(2*n+1)-u2], z, 0);
      deu[1,1] = subst(MAT[v2,(2*n+1)-u1]/z, z, 0);
      deu[2,1] = subst(MAT[v2,(2*n+1)-u2]/z, z, 0);
    }
    coc = lift(prj, deu);
    b(t) = number(-coc[1,1]);
    denom_b(t) = denominator(b(t));
    mat = 0;
    mat = mat+z;
    if (flag)
    {
      mat[(2*n+1)-v2,v1] = -b(t);
      mat[(2*n+1)-v1,v2] = b(t);
      v1 = (2*n+1)-v1;
      v2 = (2*n+1)-v2;
    }
    else
    {
      mat[v1,v1] = 1;
      mat[(2*n+1)-v1,(2*n+1)-v1] = z2;
      mat[v2,v2] = 1;
      mat[(2*n+1)-v2,(2*n+1)-v2] = z2;
      mat[v1,(2*n+1)-v2] = -z*b(t);
      mat[v2,(2*n+1)-v1] = z*b(t);
    }
    if (v1+v2<(2*n+1))
    {
      vareqX[size(vareqX)+1] = t;
    }
    MAT = mat*MAT;
  }
  MAT = MAT/z2;
  if (k==split_pos)
  {
    MAT = subst(MAT, z, z+x);
  }
}

// PREPARATION FOR RESTART
// The variables a(...) were elements of the base ring, so that we could
// compute rational functions, but now we need to regard them as polynomial
// indeterminates.
link asl = "ASCII:w auxfile";
open(asl);
str = "ideal eqnY = ";
for (i=1; i<=size(vareqY); i++)
{
  str = str+string(a(vareqY[i]))+", ";
}
str = str[1,size(str)-2]+";";
write(asl, str);
str = "ideal pre_eqnX = ";
for (i=1; i<=size(vareqX); i++)
{
  c = numerator(b(vareqX[i]));
  for (j=1; j<=totdim; j++)
  {
    d = numerator(c/denom_b(j));
    while (c!=d)
    {
       c = d;
       d = numerator(c/denom_b(j));
    }
  }
  str = str+string(c)+", ";
}
str = str[1,size(str)-2]+";";
write(asl, str);
str = "poly denom_b(1..totdim);";
write(asl, str);
for (i=1; i<=totdim; i++)
{
  str = "denom_b("+string(i)+") = "+string(denom_b(i))+";";
  write(asl, str);
}
close(asl);
kill asl;

// RESTART
kill R;
ring R = (0), (x,a(1..totdim)), rp;
link asl = "ASCII:r auxfile";
open(asl);
str = read(asl);
close(asl);
execute(str);

// SATURATING
// The ideal eqnX has the same localization than pre_eqnX with respect to
// the multiplicative set generated by the denom_b(1..totdim).
LIB "elim.lib";
ideal eqnX = pre_eqnX;
for(i=1; i<=totdim; i++)
{
  eqnX = sat(eqnX, denom_b(i))[1];
}
print("Ideals eqnY and eqnX have been computed.");
print("You can check whether eqnX is contained in eqnY + (x) by asking division(eqnX, eqnY+ideal(x))[2];");

// REMOVAL OF SUPERFLUOUS GENERATORS
ideal other;
for (i=size(eqnX); i>=1; i--)
  {
    other = ideal();
    for (j=1; j<=size(eqnX); j++)
      {
        if (j!=i)
          {
            other = other + ideal(eqnX[j]);
          }
      }
    if (reduce(eqnX[i], std(other))==0)
      {
        eqnX = other;
      }
  }