-- Cioffi F. & Sessa C.
-- Date: 30/07/2022

-----------------------------
-- MainKaLu(I,J,K,L,I1,J1) --
-----------------------------

-- AIM:
-- Main function for the computation of the
-- Kazhdan-Lusztig polynomials with respect
-- to two Schubert varieties S and S'
-- with S' contained in S. Q will be the
-- index corresponding to the S-variety S'_F
--
-- INPUT:
-- I,J,K,L - data defining S
-- I1,J1,K,L - data defining S'
--
-- LOCAL VARIABLES:
-- ADM - 4-indices list which contains all
--       0-admissible vectors T such that
--       0 <= T <= Q. ADM[H,Z] consists of
--       two lists; the latter is the vector
--       T (with respect to the essential flag
--       of 0) and the former is the vector of
--       pointers TA to its essential conditions.
--       The value H-1 represents the distance between
--       T and 0, whereas Z stands for the position
--       of T in the list ADM[H], which is ordered
--       with respect to the number of essential
--       conditions of the vectors.
--
-- SUBFUNCTIONS:
-- SchubertVariety, PQAmmissibile, Divide,
-- Actual, OrdineCondizioni, TauPQ,
-- TauPQAmmissibile, DistanzaSigma

Define MainKaLu(I,J,K,L,I1,J1)
  TopLevel SchubertVariety;
  TopLevel PQAmmissibile; 
  TopLevel Divide; 
  TopLevel Actual; 
  TopLevel TauPQ; 
  TopLevel TauPQAmmissibile; 
  TopLevel DistanzaSigma;
  TopLevel OrdineCondizioni;

  B0Q:=[];
  Omega:=Len(I);
  Omega1 := Len(I1);
  Zero:=[0 | Variabile In 1..Omega];
  Zero1 := [0 | Variabile In 1..Omega1];

-- Check if I,J,K,L give a Schubert variety

  If Not(SchubertVariety(Zero,I,J,K,L)) Then
    PrintLn "The data of S do not satisfy conditions of Remark 2.7";
  Else
    ZeroA:=Actual(Zero,I,J,K,L);
    Omega := Len(ZeroA);

-- Check if S is smooth.

    If Omega = 0 Then
      PrintLn "All conditions of S are redundant."
    Else

-- Take only the essential conditions
-- of S. Reduce the length of Zero and 
-- change I and J concordingly;
-- that is, consider only the elements
-- indexed by ZeroA.

      I0:=[I[ZeroA[H]] | H In 1..Omega];
      J0:=[J[ZeroA[H]] | H In 1..Omega];
      Zero :=[0 | H In 1..Omega];
      ZeroA := [H | H In 1..Omega];

-- Check if I1,J1,K,L give a Schubert variety

      If Not(SchubertVariety(Zero1,I1,J1,K,L)) Then
        PrintLn "The data of S' do not satisfy conditions of Remark 2.7";
      Else

-- Compute the sequences Lambda0 and
-- Lambda1 associated to S and S'

        Lambda0 := Lambdax(I0, J0, K, L);
        Lambda1 := Lambdax(I1, J1, K, L);

-- Check if S' is contained in S.

        T := 1;
        CondLambda := True;
        While T <= K And CondLambda Do
          If Lambda0[T] > Lambda1[T] Then
            CondLambda := False;
          Else
            T := T + 1;
          EndIf;
        EndWhile;

        If Not(CondLambda) Then
          PrintLn "S' is not contained in S.";
        Else

-- Replace S' with S'_F, which is
-- an S-variety of the form Delta_Q

          Q := Replacement(I0, J0, K, L, Lambda1);
          QA := Actual(Q,I0,J0,K,L);

-- Check if Q = 0

          If Q = Zero Then
            B0Q := 1;
            PrintLn "Q = 0, thus A_0Q = G_0Q = B_0Q = 1.";
          Else
            ADM:=[[[Zero,ZeroA]]];

-- Compute the list of all S varieties
-- Delta_Tau strictly contained in S
-- and strictly containing Delta QP.

            V:=TauPQ(Zero,Q); 
            ADM:=concat(ADM,TauPQAmmissibile(Zero,I0,J0,K,L,V));

