////////////////////////////////////////////////////////////////////
// ntcodes.lib : a library to implement the norm-trace codes      //
////////////////////////////////////////////////////////////////////


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


proc permute_intvec (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 iiiiv (int n,intvec v)
{
  // if n is in increasingly sorted v, returns position, else 0 
  int s=size(v);
  int i,r;
  for (i=1;i<=s;i=i+1)
  {
    if (n<v[i])
	{
	  break;
	}
	else 
	{
	  if (n==v[i])
	  {
	    r=i;
		break;
	  }
	}
  }
  return(r);
}


proc delentryintvec (int i,intvec v)
{
  // deletes from v the i-th component 
  int s=size(v);
  intvec w=v;
  if (i==1)
  {
    if (s>1)
    {
      w=v[2..s];
    }
    else 
    {
      w=0;
    }
  }
  if (i>1)
  {
    if (i<s)
	{
	  w=v[1..i-1],v[i+1..s];
	}
  }
  if (i==s)
  {
    if (s>1)
    {
      w=v[1..s-1];
    }
    else 
    {
      w=0;
    }
  }
  return(w);
}


proc NTcurve (int p,int r)
{
  // returns the lists li and ti and a ring with the rational points ranged by orbits, tau and auxiliar functions ... 
  
  // 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;
  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_intvec(lno,perm);
  ell=lno,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);

  // ideal M 
  ideal M;
  for (i=1;i<=theta;i=i+1)
  {
    M[i]=y^(p-1)-number(a^(ell[i]*(p-1)));
  }
  M[theta+1]=y^(p-1);
  export(M);
  
  // ideal F 
  ideal F;
  F[1]=1;
  for (i=2;i<=theta+1;i=i+1)
  {
    F[i]=F[i-1]*M[i-1];
  }
  export(F);
  
  // matrix B -- j from 0 in the paper ... 
  matrix B[theta][p^r];
  for (i=1;i<=p^(r-1);i=i+1)
  {
    for (j=0;j<(p^r)-1;j=j+1)
	{
	  B[i,j+1]=1;
	  for (k=1;k<p-1;k=k+1)
	  {
	    B[i,j+1]=B[i,j+1]*(y-number(a^(ell[i]+(j+k)*theta)));
	  }
	  for (k=1;k<theta;k=k+1)
	  {
	    B[i,j+1]=B[i,j+1]*(x-number(a^(t[i]+j+k*(p-1))));
	  }
    }
  }
  for (i=p^(r-1)+1;i<=theta;i=i+1)
  {
    for (j=0;j<(p^r)-1;j=j+1)
	{
	  B[i,j+1]=1;
	  for (k=1;k<p-1;k=k+1)
	  {
	    B[i,j+1]=B[i,j+1]*(y-number(a^(ell[i]+(j+k)*theta)));
	  }
    }
  }
  export(B);
  
  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;

  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 root_diagram_NTcode (list L,int m)
{
  // input m --> list DC --> poly P(t) --> ideal f --> Evaluate[i,j] 
  // checks, defs, and compute parameters ... 
  // return L with m,DC,[newring,P(t)],f ( --> algorithm 4.2 ) 

  // 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");
  }

  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 DC;
  intvec DCi,DCiNO;
  int i,j,k,s;
  int beta,gamma,counter;
  number alpha;
  for (i=1;i<=theta+1;i=i+1)
  {
    DCiNO=0;
	if (m>=(i-1)*(p^r-1))
	{
	  if (m>=(i-1)*(p^r-1)+(p-2)*theta+(p^(r-1))*(theta-1))
	  {
	    // empty 
		DCi=-1;
	  }
	  else 
	  {
	    // theorem 3.5 
		counter=0;
		for (beta=0;beta<theta;beta=beta+1)
		{
		  for (gamma=0;gamma<=p-2;gamma=gamma+1)
		  {
		    if (m>=(i-1)*(p^r-1)+gamma*theta+beta*p^(r-1))
			{
			  // no additional condition 
			  counter=counter+1;
			  alpha=a^(-beta-gamma*theta);
			  DCiNO[counter]=pardeg(alpha);
			}
		  }
		}
		DCiNO=sort(DCiNO)[1];
		// DCiNO=compress_intvec(DCiNO);
		// print("E("+string(i)+")={"+string(DCiNO)+"}");
		s=size(DCiNO);
		// DCi ... 
	    if (i<=k1)
	    {
	      DCi=0..Boxes[i]-1;
	    }
	    else 
	    {
	      if (i<=theta)
		  {
		    DCi=0;
		    for (j=2;j<=Boxes[i];j=j+1)
		    {
		      DCi=DCi,(j-1)*theta;
		    }
		  }
		  else 
		  {
		    DCi=0;
		  }
	    }
		for (j=1;j<=s;j=j+1)
		{
		  DCi=delentryintvec(iiiiv(DCiNO[j],DCi),DCi);
		}
	  }
	}
	else 
	{
	  // full 
	  if (i<=k1)
	  {
	    DCi=0..Boxes[i]-1;
	  }
	  else 
	  {
	    if (i<=theta)
		{
		  DCi=0;
		  for (j=2;j<=Boxes[i];j=j+1)
		  {
		    DCi=DCi,(j-1)*theta;
		  }
		}
		else 
		{
		  DCi=0;
		}
	  }
	}
	DC[i]=DCi;
  }
  
  // ideal P(t) ... 
  ring Rt=(p^r,a),(t),lp;
  ideal P;
  for (i=1;i<=theta+1;i=i+1)
  {
    P[i]=1;
	if (DC[i]<>-1)
	{
	  s=size(DC[i]);
	  for (j=1;j<=s;j=j+1)
	  {
        P[i]=P[i]*(t-number(a^(DC[i][j])));
	  }
	}
  }
  matrix A=coeffs(P,t); 
  int nr=nrows(A);
  if (nr<k2)
  {
    int nR=k2-nr;
	matrix AA[nR][theta+1];
	A=transpose(concat(transpose(A),transpose(AA)));
  }
  export(P);
  export(A);
  
  // ideal f ... 
  setring Rxy;
  ideal f;
  matrix A=imap(Rt,A);
  number bb,xx,yy;
  // OJO CON LOS INDICES ... 
  for (i=1;i<=theta;i=i+1)
  {
    f[i]=0;
	if (DC[i]==-1)
	{
	  s=0;
	}
	else 
	{
	  s=size(DC[i]);
	}
	s=Boxes[i];
    for (j=0;j<s;j=j+1)
	{
	  // aj == A[j+1,i];
	  // Pij == ( O[i][j+1,1] , O[i][j+1,2] ) 
	  xx=number(O[i][j+1,1]);
	  yy=number(O[i][j+1,2]);
	  bb=number(subst(B[i,j+1],x,xx,y,yy));
	  f[i]=f[i]+(A[j+1,i]*B[i,j+1])/bb;
	}
    f[i]=f[i]*F[i];
  }
  f[theta+1]=F[theta+1];
  export(f);
  
  // Evaluate[i,point]==Evaluate[i][k,j] ... 
  list Evaluate;
  // matrix evaluate[theta+1][k2];
  for (i=1;i<=theta+1;i=i+1)
  {
    matrix evaluate[theta+1][k2];
    for (k=i;k<=theta+1;k=k+1) // k=1 --> k=i 
    {
      for (j=1;j<=Boxes[k];j=j+1)
	  {
	    evaluate[k,j]=number(subst(f[i],x,O[k][j,1],y,O[k][j,2]));
	  }
    }
    Evaluate[i]=evaluate;
    kill evaluate;
  }
  export(Evaluate);
  setring Rt;
  list Evaluate=imap(Rxy,Evaluate);
  export(Evaluate);
  
  if (defring>0)
  {
    setring base_r;
  }
  
  intvec codepars=p,r,n,m,g;
  
  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_NTcode (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;
  
  intvec Boxes=L[3];
  int theta=size(Boxes)-1;
  list DC=L[5];
  def Rt=L[6];
  // intvec L[7]=p,r,n,m,g;
  int p=L[7][1];
  int r=L[7][1];
  // int n=L[7][1];
  // int g=L[7][1];
  // int m=L[7][1];
  
  int i,j,k,s;
  
  setring Rt;
  matrix GB[theta+1][theta+1];
  // algorithm 4.2 
  for (i=1;i<=theta+1;i=i+1)
  {
    if (DC[i]==-1)
	{
	  s=0;
	}
	else 
	{
	  s=size(DC[i]);
	}
	if (s<Boxes[i])
	{
	  for (k=1;k<i;k=k+1)
	  {
	    GB[i,k]=0;
	  }
	  for (k=i;k<=theta+1;k=k+1)
	  {
	    GB[i,k]=0;
	    for (j=1;j<=Boxes[i];j=j+1)
		{
		  GB[i,k]=GB[i,k]+number(Evaluate[i][k,j])*t^(j-1);
		}
	  }
	}
	else 
	{
	  GB[i,i]=t^(Boxes[i])-1;
	}
  }
  export(GB);
  
  if (defring>0)
  {
    setring base_r;
  }
  
  // update CODE 
  CODE[1]=L[1];
  CODE[2]=L[2];
  CODE[3]=Boxes;
  CODE[4]=L[4];
  CODE[5]=DC;
  CODE[6]=Rt;
  CODE[7]=L[7];

  return(CODE);
}

;