Package $contrib/typevectors

Alias TV := $contrib/typevectors,
      HP := $hp;

/*--------------------------------------------*\

   Run CoCoA and type:

     Alias TV := $contrib/typevectors;
     TV.About();
     TV.Man();

\*--------------------------------------------*/

Define About()
  PrintLn "
    KeyWords : type vectors, Hilbert functions, points.
    Author   : E. Carlini, M. Stewart
    Version  : CoCoA-4.2
    Date     : 8 May 2001 (minor revision 3 July 2002 and 31 July 2002)
"
EndDefine; -- About

------[   Manual   ]--------

Define Man()
  PrintLn "
Suggested alias for this package:
    Alias TV := $contrib/typevectors;

SYNTAX
    TV.FromHF( HF: TAGGED('$hp.Hilbert')): LIST
    TV.ToHF( TV: LIST ): TAGGED('$hp.Hilbert')
    TV.PrintRes( TV: LIST )
    TV.Points( TV: LIST): LIST

    TV.Alpha( TV: LIST ): INTEGER
    TV.Sigma( TV: LIST ): INTEGER
    TV.Type( TV: LIST ):  INTEGER
    TV.IsValid( TV: LIST ): BOOL

DESCRIPTION

This package enables one to compute with ""type vectors"" as described in
the article by A.Geramita, T.Harima, Y.Shin entitled ""An alternative to
the Hilbert function for the ideal of a finite set of points in P^n"",
published in Illinois J.Math. vol.45 (2001), no. 1, pages 1--23.

Let T be an n-type vector, we represent a 2-type vector as [[d_1],...,[d_r]]
so that the representation of T is obtained by recursion.

   The function
           TV.FromHF(H)
   returns a LIST representing the type vector corresponding to the Hilbert
   funtion H.

   The function
           TV.ToHF(TV)
   returns a TAGGED('$hp.Hilbert') object that is the Hilbert function
   associated to TV.

   The function
           TV.Points(TV)
   returns a LIST of projective points having Hilbert function
   associated to TV.

   The function
           TV.PrintRes(TV)
   prints the resolution of the ideal of a K-configuration described by TV.

   The functions
           TV.Alpha(TV)
           TV.Sigma(TV)
   return an INTEGER representing the alpha and sigma quantities of TV.

   The function
           TV.Type(TV)
   returns the INTEGER n such that TV is a n-type vector.

   The function
           TV.IsValid(TV)
   returns TRUE if TV is a valid type vector, FALSE otherwise.

>EXAMPLE<

   Use P ::= Q[x,y,z,t];
   HF := Hilbert(R/Ideal(x, y^3, z^2-tx));  HF;
   T := TV.FromHF(HF);  T;
   T := TV.FromHF([[1, 3, 5], 6]);  T;
   TV.ToHF([[[2], [4]]]);
   Pts := TV.Points(T);  Pts;
   Hilbert(P/IdealOfProjectivePoints(Pts));
   TV.PrintRes(T);

   TV.Alpha(T);
   TV.Sigma(T);
   TV.IsValid(T);
   TV.Type(T);
";
EndDefine; -- Man

----------------------------------------------------------------------

----------------------------[ Type Vectors ]--------------------------

Define FromHF(HF)
  S := TV.HFtoMyHF(HF);
  N := NumIndets();
  If Not(TV.MyIsDiffOS(S)) Then
    Error("TV.FromHF: Not a differentiable O-sequence: ", HF)
  EndIf;
  Return TV.MakeNList(N-1, TV.PTVofHF(S,0));
EndDefine; -- FromHF

Define ToHF(TV)
  L := [1];
  LD := TV.ListDepth(TV);
  I := 1;
  While L[I]<>TV.HFNT(LD, TV, I) Do
    Append(L, TV.HFNT(LD, TV, I) );
    I := I+1;
  EndWhile;
  Return TV.MyHFtoHF(TV.Reduce([Tail(L), L[I]]));
EndDefine; -- ToHF

Define PrintRes(TV);
  ResL := TV.Betti(TV);
  Print '0 --> ';
  For I := 1 To Len(ResL) Do
    TV.PrintShiftedR(ResL[I]);
    Print ' --> ';
  EndFor;
  Print RingEnv();
EndDefine; -- PrintRes