-- Change the order of the elements of ADM[H]
-- (for all H) so that they are ordered with
-- respect to the (increasing) number of
-- essential condition.

            For H:=2 To len(ADM) Do
              SortBy(ref ADM[H],OrdineCondizioni);
            EndFor;

-- Add Q as the last element of ADM

            ADM:=concat(ADM,[[[Q,QA]]]);

-- Print ADM one row at a time

            PrintLn "ADM := ";
            For H := 1 To Len(ADM) Do
              PrintLn ADM[H];
            EndFor;

-- Compute B0Q

            B0Q:=DistanzaSigma(ADM,I0,J0,K,L);
          EndIf;
        EndIf;
      EndIf;
    EndIf;
  EndIf;
  Return B0Q;
EndDefine;

----------------------
-- Lambdax(I,J,K,L) --
----------------------

-- AIM:
-- Computation of the sequence Lambda associated
-- to the Schubert variety S given by I, J, K, L.
--
-- INPUT:
-- I, J, K, L - Data defining S.
--

Define Lambdax(I,J,K,L)

  Omega := Len(I);
  Lambda := [L - K - J[1] + I[1] | H In 1..I[1]];
  For T := 2 To Omega Do
    Temp := [L - K - J[T] + I[T] | H In (I[T - 1] + 1)..I[T]];
    Lambda := Concat(Lambda, Temp);
  EndFor;

  If I[Omega] < K Then
    Temp := [0 | H In (I[Omega] + 1)..K];
    Lambda := Concat(Lambda, Temp);
  EndIf;
  Return Lambda;
EndDefine;

----------------------------
-- Replacement(J, Lambda) --
----------------------------

-- AIM:
-- Replace S' with the S-variety S'_F.
--
-- INPUT:
-- I, J - Vectors consisting of the 
--        conditions and the dimensions
--        of the vector spaces of the flag
--        of S, respectively.
-- Lambda - sequence associated to S'.
--
-- OUTPUT:
-- Q - vector such that S'_F = Delta_Q.
--
-- LOCAL VARIABLES:
-- Vect - vector needed for the computation
--        of the conditions of S' with
--        respect to the flag of S.
-- I1 - vector of the pair with flag F
--      defining Delta_Q
--


Define Replacement(I, J, K, L, Lambda)

  Omega := Len(J);
  I1 := [];
  Vect := [L - K + H - Lambda[H] | H In 1..K];
  Append(ref Vect, L + 1);

-- For each T, find U such that
-- J[T] <= Vect[U] - 1

  T := 1;
  U := 1;
  While T <= Omega And U <= K + 1 Do
    If J[T] <= Vect[U] - 1 Then
      Append(ref I1, U - 1);
      T := T + 1;
    Else
      U := U + 1;
    EndIf;
  EndWhile;

  Q := I1 - I;
  Return Q;
EndDefine;

-----------------------------
-- CalcoloMTX(T,I,J,K,L,X) --
-----------------------------

-- AIM:
-- computation of dim Delta_T - dim Delta_X by
-- means of the formula obtained by interpreting
-- such property by means of the Ferrer's
-- diagrams.
--
-- INPUT:
-- I, J, K, L - data defining the Schubert
--              variety S
-- T, X - vectors representing Delta_T and Delta_X
--
-- OUTPUT:
-- MTX - it is dim Delta_T - dim Delta_X

Define CalcoloMTX(T,I,J,K,L,X)

  TopLevel Lambdax;

  LambdaX := Lambdax(I + X,J,K,L);
  LambdaT := Lambdax(I + T,J,K,L);
  MTX:=sum(LambdaX-LambdaT);
  If MTX<0 Then
    println "WARNING: negative MTX";
  EndIf;
  Return MTX;
EndDefine;

------------------------
-- Binomiale(Omega,D) --
------------------------

-- AIM:
-- For any N <= Omega and K < D, the number
-- Bin(N-1+K,K) of terms in N variables of
-- degree K is computed.
--
-- INPUT:
-- Omega, D - integers 
--
-- OUTPUT:
-- B - list of all binomials Bin(N-1+K,K).

