////////////////////////////////////////////////////////////////////
// ntbuch.lib : norm-trace codes with Buchberger                  //
////////////////////////////////////////////////////////////////////


//LIB "general.lib";
// LIB "matrix.lib";
// LIB "brnoeth.lib";


proc permute_intvec2 (intvec L,intvec P)
{
  int s=size(L);
  int n=size(P);
  intvec Lord=L;
  if (s<n)
  {
    Lord[n]=0;
    s=size(Lord);
  }
  int i;
  for (i=1;i<=n;i=i+1)
  {
    Lord[i]=L[P[i]];
  }
  return(Lord);
}


proc NTcurve2 (int p,int r)
{
  // check correctness of parameters 
  if (r<2)
  {
    ERROR("field extension must have degree 2 or higher");
  }
  if (prime(p)<p)
  {
    ERROR("field characteristic must be prime");
  }
  
  // check if basering is defined and preserve it 
  int defring=defined(basering);
  if (defring>0)
  {
    def base_r=basering;
  }

  // local variables 
  int i,j,s;
  int ino,iyes;
  int li,k;
  intvec lno,lyes;
  intvec ell,t;
  intvec perm;
  intvec Boxes,aux;
  list L;
  list S;

  // general parameters 
  int n=p^(2*r-1);
  int k1=p^(r-1);
  int k2=(p^r)-1;
  int k3=k2-k1;
  int theta=k2/(p-1);

  // Boxes 
  for (i=1;i<=k1;i=i+1)
  {
    Boxes[i]=k2;
  }
  int k4=(k1-1)/(p-1);
  for (i=1;i<=k4;i=i+1)
  {
    aux[i]=p-1;
  }
  Boxes=Boxes,aux;
  Boxes=Boxes,1;
  
  // compute li's 
  ring rT=(p^r,a),(T),lp;
  poly P=T;
  for (i=1;i<r;i=i+1)
  {
    P=P+T^(p^i);
  }
  for (i=0;i<(p^r)-1;i=i+1)
  {
    if (subst(P,T,a^i)==0) 
    {
      iyes=iyes+1;
      lyes[iyes]=i;
    }
    else 
    {
      ino=ino+1;
      lno[ino]=i;
    }
  }

  // compute the ti's 
  poly pi;
  number b;
  for (i=1;i<=k3+1;i=i+1)
  {
    li=lno[i];
    b=a^li;
    pi=T^(theta)-b;
    for (j=1;j<r;j=j+1)
    {
      b=b^p;
      pi=pi-b;
    }
    for (j=0;j<(p^r)-1;j=j+1)
    {
      if (subst(pi,T,a^j)==0)
      {
        break;
      }
    }
    t[i]=j;
  }

  // sort t & permute lno 
  S=sort(t);
  perm=S[2];
  t=S[1];
  lno=permute_intvec2(lno,perm);
  ell=lno,lyes;

  if (lyes<>0)
  {
    s=size(lyes);
  }

  // auxiliary ring 
  ring Rxy=(p^r,a),(x,y),lp;
  // just in case this is useful ... 
  map tau=basering,a*x,(a^theta)*y;
  export(tau);
  // equation of the norm-trace curve 
  poly eq=x^(theta);
  for (i=0;i<r;i=i+1)
  {
    eq=eq-y^(p^(i));
  }
  export(eq);

  // list of points ranged by orbits 
  list O;
  matrix A[k2][2];
  for (i=1;i<=k1;i=i+1)
  {
    A[1,1]=number(a^(t[i]));
    A[1,2]=number(a^(lno[i]));
    for (j=2;j<=k2;j=j+1)
    {
      A[j,1]=number(a)*A[j-1,1];
      A[j,2]=number(a^(theta))*A[j-1,2];
    }
    O[i]=A;
  }
  kill A;
  matrix A[p-1][2];
  for (j=1;j<p;j=j+1)
  {
    A[j,1]=number(0);
  }
  for (i=1;i<=k4;i=i+1)
  {
    A[1,2]=number(a^(lyes[i]));
    for (j=2;j<p;j=j+1)
    {
      A[j,2]=number(a^(theta))*A[j-1,2];
    }
    O[k1+i]=A;
  }
  kill A;
  matrix A[1][2];
  A[1,1]=number(0);
  A[1,2]=number(0);
  O[theta+1]=A;
  export(O);

  if (defring>0)
  {
    setring base_r;
  }
  
  kill rT;

  int g=(p^(r-1)-1)*((p^r-1)/(p-1)-1)/2;
  intvec codepars=p,r,n,0,g,theta,s;

  L[1]=ell;
  L[2]=t;
  L[3]=Boxes;
  L[4]=Rxy;
  // L[5]=DC;
  // L[6]=Rt;
  L[7]=codepars;

  print("n="+string(n));
  print("g="+string(g));
  
  return(L);
}