Define Points(TV)
  N := TV.ListDepth(TV);
  If N = 0 Then Return [[1]]; EndIf;
  If N = 1 Then Return [[X-1,1] | X In 1..TV[1]]; EndIf;
  TV := Reversed(TV);
  Result := [];
  For N := 0 To Len(TV)-1 Do
    X := TV.Points(TV[N+1]);
    TV.PrependN(X,N);
    Result := Concat(Result,X);
  EndFor;
  Return Result;
EndDefine; -- Points

Define Alpha(TV)
  N:=TV.ListDepth(TV);
  If N = 0 Then Return -1;      EndIf;
  If N = 1 Then Return TV[1];   EndIf;
  If N >=2 Then Return Len(TV); EndIf;
EndDefine; -- Alpha

Define Sigma(TV)
  N:=TV.ListDepth(TV);
  If N = 0 Then Return 1;     EndIf;
  If N = 1 Then Return TV[1]; EndIf;
  If N >=2 Then Return TV.Sigma(TV[Len(TV)]); EndIf;
EndDefine; -- Sigma

Define IsValid(TV)
  N:=TV.ListDepth(TV);
  If N = 0 Then Return TV = 1;  EndIf;
  If N = 1 Then Return TV[1]>0; EndIf;
  If N >=2 Then
    Foreach TT In TV Do
      If Not TV.IsValid(TT) Then Return FALSE; EndIf;
    EndForeach;
    For I:=1 To Len(TV)-1 Do
      If TV.Alpha(TV[I+1])<=TV.Sigma(TV[I]) Then Return FALSE; EndIf;
    EndFor;
  EndIf;
  Return TRUE;
EndDefine; -- IsValid

Define Type(TV)
  If Not TV.IsValid(TV) Then
    Error('TV.Type: not a valid type vector: ', TV);
  EndIf;
  Return TV.ListDepth(TV);
EndDefine; -- Type

Define IsValid(TV)
  N:=TV.ListDepth(TV);
  If N = 0 Then Return TV = 1;  EndIf;
  If N = 1 Then Return TV[1]>0; EndIf;
  If N >=2 Then
    Foreach TT In TV Do
      If Not TV.IsValid(TT) Then Return FALSE; EndIf;
    EndForeach;
    For I:=1 To Len(TV)-1 Do
      If TV.Alpha(TV[I+1])<=TV.Sigma(TV[I]) Then Return FALSE; EndIf;
    EndFor;
  EndIf;
  Return TRUE;
EndDefine; -- IsValid

----------------------------[ AUXILIARY ]----------------------

----------------------------[ OSequence ]----------------------

Define IsOS(OS)
  S:=TV.OStoMyOS(OS);
  S := TV.Reduce(S);
  S1 := S[1];
  If S1=[] Then Return TRUE EndIf;
  Foreach X In S1 Do If X<0 Then Return FALSE EndIf EndForeach;
  If S[2]<>0 And Deg(S[2])>0 Then
    Error("TV.IsOS only works with 0 dimensional O-Sequence")
  EndIf;
  L := Len(S1);
  For I := 1 To L-1 Do
    If EvalBinExp(BinExp(S1[I],I),1,1) < S1[I+1] Then Return FALSE EndIf;
  EndFor;
  If EvalBinExp(BinExp(S1[L],L),1,1) < LC(S[2]) Then Return FALSE EndIf;
  Return TRUE;
EndDefine; -- IsOS