Define Binomiale(Omega,D)
  B:=[[ 0 | K In 1..(D+1)] | N In 1..Omega];
  For N:=1 To Omega Do
    B[N,1]:=1;
    For K:=1 To D Do
      B[N,K+1]:=(B[N,K]*(K+N-1))/K ;
    EndFor;
  EndFor;
  Return B;
EndDefine;

-----------------------------
-- OrdineCondizioni(T1,T2) --
-----------------------------
-- AIM:
-- Given two lists correspoding to
-- two S-varieties Delta_P and Delta_Q and
-- their essential flags (i.e. each list
-- contains the vectors of increments P and Q
-- and the lists PA and QA of pointers to
-- their essential conditions), we check if the
-- essential flag of Delta_Q is shorter than the
-- one of Delta_P.
--
-- INPUT:
-- T1, T2: lists as described above.
--
-- OUTPUT:
-- B: True if len(T1[2])<=len(T2[2])

Define OrdineCondizioni(T1,T2)
  B:=False;
  If len(T1[2])<=len(T2[2]) Then
    B:=True;
  EndIf;
  Return B;
EndDefine;

-------------------
-- Divide(L1,L2) --
-------------------

-- AIM:
-- Given two vectors L1 and L2 of the same
-- length, we check if L1[T] <= L2[T] for
-- all T.
--
-- INPUT:
-- L1, L2 - lists of integers
--
-- OUTPUT:
-- A - logical variable which is TRUE if
--     L1[T] <= L2[T] for all T.

Define Divide(L1,L2)
  A:=True;
  K:=1;
  While K<=len(L1) And A Do
    If L1[K]>L2[K] Then
      A:=False;
    Else 
      K:=K+1;
    EndIf;
  EndWhile;
  Return A;
EndDefine;

----------------
-- TauPQ(P,Q) --
----------------
-- AIM:
-- Given two vectors P and Q, we find all
-- vectors V with P < V < Q.
--
-- INPUT:
-- P, Q - vectors corresponding to
--        S-admissible varieties
-- OUTPUT:
-- T - list of all vectors between P and Q.
--     For all H and K, T[H,K] is the H-th
--     vector of degree K + Sum(P), with
--     respect to the degLex order.
--
-- SUBFUNCTIONS:
-- Binomiale, Divide

Define TauPQ(P,Q)
  TopLevel Binomiale;
  TopLevel Divide;

-- Compute all vectors whose distance from S 
-- is less than the one of Q. They will be 
-- saved w.r.t. degLex order.

  Omega:=len(P);
  D:=sum(Q);

  B:=Binomiale(Omega,D-1);

  E:=[[0 | K In 1..Omega] | H In 1..Omega];
  For K:=1 To Omega Do
    E[K,K]:=1;
  EndFor;
  Tau:=[ 0 | K In 1..(D-1)];

  If D>1 Then
    Tau[1]:=E;

    For K:=2 To D-1 Do
      TK:=[Tau[K-1,H]+E[1] | H In 1..B[Omega,K]];
      For J:=2 To Omega Do
        I:=sum([B[Omega-H+1,K-1] | H In 1..(J-1)])+1;
        TK:=concat( TK,[Tau[K-1,C]+E[J] | C In I..(B[Omega,K])]);
      EndFor;
      Tau[K]:=TK;
    EndFor;
  EndIf;

-- Take only the terms whose distance from S is greater
-- than the one of P

  T:=[[] | K In (sum(P)+1)..(D-1)];
  For K:=sum(P)+1 To D-1 Do
    LTK:=Len(Tau[K]);
    For H:=1 To LTK Do
      If Divide(P,Tau[K,H]) And Divide(Tau[K,H],Q) Then
        T[K-sum(P)]:=concat(T[K-sum(P)],[Tau[K,H]]);
      EndIf;
    EndFor;
  EndFor;
  Return T;
EndDefine;

--------------------------------
-- SchubertVariety(P,I,J,K,L) --
--------------------------------

-- AIM:
-- Given a Schubert variety S, we check if
-- a vector P corresponds to an S-variety.
-- In other words, we check conditions (4),
-- written in the form of Remark 2.6.
--
-- INPUT:
-- P - vector of integers
-- I, J, K, L - data defining the Schubert
--              variety S
--
-- LOCAL VARIABLES:
-- M, N - conditions on the components of P,
--        except for P[1].
--
-- OUTPUT:
-- A - logical variable which is TRUE if
--     P corresponds to an S-variety.