proc NTmodule_basis (list L,int m)
{
  // check if basering is defined and preserve it 
  int defring=defined(basering);
  if (defring>0)
  {
    def base_r=basering;
  }
  
  list CODE;

  intvec Boxes=L[3];
  int theta=size(Boxes)-1;
  def Rxy=L[4];
  setring Rxy;
  int p=char(basering);
  int N=theta*(p-1);
  N=N+1;
  int r=1;
  while (p<N)
  {
    N=N/p;
    r=r+1;
  }

  // check m 
  int n=p^(2*r-1);
  int g=(p^(r-1)-1)*((p^r-1)/(p-1)-1)/2;
  if (m>=n)
  {
    ERROR("m is not smaller than the number of rational points");
  }
  if (m<=2*g-2)
  {
    ERROR("m is not larger than the Euler characteristic");
  }
  
  // s=#B 
  int s=L[7][7];
  
  int k1=p^(r-1);
  int k2=(p^r)-1;
  int k3=k2-k1;
  // int theta=k2/(p-1);
  int k4=(k1-1)/(p-1);
  
  // list auxL;
  int i,j,k,d,h,counter;
  poly f,ff;

  ideal b;
  counter=1;
  for (i=0;i<theta;i=i+1)
  {
    d=m-i*p^(r-1);
    if (d<0)
    {
      d=0;
    }
    else 
    {
      d=d/theta;
    }

    // s==1 & h==0..d 

    ff=x^i;
    for (j=0;j<=d;j=j+1)
    {
      f=ff*(y^j);
      b[counter]=f;
      counter=counter+1;
    }
  }
  export(b);

  // transform b into a module basis (section 2.3) 
  ring Rt=(p^r,a),t,(c,lp);
  module MB;
  int exponumber;
  counter=counter-1;
  for (i=1;i<=counter;i=i+1)
  {
    vector V;
    for (j=1;j<=theta+1;j=j+1)
    {
      poly H;
      for (k=1;k<=Boxes[j];k=k+1)
      {
        setring Rxy;
        number eva=number(subst(b[i],x,O[j][k,1],y,O[j][k,2]));
        exponumber=pardeg(eva);
        kill eva;
        setring Rt;
        H=H+(a^exponumber)*t^(k-1);
      }
      V=V+H*gen(j);
      kill H;
    }
    MB[i]=V;
    kill V;
  }
  export(MB);
  
  if (defring>0)
  {
    setring base_r;
  }
  
  intvec codepars=p,r,n,m,g,theta,s;
  
  CODE[1]=L[1];
  CODE[2]=L[2];
  CODE[3]=Boxes;
  CODE[4]=Rxy;
  //CODE[5]=DC;
  CODE[6]=Rt;
  CODE[7]=codepars;
  
  return(CODE);
}


proc GB_NTbuchberger (list L)
{
  // requires: NTcode data (DC, rational points, p -or q-, r, theta, m, rings ... 
  // (Evaluate) --> algorithm 4.2 
  // return groebner basis G (Matrix) 
  
  // check if basering is defined and preserve it 
  int defring=defined(basering);
  if (defring>0)
  {
    def base_r=basering;
  }

  list CODE=L;
  
  def Rt=L[6];
  setring Rt;
  module GB=groebner(MB);
  export(GB);
  
  if (defring>0)
  {
    setring base_r;
  }
  
  // update CODE 
  CODE[6]=Rt;

  return(CODE);
}

;