Define DiffOS(OS)
Help'
 Given an O-Sequence  S  this function computes its first difference.

 Syntax= DiffOS( [ ListOfNonNegativeIntegers , APolynomialOfR[t] )
';
  Using Qt Do
    S:=TV.OStoMyOS(OS);
    S := TV.Reduce(S);
    Return TV.MyOStoOS([[TV.Value(S,I+1)-TV.Value(S,I) | I In 0..Len(S[1])],
		      Subst(S[2], t, t+1) - S[2] ]);
  EndUsing;
EndDefine; -- DiffOS

Define IsDiffOS(S)
  Return TV.IsOS(TV.DiffOS(S))
EndDefine; -- IsDiffOS

Define Reduce(S)
  Help'
Given an O-Sequence S := [ [b_1,..,b_r] , P(t) ],
Reduce(S) := [ [b_1,..b_s] , P(t) ] is such that P(s) /neq b_s ( s<=r )
';
  LenS1 := Len(S[1]);
  If LenS1=0 Or S[1,LenS1]<>Subst(S[2], Indet(1), LenS1) Then
    Return S
  EndIf;
  L := [ I In 1..LenS1 | S[1,I]<>Subst(S[2], Indet(1), I) ];
  If L=[] Then M := 0 Else M := Max(L) EndIf;
  Return [ First(S[1], M), S[2] ];
EndDefine; -- Reduce

Define Integrate(L)
  If Len(L)=1 Then Return L EndIf;
  If Len(L)=0 Then Return L EndIf;
  M := [];
  Append(M,L[1]);
  For I := 2 To Len(L) Do
    Append(M,(M[I-1]+L[I]));
  EndFor;
  Return M;
EndDefine; -- Integrate

----------------------------[ Various ]------------------

Define Value(S,I)
  If I=0 Then Return 1 EndIf;
  If I<0 Then Return 0 EndIf;
  If I<=Len(S[1]) Then Return S[1][I] EndIf;
  Return Qt :: Eval(S[2], [I]);
EndDefine; -- Value

Define ListDepth(L)
  If Type(L)<>LIST Then Return 0 EndIf;
  If Flatten(L)=L  Then Return 1 EndIf;
  Return Max([ TV.ListDepth(LL) | LL In L ])+1;
EndDefine; -- ListDepth

Define MakeNList(N, L)
  For I:=1 To N-TV.ListDepth(L) Do L := [L] EndFor;
  Return L;
EndDefine; -- MakeNList

Define PTVofHF(S, Idt)
  N := TV.Value(S,1)-1;
  If N=0 Then Return 1 EndIf;
  If N=1 Then Return([S[2]]) EndIf;
  If N=2 Then
    DS := TV.MyDiffOS(S);
    DS := Concat([1], DS[1]);
    L := [];
    While Len(DS)<>0 Do
      L := Concat([[Len(DS)]], L);
      DS := Diff([(DS[I]-1) | I In 1..Len(DS) ], [0]);
    EndWhile;
    Return L;
  EndIf;
  TV1 := TV.PTVofHF(TV.GMRS1(S), Idt+1);
  TV2 := TV.PTVofHF(TV.GMRS2(S), Idt+1);
  If TV.ListDepth(TV1)=N Then Return Concat(TV1, [TV2]) EndIf;
  Return [TV.MakeNList(N-1,TV1), TV2];
EndDefine; -- PTVofHF

Define GMRS1(S)
  L := TV.ProS1(S);
  Return TV.Reduce(L[1]);
EndDefine; -- GMRS1

Define ProS1(S)
  S := TV.Reduce(S);
  N := TV.Value(S,1);
  If N=2 Then Return [ [S[1], S[2]-1], -1] EndIf;
  L := [];
  D := 2;
  While TV.Value(S,D) - TV.Value([L,0], D-2) >= Bin(N+D-2, D)  Do
    Append(L, TV.Value(S,D)-Bin(N+D-2, D) );
    D := D+1;
  EndWhile;
  Return [ [L, TV.Value([L,0], D-2)], D-1];
EndDefine; -- ProS1

Define GMRS2(S)
  S := TV.Reduce(S);
  N := TV.Value(S,1);
  S1 := TV.GMRS1(S);
  H := TV.GMRH(S);
  If H=-1 Then Return([[], Bin(N-2+Indet(1), N-2)])  EndIf;
  C := [ [Bin(N-2+I, N-2) | I In 1..H ], 0];
  For J := H+1 To Len(S[1]) Do
    Append(C[1], TV.Value(S,J)-TV.Value(S1,H-1));
  EndFor;
  C[2] := S[2]-TV.Value(S1,H-1);
  Return TV.Reduce(C);
EndDefine; -- GMRS2

Define GMRH(S)
  L := TV.ProS1(S);
  Return L[2];
EndDefine; -- GMRH

Define HFNT(LD, TV, J)
  If J<0 Then Return 0 EndIf;
  If J=0 Then Return 1 EndIf;
  LenTV := Len(TV);
     If LD=1 Then
       Return Min(TV[LenTV], J+1);
     EndIf;
     Return Sum([ TV.HFNT(LD-1, TV[I], J-LenTV+I) | I In 1..LenTV ]);
EndDefine; -- HFNT

Define PrintShiftedR(L)
  While Not L=[] Do
    D := L[1];
    J := 0;
    While Not L=[] And L[1]=D Do
      L := Tail(L);
      J := J+1;
    EndWhile;
    If J=1  Then Print RingEnv(), '(-', D ,')'  EndIf;
    If J<>1 Then Print RingEnv(), '^', J, '(-', D, ')'  EndIf;
    If Not L=[] Then Print '(+)' EndIf;
  EndWhile;
EndDefine; --  PrintShiftedR

Define Betti(TV)
  If Len(TV)=1 Then
    If TV.ListDepth(TV)=2 Then
      Return([[Comp(TV[1],1)+1],[Comp(TV[1],1),1]]);
    EndIf;
    Rold := TV.Betti(TV[1]);
    Rnew := [];
    Append(Rnew,[I+1 | I In Rold[1]]);
    J := 1;
    While J<Len(Rold) Do
      Append(Rnew,Sorted(Concat(Rold[J],[I+1|I In Rold[J+1]])));
      J := J+1;
    EndWhile;
    Append(Rnew,Sorted(Concat(Rold[Len(Rold)],[1])));
    Return Rnew;
  EndIf;
  ResY := TV.Betti(Diff(TV,[TV[Len(TV)]]));
  ResX := TV.Betti([TV[Len(TV)]]);
  ResYuX := [];
  For J := 1 To Len(ResX)-1 Do
    Append(ResYuX,Sorted(Concat([I+1| I In ResY[J]],ResX[J])));
  EndFor;
  J := Len(ResX)-1;
  Append(ResYuX,Sorted(Concat([I+1| I In ResY[J+1]],Diff(ResX[J+1],[1]))));
  Return ResYuX;
EndDefine; -- Betti

Define PrependN(Var L, N)
  For I  :=  1 To Len(L) Do
    Insert(L[I],1,N);
  EndFor;
EndDefine; -- PrependN

----------------------------[ MyOSequence ]----------------------

Define MyIsOS(S)
  S := TV.Reduce(S);
  S1 := S[1];
  If S1=[] Then Return TRUE EndIf;
  Foreach X In S1 Do If X<0 Then Return FALSE EndIf EndForeach;
  If S[2]<>0 And Deg(S[2])>0 Then
    Error("This function only works with 0 dimensional O-Sequence")
  EndIf;
  L := Len(S1);
  For I := 1 To L-1 Do
    If EvalBinExp(BinExp(S1[I],I),1,1) < S1[I+1] Then Return FALSE EndIf;
  EndFor;
  If EvalBinExp(BinExp(S1[L],L),1,1) < LC(S[2]) Then Return FALSE EndIf;
  Return TRUE;
EndDefine; -- MyIsOS

Define MyDiffOS(S)
Help'
 Given an O-Sequence  S  this function computes its first difference.

 Syntax= MyDiffOS( [ ListOfNonNegativeIntegers , APolynomialOfR[t] )
';
  T := Indet(1);
  S := TV.Reduce(S);
  Return [ [TV.Value(S,I+1)-TV.Value(S,I) | I In 0..Len(S[1])],
	  Subst(S[2], T, T+1) - S[2] ];
EndDefine; -- MyDiffOS

Define MyIsDiffOS(S)
  Return TV.MyIsOS(TV.MyDiffOS(S))
EndDefine; -- MyIsDiffOS

Define MyIntegrate(L)
  If Len(L)=1 Then Return L EndIf;
  If Len(L)=0 Then Return L EndIf;
  M := [];
  Append(M,L[1]);
  For I := 2 To Len(L) Do
    Append(M,(M[I-1]+L[I]));
  EndFor;
  Return M;
EndDefine; -- MyIntegrate

----------------------------[ Aggiustamenti ]------------------

Define MyOStoOS(S)
  If Len(S[1])>0 And S[1][1]=1 Then  Return S; EndIf;
  Return [Concat([1],S[1]),S[2]];
EndDefine; -- MyOStoOS

Define OStoMyOS(S)
  S1:=S[1];
  If Len(S1)<1 Then Return S; EndIf;
  If S1[1]<>1  Then Return S; EndIf;
  Return [Tail(S[1]),S[2]];
EndDefine;

Define HFtoMyHF(H)
  Return TV.OStoMyOS(@H);
EndDefine; -- HFtoMyHF

Define MyHFtoHF(H)
  Return Tagged(TV.MyOStoOS(H),"$hp.Hilbert");
EndDefine; -- MyHFtoHF

EndPackage;