Define SchubertVariety(P,I,J,K,L)
  A:=True;
  Omega:=len(I);
  If 0 <= P[1] And P[1] <= J[1]-I[1] Then
    Alpha:=2;
    While Alpha<Omega And A Do
      M:=max(0,I[Alpha-1]+P[Alpha-1]-I[Alpha]);
      N:=min(J[Alpha]-J[Alpha-1]+I[Alpha-1]+P[Alpha-1]-I[Alpha], J[Alpha]-I[Alpha]);
      If M<= P[Alpha] And P[Alpha]<=N Then
        Alpha:=Alpha+1;
      Else
        A:=False;
      EndIf;
    EndWhile;
    If A Then
      M:=max(0,I[Alpha-1]+P[Alpha-1]-I[Alpha],K-L+J[Alpha]-I[Alpha]);
      N:=min(J[Alpha]-J[Alpha-1]+I[Alpha-1]+P[Alpha-1]-I[Alpha], J[Alpha]-I[Alpha], K-I[Alpha]);
      If M > P[Alpha] Or P[Alpha] > N Then
        A:=False;
      EndIf;
    EndIf;
  Else
    A := False;
  EndIf;
  Return A;
EndDefine;

-----------------------
-- Actual(P,I,J,K,L) -- 
-----------------------
-- Note: this function should be called
-- "Essential", but we agreed not to
-- update its name.
-----------------------

-- AIM:
-- Given an S-variety Delta_P, we
-- determine its essential conditions.
--
-- INPUT:
-- I, J, K, L - data of the Schubert
--              variety S
-- P - vector corresponding to the
--     S-variety Delta_P
--
-- OUTPUT:
-- PA - vector containing the integers
--      corresponding to the numbers of
--      the essential conditions (e.g.
--      1 in PA means that the 1st
--      condition of P is essential).

Define Actual(P,I,J,K,L)
  Omega:=len(P);
  If Omega>=2 Then
    PA:=[];
    PI:=P+I;
    If PI[2]-PI[1]<J[2]-J[1] Then
      PA:=concat(PA,[1]);
    EndIf;
    For Alpha:=2 To Omega-1 Do
      If PI[Alpha+1]-PI[Alpha]<J[Alpha+1]-J[Alpha] And PI[Alpha-1]<PI[Alpha] Then
        PA:=concat(PA,[Alpha]);
      EndIf;        
    EndFor;
    If K+J[Omega]<L+PI[Omega] And PI[Omega-1]<PI[Omega] Then
      PA:=concat(PA,[Omega]);
    EndIf;     
  Else
   PA:=[1];
  EndIf;
  Return PA;
EndDefine;

--------------------------
-- PQAmmissibile(PA,QA) --
--------------------------

-- AIM:
-- Given P and Q, we check whether
-- Delta_Q is a Delta_P variety.
--
-- INPUT:
-- QA, PA - lists containing the pointers
--          to the essential conditions
--          of Q and P
--
-- OUTPUT:
-- B - logical variable which is TRUE if
--     Q is P-admissible

Define PQAmmissibile(PA,QA)
  B:=True;
  OmegaPA:=len(PA);
  OmegaQA:=len(QA);
  If OmegaQA<=OmegaPA Then
    Z:=1;
    While B And Z<=OmegaQA Do
      If QA[Z] IsIn PA Then
        Z:=Z+1;
      Else
        B:=False;
      EndIf;
    EndWhile;
  Else 
    B:=False;
  EndIf;
  Return B;
EndDefine;

-----------------------------------
-- TauPQAmmissibile(P,I,J,K,L,T) -- 
-----------------------------------

-- AIM:
-- Given the list of all vectors strictly
-- between P and Q, we construct the one
-- containing only the P-admissible
-- vectors.
--
-- INPUT:
-- P - vector corresponding to Delta_P
-- I, J, K, L - data defining the Schubert
--              variety S
-- T - list of all vectors between P and
--     Q (P, Q are excluded)
--
-- OUTPUT:
-- NT -- list of all P-admissible vectors
--       strictly containing Delta_Q and
--       different from P. NT[C,H] is the
--       H-th vector of distance from S equal 
--       to C + Sum(P) with respect to the 
--       degLex order.
--
-- SUBFUNCTIONS:
-- SchubertVariety, Actual

