X := Poincare(CurrentRing());

Define IdealPP(Points)
  Return IdealOfProjectivePoints(Points)
End;

Define OSofPP(Points)
Help '
 A projective point P is a list P := [a_0,..,a_n]  where the a_i are not all
zero.

 Points is a list of projective points Points := [P_1,..,P_m] .

 An O-Sequence is represented as a list [ FirstValuesList , Polynomial] ,
e.g. the O-Sequence (1,2,3,4,5,6,6,6....) is represented as [ [2,3,4,5],6].

 OSofPP is the O-Sequence whose values are the Hilbert function of the
coordinate ring of the set of points Points.
';

 I := IdealPP(Points);
 H := @Hilbert(CurrentRing()/I);
 If H[1]=[] Then Return [[], H[2]] End;
 Return [Tail(H[1]), H[2]];
End;

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

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

Define IsDiffOS(S)
  Return IsOS(DiffOS(S))
End;

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
  End;
  L := [ I In 1..LenS1 | S[1,I]<>Subst(S[2], Indet(1), I) ];
  If L=[] Then M := 0 Else M := Max(L) End;
  Return [ First(S[1], M), S[2] ];
End;

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

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

 
Define ProS1Anna(S)
  S := Reduce(S);
  N := Value(S,1);
  If N=2 Then Return [ [S[1], S[2]-1], -1]  End;
  L := [];
  D := 2;
  While Value(S,D) - Value([L,0], D-2) >= Bin(N+D-2, D)  Do
    Append(L, Value(S,D)-Bin(N+D-2, D) );
    D := D+1;
  End;
  Return [ [L, Value([L,0], D-2)], D-1];
End; -- ProS1Anna

Define GMRH(S)
  L := ProS1Anna(S);
  Return L[2];
End;

Define GMRS1(S)
  L := ProS1Anna(S);
  Return Reduce(L[1]);
End; -- GMRS1

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

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

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

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

Define Indent(N)
  If N <= 0 Then Return "" End; 
  Return Sum(["   " | I In 1..N ]) 
End; 

Define PTVofHF(S, Idt)
     N := Value(S,1)-1;
    If N=0 Then Return 1 End;
  If N=1 Then Return([S[2]]) End;          
  If N=2 Then 
         DS := DiffOS(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]);
    End;
    Return L;
  End; 
  TV1 := PTVofHF(GMRS1(S), Idt+1);
     TV2 := PTVofHF(GMRS2(S), Idt+1);
     If ListDepth(TV1)=N Then Return Concat(TV1, [TV2]) End;
  Return [MakeNList(N-1,TV1), TV2];
End; -- PTVofHF

Define TVofHF(S)
  N := NumIndets();
  If Value(S,1)>N Then Error("Invalid Hilbert function for this ring") End;
  If Not(IsDiffOS(S)) Then Error("Not a differentiable OS") End;
  Return MakeNList(N-1, PTVofHF(S,0));
End;

Define HFofTV(TV)
  L := [1];
  LD := ListDepth(TV);
  I := 1;
    While L[I]<>HFNT(LD, TV, I) Do
         Append(L, HFNT(LD, TV, I) );
    I := I+1;
  End; 
  Return Reduce([Tail(L), L[I]]);
End;

Define PrettyTV(T);
  For I := 1 To Len(T) Do PrintLn NewLine, T[I]  End;
End;

Define BETTIofTV(TV)
  If Len(TV)=1 Then
    If ListDepth(TV)=2 Then
      Return([[Comp(TV[1],1)+1],[Comp(TV[1],1),1]]);
    End;
    Rold := BETTIofTV(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;
    End;
    Append(Rnew,Sorted(Concat(Rold[Len(Rold)],[1])));
    Return Rnew;
  End;
  ResY := BETTIofTV(Diff(TV,[TV[Len(TV)]]));
  ResX := BETTIofTV([TV[Len(TV)]]);
  ResYuX := [];
  For J := 1 To Len(ResX)-1 Do
    Append(ResYuX,Sorted(Concat([I+1| I In ResY[J]],ResX[J])));
  End;
  J := Len(ResX)-1;
  Append(ResYuX,Sorted(Concat([I+1| I In ResY[J+1]],Diff(ResX[J+1],[1]))));
  Return ResYuX;
End;

Define ShiftedR(L)
  While Not L=[] Do
    D := L[1];
    J := 0;
    While Not L=[] And L[1]=D Do 
      L := Tail(L);
      J := J+1;
    End;
    If J=1  Then Print         ' R(-', D ,')'  End;
    If J<>1 Then Print ' R^', J, '(-', D, ')'  End;
    If Not L=[] Then Print ' (+)' End;
  End;
End;

Define RESofTV(TV);
  PrintLn;
  Print('0 ---->');
  ResL := BETTIofTV(TV);
  For I := 1 To Len(ResL) Do
    ShiftedR(ResL[I]);
    Print(' ---->');
  End;
  Print(' I ----> 0');
End;


Define RESofBETTI(B)
  If B=[] Then Return[]; End;
  Print('0 ----> ');
  LenB := Len(B);
  For I := 1 To Len(B) Do
    M := B[I];
    While M<>[] Do
      Exponent := 1;
      For J := 2 To Len(M) Do
        If M[J]=M[1] Then 
          Exponent := Exponent + 1; 
        End;
      End;
      If Exponent=1 Then 
        Print('R(',-M[1],')'); 
      End;
      If Exponent<>1 Then 
        Print('R^',Exponent,'(',-M[1],')'); 
      End;
      M := Last(M,(Len(M)-Exponent));
      If M<>[] Then
        Print(' (+) ');
      End;
    End;
    If I<>Len(B) Then
      Print(' ----> ');
    End;
  End;
  Print(' ----> I ----> 0');
End;

Define PrependN(Var L, N) -- add the number N to the beginning
-- of each element in a list of lists
  For I  :=  1 To Len(L) Do
    Insert(L[I],1,N);
  End;
End;

Define PointsFromType(T)
  N  :=  ListDepth(T);
  If N = 0 Then Return [[1]] End;
  If N = 1 Then
    Return [[X-1,1] | X In 1..T[1]];
  End;
  T  :=  Reversed(T); 
  Result  :=  [];
  For N  :=  0 To Len(T)-1 Do
    X  :=  PointsFromType(T[N+1]);
    PrependN(X,N); -- this function defined above
    Result  :=  Concat(Result,X);
  End;
  Return Result;
End;

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

Define DiffHFofTV(TV)
  M := HFofTV(TV);
  K := DiffOS(M);
  Return K;
End;

Define TVofDiffHF(S)
  A := Concat([1],S[1]);
  B := Integrate(A);
  N := B[Len(B)];
  B := Tail(B);
  B := First(B,(Len(B)-1));
  C := Concat([B],[N]);
  D := TVofHF(C);
  Return D;
End;