Define TauPQAmmissibile(P,I,J,K,L,T)
  TopLevel SchubertVariety; 
  TopLevel Actual;

  D:=len(T);
  NT:=[[] | C In 1..D];

  PA:=Actual(P,I,J,K,L);
  OmegaP:=len(PA);
  For C:=1 To D Do
    LTC:=len(T[C]);
    For H:=1 To LTC Do
      TA:=Actual(T[C,H],I,J,K,L);
      OmegaT:=len(TA);
      If SchubertVariety(T[C,H],I,J,K,L) Then
        If PQAmmissibile(PA,TA) Then
          NT[C]:=concat(NT[C],[[T[C,H],TA]]);
        EndIf;
      EndIf;
    EndFor;
  EndFor;
  Return NT;
EndDefine;

--------------------------------
-- DistanzaSigma(ADM,I,J,K,L) --
--------------------------------

-- AIM:
-- Given the list ADM consisting of P, of
-- Q and of all P-admissible vectors bewteen
-- them, we compute, whenever it makes sense,
-- the polynomials A_TX, G_TX and B_TX with T
-- and X in ADM.
--
-- INPUT:
-- ADM - 4-indices list which contains all
--       P-admissible vectors T such that
--       P <= T <= Q. ADM[H,Z] consists of
--       two lists; the latter is the vector
--       T (with respect to the essential flag
--       of P) and the former is the vector of
--       pointers TA to its essential conditions.
--       The value H-1 represents the distance between
--       T and P, whereas Z stands for the position
--       of T in the list ADM[H], which is ordered
--       with respect to the number of essential
--       conditions of the vectors.
-- I,J,K,L - data defining the Schubert variety S
-- 
-- LOCAL VARIABLES:
-- LADM - length of ADM. In particular, LADM-1 is
--        the distance between Q and P.
-- LRigheADM - list containing the lenght of ADM[H],
--             as H varies through {1, ..., LADM}.
-- LC - number of essential conditions of the first
--      element of ADM[H], for any H
-- IndiciInizio - list such that IndiciInizio[H,C] is the
--                index in ADM[H] of the first vector
--                having C essential conditions. It is
--                set to 0 if there are not such vectors
--                in ADM[H].
--
-- OUTPUT:
-- A, G, B - lists such that A[T1,T2,X1-T1,X2],
--           G[T1,T2,X1-T1,X2] and B[T1,T2,X1-T1,X2]
--           are the polynomials A_TX, G_TX and
--           B_TX with T = ADM[T1,T2] and
--           X = ADM[X1,X2].
--
-- SUBFUNCTIONS:
-- Divide, Apol, PQAmmissibile, CalcoloMTX,
--         Utilde, CalcoloXTA.

Define DistanzaSigma(ADM,I,J,K,L)
  TopLevel Divide;
  TopLevel Apol;
  TopLevel PQAmmissibile;
  TopLevel CalcoloMTX;
  TopLevel Utilde; 
  TopLevel CalcoloXTA;

  Omega:=Len(ADM[1,1,1]);
  LADM:=Len(ADM);
  LRigheADM:=[Len(ADM[H]) | H In 1..LADM];

-- Computation of IndiciInizio

  IndiciInizio:=[ [0 | Variabile In 1..Omega] | H In 1..LADM ];
  For H:=1 To LADM Do
    LC:=len(ADM[H,1,2]);
    IndiciInizio[H,LC]:=1;
    C:=2;
    While C<=LRigheADM[H] Do
      If len(ADM[H,C,2])=LC Then
        C:=C+1;
      Else
        LC:=len(ADM[H,C,2]);
        IndiciInizio[H,LC]:=C;
        C:=C+1;
      EndIf;
    EndWhile;
  EndFor;

-- Initialize A, G, B. IMPORTANT: at the end, A_TX = 0 if Delta_T
-- does not contain Delta_X; G_TX = -1 if Delta_X is not a Delta_T variety.

  A:=[ [ [ [ 0 | W In 1..LRigheADM[VV]] | VV In (H+1)..LADM  ] | Z In 1..LRigheADM[H] ] | H In 1..(LADM-1) ];

  B:=[ [ [ [ -1 | W In 1..LRigheADM[VV]] | VV In (H+1)..LADM  ] | Z In 1..LRigheADM[H] ] | H In 1..(LADM-1) ];

  G:=[ [ [ [ -1 | W In 1..LRigheADM[VV]] | VV In (H+1)..LADM  ] | Z In 1..LRigheADM[H] ] | H In 1..(LADM-1) ];

-- The following nidification of for cycles is
-- given with respect to the order of the
-- mutual distances between pairs of S-varieties

-- Computation of A

  For Sigma:=1 To LADM-1 Do
    For H:=1 To LADM-Sigma Do
      For Z:=1 To LRigheADM[H] Do
        V:=Sigma;
        For W:=1 To LRigheADM[V+H] Do
          X:=ADM[V+H,W];
          T:=ADM[H,Z];
          If Divide(T[1],X[1]) Then
            A[H,Z,V,W]:=Apol(T,I,X[1]);
          EndIf;
        EndFor;
      EndFor;
    EndFor;
  EndFor;

-- Print A

  PrintLn "A := ";
    For H := 1 To Len(A) Do
    PrintLn A[H];
    EndFor;

-- Computation of G and B

  PrintLn "Sigma := 1;";
  V:=1; 
  For H:=1 To LADM-1 Do
    For Z:=1 To LRigheADM[H] Do
      For W:=1 To LRigheADM[V+H] Do
        X:=ADM[V+H,W];
        T:=ADM[H,Z];
        PrintLn "T := ", T, "; X := ", X, ";";

-- If A[H,Z,V,W]<>0, we know that Delta_T contains
-- Delta_X. Moreover, if X is T-admissible, we can
-- compute G_TX and B_TX; otherwise, we have to
-- replace X with X^T and put B_TX = B_TX^T
-- (G_TX does not have sense in this case).
 
        If A[H,Z,V,W]<>0 Then 
          If PQAmmissibile(T[2],X[2]) Then
            MTX:=CalcoloMTX(T[1],I,J,K,L,X[1]);
            PrintLn "MTX := ", MTX, ";";
            G[H,Z,V,W]:= Utilde(MTX,A[H,Z,V,W]);
            B[H,Z,V,W]:=A[H,Z,V,W]-G[H,Z,V,W];
          Else

-- If X is not T-admissible, X^T = T. B_TX = 1
-- and G_TX does not have sense in this case.

            B[H,Z,V,W]:=1;
          EndIf;
        EndIf;
      EndFor;
    EndFor;
  EndFor;

  For Sigma:=2 To LADM-1 Do
  PrintLn "Sigma := ", Sigma, ";";
    For H:=1 To LADM-Sigma Do
      For Z:=1 To LRigheADM[H] Do
        V:=Sigma;
        For W:=1 To LRigheADM[V+H] Do
          X:=ADM[V+H,W];
          T:=ADM[H,Z];
          PrintLn "T := ", T, "; X := ", X, ";";
          If A[H,Z,V,W]<>0 Then
            If PQAmmissibile(T[2],X[2]) Then
              MTX:=CalcoloMTX(T[1],I,J,K,L,X[1]);
              PrintLn "MTX := ", MTX, ";";

-- Computation of RTX

              RTX:=A[H,Z,V,W];
              For MU:=H+1 To V+H-1 Do
                For Col:=1 To LRigheADM[MU] Do
                  Y:=ADM[MU,Col];

-- Check if Delta_T contains Delta_Y, Delta_Y 
-- contains Delta_X and Y is T-admissibile
-- (the last one is equivalent to
-- G[H,Z,MU-H,Col]<>-1).

                  If A[H,Z,MU-H,Col]<>0 And A[MU,Col,V+H-MU,W]<>0 And G[H,Z,MU-H,Col]<>-1 Then
                    RTX:=RTX-G[H,Z,MU-H,Col]*B[MU,Col,V+H-MU,W];
                  EndIf;
                EndFor;
              EndFor;
              G[H,Z,V,W]:=Utilde(MTX,RTX);
              B[H,Z,V,W]:=RTX-G[H,Z,V,W];
            Else

-- Find the essential conditions XTA:=(X^T)A of X^T

              XTA:=CalcoloXTA(T,I,J,K,L,X);

-- Find the variety [X[1],XTA] in ADM[C], with C >= H
-- because it is contained in Delta_T and contains
-- Delta_X by construction.

              NC:=Len(XTA);
              Logica:=False;
              HH:=H;
              While Not(Logica) And HH < V+H Do
                N1:=IndiciInizio[HH,NC];
                If N1>0 Then
                  If NC=Omega Then
                    N2:=LRigheADM[HH]+1;
                  Else
                    C:=1;
                    While NC+C<=Omega And IndiciInizio[HH,NC+C]=0 Do
                      C:=C+1;
                    EndWhile;
                    If NC+C>Omega Then
                      N2:=LRigheADM[HH]+1;
                    Else
                      N2:=IndiciInizio[HH,NC+C];
                    EndIf;
                  EndIf;
                  While N1<N2 And Not(Logica) Do
                    App:=ADM[HH,N1];        
                    If [ X[1,XTA[N]] | N In 1..NC] = [ App[1,XTA[N]] | N In 1..NC] Then
                      Logica:=True;
                      If HH<>H Then
                        B[H,Z,V,W]:=B[H,Z,HH-H,N1];
                      Else
                        B[H,Z,V,W]:=1;
                      EndIf;
                    EndIf;
                    N1:=N1+1;
                  EndWhile;
                EndIf;
                HH:=HH+1;
              EndWhile;
            EndIf;
          EndIf;
        EndFor;
      EndFor;
    EndFor;
  EndFor;

-- Print G and B

  PrintLn "G := ";
  For H := 1 To Len(G) Do
    PrintLn G[H];
  EndFor;

  PrintLn "B := ";
  For H := 1 To Len(B) Do
    PrintLn B[H];
  EndFor;

  Return B[1,1,LADM-1,1];
EndDefine;

-----------------------------
-- CalcoloXTA(T,I,J,K,L,X) --
-----------------------------

-- AIM:
-- Given two S-varieties Delta_T and Delta_X,
-- with X T-admissible, we describe X by means
-- of the essential flag of T.
--
-- INPUT:
-- I, J, K, L - data defining the Schubert
--              variety S
-- T, X - lists with T[1] and X[1] increments
--        with respect to S and T[2] and X[2]
--        lists of pointers to the essential
--        flags of Delta_T and Delta_X.
--
-- OUTPUT:
-- XTA - list containing the vector
--       of increments X^T (with
--       respect to S) and the pointers
--       to the essential flag of S.

Define CalcoloXTA(T,I,J,K,L,X);
  Ome:=Len(T[2]);
  If Ome>=2 Then
    XTA:=[];
    XTI:=X[1]+I;
    If XTI[T[2,2]]-XTI[T[2,1]]<J[T[2,2]]-J[T[2,1]] Then
      XTA:=concat(XTA,[T[2,1]]);
    EndIf;

    For Alpha:=2 To Ome-1 Do
      If XTI[T[2,Alpha+1]]-XTI[T[2,Alpha]]<J[T[2,Alpha+1]]-J[T[2,Alpha]] And XTI[T[2,Alpha-1]]<XTI[T[2,Alpha]] Then
        XTA:=concat(XTA,[T[2,Alpha]]);
      EndIf;        
    EndFor;

    If K+J[T[2,Ome]]<L+XTI[T[2,Ome]] And XTI[T[2,Ome-1]]<XTI[T[2,Ome]] Then
      XTA:=concat(XTA,[T[2,Ome]]);
    EndIf;     
  Else
    XTA:=[T[2,1]];
  EndIf;

  Return XTA;
EndDefine;

-------------
-- Hpol(B) --
-------------

-- AIM:
-- Given an integer B, we compute the sum of
-- even powers X[1]^2A, where X[1] is the
-- first unknown of the current ring (by
-- default, it is Q[x,y,z]).
--
-- INPUT:
-- B - integer
--
-- OUTPUT:
-- H - the polynomial described above.

Define Hpol(B)
  TopLevel CurrentRing;
  X:=indets(CurrentRing);
  If B<0 Then
    H:=0;
  Else 
    H:=sum([X[1]^(2*A) | A In 0..B]);
  EndIf;

  Return H;
EndDefine;

------------
-- Pol(B) --
------------

-- AIM:
-- Given an integer B, we compute the product
-- of all polynomials Hpol(A) (see the
-- corresponding function) as A runs through
-- {1, ..., B - 1}.
--
-- INPUT:
-- B - integer
--
-- OUTPUT:
-- P - polynomial as described above
--
-- SUBFUNCTIONS:
-- Hpol

Define Pol(B)
  TopLevel Hpol;
  If B<0 Then
    P:=0;
  Else
    P:=1;
  EndIf;

  If B>1 Then
    P:=product([Hpol(A) | A In 0..(B-1)]);
  EndIf;

  Return P;
EndDefine;

-----------------
-- Apol(T,I,X) --  
-----------------

-- AIM:
-- Given two S-varieties Delta_T and Delta_X,
-- with Delta_X contained in Delta_T, we compute
-- the polynomial A_TX of the fibre of Pi_T at
-- (any point of) Delta_X.
--
-- INPUT:
-- T - pair of lists: the former contains the vector
--     defining Delta_T, while the latter is the list 
--     of pointers to the essential conditions of Delta_T
-- I - some of the data defining the Schubert variety S
-- X - vector defining another S-variety
--
-- OUTPUT:
-- A - polynomial of the fiber F_TX
--
-- SUBFUNCTIONS:
-- Pol

Define Apol(T,I,X)
  TA:=T[2];
  OmegaT:=len(TA);
  A:=Pol(I[TA[1]]+X[TA[1]])/( Pol(I[TA[1]]+T[1,TA[1]]) * Pol(X[TA[1]]-T[1,TA[1]]) );
  If OmegaT > 1 Then
    L:=[ Pol(I[TA[Alpha]]+X[TA[Alpha]] - I[TA[Alpha-1]]-T[1,TA[Alpha-1]])/( Pol(I[TA[Alpha]]+T[1,TA[Alpha]] - I[TA[Alpha-1]]-T[1,TA[Alpha-1]]) * Pol(X[TA[Alpha]]-T[1,TA[Alpha]]) ) | Alpha In 2..OmegaT];
    A:=A*product(L);
  EndIf;

  Return A;
EndDefine;

-------------------
-- Utilde(N,F) ----
-------------------

-- AIM:
-- Given a polynomial F and an integer N,
-- we delete all monomials of degree < N
-- of F and make the resulting polynomial
-- symmetric.
--
-- INPUT:
-- N - integer
-- F - polynomial
--
-- OUTPUT:
-- H - polynomial

Define Utilde(N,F)
  TopLevel CurrentRing;
  X:=indets(CurrentRing);
  If F <> 0 Then
    M := deg(F);
  Else
    M := 0;
  EndIf;

  If M = 0 Then
    H := 0;
  Else
    If N > M Then
      H := 0;
    Else
      If N = 0 Then
        H := Eval(F, [0]);
        Else
        If N > 0 Then
          H := F;

-- Truncate H

          For T := 1 To N Do
            H := (H - eval(H, [0]))/X[1];
          EndFor;

-- Save the coefficients we need to symmetrize H

          CO := [eval(H, [0])];
          G := H;
          M := deg(H);
          For T := 1 To M Do
            G := (G - CO[T])/X[1];
            append(ref CO, eval(G, [0]));
          EndFor;

-- Now we symmetrize H

          H := H*(X[1]^N);

-- Delete the coefficient of the middle term

          remove(ref CO, 1); 
          For T := 1 To M Do
            H := H + CO[T]*(X[1]^(N - T));
          EndFor;
        EndIf;
      EndIf;
    EndIf;
  EndIf;
  Return H;
EndDefine;

--------------
-- Examples --
--------------

-- MainKaLu(I,J,K,L,I1,J1);

-- I := [1,3,5]; J := [8,12,17]; K := 7; L := 20;
-- I1 := [2,5,6]; J1 := [6, 11, 16];