with(ListTools):
with(LinearAlgebra):
with(Groebner):
truediff:=false:
kernelopts(printbytes=false):

################################################################
# This file implements several algorithms for computing linear #
# recurrence relations satisfied by a sequence.                #
################################################################

########################
# Secondary procedures #
########################

SeqVector:= proc (Tab,d)
description "Returns Tab as a vector of size d";
    return Vector[row](d,i->Tab(i-1));
end proc:

SeqMatrix:= proc (Tab,d1,d2:=d1)
# description "Returns Tab as a matrix of size d1 times d2 with the first "
#     "index increasing on the rows and the second one on the columns";
#     return Matrix(d1,d2,(i,j)->Tab(i-1,j-1));
description "Returns Tab as a matrix of size d2 times d1 with the first "
    "index increasing on the columns and the second one decreasing on the rows";
    return Matrix(d2,d1,(i,j)->Tab(j-1,d2-i));
end proc:

# Print Tab as a vector of matrix with i increasing on the columns
# j increasing on the rows and k increasing in the vector
Seq3Tensor:= proc (Tab,d1,d2:=d1,d3:=d1)
description "Returns Tab as a vector of size d3 of matrices of size d1 times d2 "
    "with the first index increasing on the rows, the second one on the columns "
    "and the third one on the vector.";
    return Vector[row] (d3,k->SeqMatrix ((i,j)->Tab(i,j,k-1),d1,d2));
end proc:

redp:= proc (x,p:=mychar)
description "Reduces x mod p if p <> 0, does nothing otherwise.";
    if (p = 0) then return expand (x);
    else return expand (x) mod p;
    end if;
end proc:

listVars:= proc (n, s:='x')
description "Makes a list of n variables from string s.";
local i:
    ASSERT (n > 0, "n must be positive");
    if (n = 1) then return [s];
    else return [seq (cat (s,i), i=1..n)];
    end if:
end proc:

orderComparison:= proc (mon_ord)
description "Makes a monomial comparison function from monomial ordering "
    "mon_ord.";
    return (a,b) -> TestOrder (a,b,mon_ord);
end proc:

sortListMon:= proc (L, mon_ord)
description "Sorts the list of monomials L for mon_ord.";
    return sort (L, orderComparison (mon_ord));
end proc:

listMonLessOrEqualThan:= proc (vars,mon_ord,max_mon,gen:=vars)
description "Makes a sorted list of all the monomials in vars up to max_mon for "
    "mon_ord. \n"
    "Parameter mon_ord should not be an elimination ordering!";
local L,LL,M,i,cmp,m,x:
    M:= Matrix (MatrixOrder (mon_ord,vars));
    for i from 1 to nops(vars) do
        ASSERT (M[1,i] > 0, "No elimination ordering!");
    end do:
    L  := [];
    LL := [1];
    cmp:= orderComparison (mon_ord);
    while (LL <> {}) do
        L := {op(L),op(LL)};
        LL:= {seq (seq (`if` (cmp (m*x,max_mon),m*x,NULL), x in gen),
                   m in LL)}:
    end do:
    return sortListMon ([op(L)], mon_ord);
end proc:

listMonOfDegreeLessOrEqualThan:= proc (vars,mon_ord,max_deg,gen:=vars)
description "Makes a sorted list of all the monomials in vars of degree "
    "up to max_deg for mon_ord.";
local L,LL,i,cmp,m,x:
    L  := [];
    LL := [1];
    while (LL <> {}) do
        L := {op(L),op(LL)};
        LL:= {seq (seq (`if` (degree (m,vars)<max_deg,m*x,NULL), x in gen),
                   m in LL)}:
    end do:
    return sortListMon ([op(L)], mon_ord);
end proc:

numberOfMonomials:= proc (n,d)
description "Returns the number of monomials of degree at most d in n variables";
    return binomial (n+d,d);
end proc:

# listMonDividing:= proc (vars,mon_ord,max_mon)
# description "Makes a sorted list of all the monomials in vars dividing max_mon for "
#     "mon_ord.";
# local L,vars2,v,W:
#     vars2:= sortListMon (vars, mon_ord);
#     L := [1];
#     for v in vars2 do
#         W:= [seq (v^d,d=1..degree (max_mon,v))];
#         L:= [op(L), seq (seq (l*w, w in W), l in L)];
#     end do:
#     return sortListMon (L, mon_ord);
# end proc:

mergeOrderedLists:= proc (L1,L2,mon_order)
description "Merges two ordered lists L1 and L2 without repetition wrt the ordering "
    "given by the monomial ordering mon_order.";
local L,LL,i,j;
    L:= [];
    i:= 1; j:= 1;
    while (i <= nops (L1) and j <= nops (L2)) do
        if (orderComparison (mon_order) (L1[i],L2[j])) then
            L:=[op(L),L1[i]];
            if (L1[i] = L2[j]) then j:= j+1;
            end if:
            i:= i+1;
        else
            L:=[op(L),L2[j]];
            j:= j+1;
        end if:
    end do:
    if (i > nops (L1))
    then L:= [op(L),op(j..-1,L2)];
    else L:= [op(L),op(i..-1,L1)];
    end if:
    return L;
    # LL:= [L[1]];
    # i:= 2;
    # while (i <= nops (L)) do
    #     if (L[i] <> L[i-1]) then
    #         LL:= [op(LL),L[i]];
    #     end if;
    #     i:= i+1;
    # end do;
    # return LL;
end proc:


MinkowskiSum:= proc (mon_ord,set_mon1,set_mon2:=set_mon1)
description "Makes the sorted for mon_ord MinkowskiSum of the two sets of terms "
    "set_mon1 and set_mon2.";
local L,t,u;
    L:= {seq (seq (t*u, u in set_mon2), t in set_mon1)};
    return sortListMon ([op(L)], mon_ord);
end proc:

isGb:= proc (G,mon_ord,p:=mychar)
description "Tests if G is a reduced Gröbner basis for mon_ord in characteristic "
    "p.";
local H;
    H:= {op(Basis (G,mon_ord,characteristic=p))};
    if (p <> 0)
    then
        return ({op(map(g->g/LeadingCoefficient (g,mon_mon_ord) mod p,G))} =
                {op(map(g->g/LeadingCoefficient (g,mon_mon_ord) mod p,H))});
    else
        return ({op(map(g->g/LeadingCoefficient (g,mon_mon_ord),G))} =
                {op(map(g->g/LeadingCoefficient (g,mon_mon_ord),H))});
    end if:
end proc:

firstIndependentRows:= proc (M,p:=mychar,max_rank:=min(Dimension(M)))
option remember:
description "Returns the first independent rows of M with a naive "
    "complexity in O(n^3).";
local N, i, j, k, col_pivot, inv_pivot, m, n, r, rows;
    userinfo (1,firstIndependentRows,`*** Matrix`, convert (M,string),
              `in charistic`, p, `***`);
    m,n:= Dimension (M);
    N:= redp (Matrix (M),p);
    col_pivot:= Array ([seq (0, i=1..m)]);
    r:= 0; i:=1;
    printf ("0%%/");
    while (r <> max_rank and i <= m) do
        if (i mod 100 = 0) then
            printf ("%.1f%%/",100*i/m);
        end if:
        j:= 1;
        while (j <= n and N[i,j] = 0) do j:= j+1;
        end do:
        if (j <= n) then
            userinfo (2,firstIndependentRows,`***\t New independent row`, i,
                      `***`);
            col_pivot[i]:= j; inv_pivot:= redp (N[i,j]^(-1),p);
            r:= r+1;
            for k from i+1 to m do
                if (N[k,j] <> 0) then
                    N[k,j..n]:= redp (N[k,j..n]-(inv_pivot*N[k,j])*N[i,j..n],p);
                end if:
            end do:
        end if:
        i:= i+1;
    end do:
    printf ("\n");
    rows:= Array ([seq (`if` (col_pivot[i]=0,NULL,i), i=1..m)]);
    userinfo (1,firstIndependentRows,`*** Rank`, r, `with rows`,
              convert (rows,list));
    return rows,r;
end proc:

firstIndependentColumns:= proc (M,p:=mychar,max_rank:=min(Dimension(M)))
description "Returns the first independent columns of M with a naive "
    "complexity in O(n^3).";
option remember:
    return firstIndependentRows (Transpose (M),p,max_rank);
end proc:

staircaseStabilization:= proc (S,vars,mon_ord)
description "Returns the full minimal staircase containing S, "
    "i.e. the set of all divisors of elments in S, sorted for mon_ord.";
local Sp, s, T, v,t,d;
    Sp:= {};
    for s in S do
        T:= {1};
        for v in vars do
            T:= {seq (seq (t*v^d, d=0..degree (s,v)), t in T)};
        end do:
        Sp:= Sp union T;
    end do:
    Sp:= [op(Sp)];
    return sortListMon (Sp, mon_ord);
end proc:
sparseStaircaseStabilization:= proc (S,set_mon,lattice_mon,mon_ord)
description "Returns the full minimal staircase containing S, "
    "i.e. the set of all divisors of elments in S, sorted for mon_ord.";
local Sp, s, T, m,q;
    Sp:= {};
    for s in S do
        T:= {op(select (m->divide (s,m,'q') and q in lattice_mon,set_mon))};
        Sp:= Sp union T;
    end do:
    Sp:= [op(Sp)];
    return sortListMon (Sp, mon_ord);
end proc:

staircaseMaximum:= proc (S,vars,mon_ord)
description "Returns the maximal for | elements of the staircase S.";
local maxStaircase, maxSpan,s,t,isMax,i,cmp,q;
    # print("S=",S);
    maxStaircase:= [];
    maxSpan     := [];
    cmp         := orderComparison (mon_ord);
    for s in S do
        isMax:= true;
        for t in S do
            if (s[2] <> 0 and divide (t[2],s[2],'q') and q <> 1) then
                isMax:= false; break;
            end if:
        end do:
        if isMax then
            if not s[2] in maxSpan then
                maxStaircase:= [op(maxStaircase), s];
                maxSpan     := [op(maxSpan), s[2]];
            # else
            #     for i from 1 to nops(maxStaircase) do
            #         if (cmp(s[1],maxStaircase[i][1])) then
            #             maxStaircase[i]:=s; break;
            #         end if:
            #     end do:
            end if:
        end if:
    end do:
    return maxStaircase;
end proc:

fromStaircaseToGb:= proc (S,vars,mon_ord)
local G,stabilizedS,s,v,m,x,g;
    if (S = []) then return [1];
    end if:
    stabilizedS:= staircaseStabilization (S,vars,mon_ord);
    G:= [op({seq (seq (v*s, v in vars), s in stabilizedS)} minus {op(stabilizedS)})];
    return InterReduce (G,mon_ord);
end proc:

fromGbToStaircase:= proc (G,vars,mon_ord)
local L,LL,i,cmp,m,x,g:
    L  := [];
    LL := [1];
    cmp:= orderComparison (mon_ord);
    while (LL <> {}) do
        L := {op(L),op(LL)};
        LL:= {seq (seq (`if` (`or` (seq (divide (m*x,g), g in G)),NULL,m*x), x in vars),
                   m in LL)}:
    end do:
    return sortListMon ([op(L)], mon_ord);
end proc:

getRows_mon:= proc (all_mon,cols_mon)
description "Returns the maximal ordered set U such that cols_mon "
    "+ U is included in all_mon";
local rows_mon, AM;
    AM      := {op(all_mon)};
    rows_mon:= select (m->{op(expand (m*cols_mon))} subset AM,all_mon);
    return rows_mon;
end proc:

# isRelation:= proc (R,all_mon,cols_mon,lcm_mon,mon_ord)
# local rows_mon;
#     rows_mon:= getRows_mon (all_mon,cols_mon);
#     if (rows_mon = []) then return true;
#     else return not orderComparison (mon_ord)(lcm_mon/rows_mon[-1],R[1]);
#     end if:
# end proc:

nfLeftPart:= proc (R,vars,l,mon_ord,p:=mychar)
local tmp,q,i,r,L;
    r:= redp (expand (R),p):
    L:= redp (expand (l),p):
    if (L = []) then
        if (r[2] <> 0) then return redp (expand (r/LeadingCoefficient (r[2],mon_ord)),p);
        else                return redp (expand (r),p);
        end if:
    end if:
    NormalForm (r[1],map (l->l[1],L),mon_ord,'q',characteristic=p):
    tmp:= redp (expand (r - add (q[i]*L[i],i=1..nops(L))),p);
    if (tmp[2] <> 0) then
        tmp:=redp (expand (tmp/LeadingCoefficient (tmp[2],mon_ord)),p);
    end if:
    try:
        tmp[1]:= sort (tmp[1],order=mon_ord);
        tmp[2]:= sort (tmp[2],order=mon_ord);
    catch:
    end try:
    return tmp;
end proc:

nfRightPart:= proc (R,vars,L,mon_ord,p:=mychar)
local res,r;
    r:= redp (expand (R),p):
    res:= nfLeftPart ([r[2],r[1]],vars,map (l->[l[3],l[2]],L),mon_ord,p);
    return [res[2],res[1]];
end proc:

leftHigherPart:= proc (R,vars,M,G,mon_ord,p)
description "Remove all the monomials of R that can divide monomials M/g for g in G";
local tmp,r;
    r:= redp (expand (R),p):
    if (G = []) then return r;
    end if:
    tmp:= expand (mirrorPolynomial (r,vars,M));
    tmp:= NormalForm (tmp,G,mon_ord);
    return expand (mirrorPolynomial (tmp,vars,M));
end proc:

nfLeftHigherPart:= proc (R,vars,L,M,G,mon_ord,p)
local tmp,rHP,LHP,q,i,r;
    r:= redp (expand (R),p):
    if (L = []) then
        if (r[2] <> 0) then return redp (expand (r/LeadingCoefficient (r[2],mon_ord)),p);
        else                return redp (expand (r),p);
        end if:
    end if:
    rHP:= [leftHigherPart (r[1],vars,M,G,mon_ord,p),r[2]];
    LHP:= map (l->[`if`(l[2]=0,l[1],leftHigherPart (l[1],vars,M,G,mon_ord,p)),
                   l[2]],L);
    NormalForm (expand (rHP)[1],map (l->l[1],expand (LHP)),mon_ord,'q',characteristic=p):
    tmp:= redp (expand (r - add (q[i]*L[i],i=1..nops(L))),p);
    if (tmp[2] <> 0) then
        tmp:=redp (expand (tmp/LeadingCoefficient (tmp[2],mon_ord)),p);
    end if:
    try:
        tmp[1]:= sort (tmp[1],order=mon_ord);
        tmp[2]:= sort (tmp[2],order=mon_ord);
    catch:
    end try:
    return tmp;
end proc:

# return 1st  monomial s such that s*R fails, -1 if no such monomial exists.
# return last monomial s such that s*R has been tested.
FailingWhenShiftedBy:= proc (R,all_mon,lcm_mon,mon_ord,B,rows_mon)
local xvar, cols_mon, lm1, lm2, shift_rows,i,j,c,m;
    xvar      := [op(indets (lcm_mon))];
    lm1       := LeadingMonomial (R[1],mon_ord);
    lm2       := LeadingMonomial (R[2],mon_ord);
    coeffs (R[2],xvar,'cols_mon');
    cols_mon  := [cols_mon];
    # shift_rows:= map (a->a/lm2,select (b->divide (b,lm2),all_mon));
    shift_rows:= map (a->ifelse (divide (a,lm2,'q'),q,NULL),all_mon);
    shift_rows:= select (a->map (b->b*a,{op(cols_mon)}) subset {op(all_mon)},
                         shift_rows);
    shift_rows:= sort ([op({op(shift_rows)} union {op(rows_mon)})],
                       (a,b)->TestOrder (a,b,mon_ord));
    shift_rows:= [0,op(shift_rows)];
    for i from 2 to nops(shift_rows) do
        m:= lcm_mon/shift_rows[i];
        c:= R[1];
        for j from 1 to nops(xvar) do c:= coeff (c,xvar[j],degree (m,xvar[j]));
        end do:
        if (c <> 0) then
            return shift_rows[i],shift_rows[-1];
        end if:
    end do:
    # printf ("%a\n", rows[-1]);
    return -1,shift_rows[-1];
end proc:

isInTheCone:= proc (m,vars,cone_gen)
option remember,q;
local g;
    if (degree (m,vars) = 0) then return true;
    end if:
    for g in cone_gen do
        if (divide (m,g,'q') and isInTheCone (q,vars,cone_gen)) then
            return true;
        end if:
    end do:
    return false;
end proc:


####################
# Table procedures #
####################

fromMonToTabElem:= proc (Tab,vars,m)
description "Returns the table element of Tab corresponding to monomial "
    "m.\n"
    "For instance, if m=x^i*y^j and vars=[x,y,z], it will return Tab(i,j,0).";
    local v:
    return Tab(seq (degree (m,v), v in vars));
end proc:

fromPolToTabRelation:= proc (Tab,vars,P,p)
description "Evaluates the polynomial P at the table Tab by "
    "replacing every monomial by its corresponding table element.\n"
    "See fromMonToTabElem.";
local c,m,size,i;
    c:= [coeffs (expand (P),vars,'m')];
    m:= [m];
    size:= nops(c);
    return redp (add (c[i]*fromMonToTabElem (Tab,vars,m[i]), i=1..size),p);
end proc:

testRelation:= proc (Tab,vars,mon_ord,P,shift_mon,p:=mychar,errors:=false)
description "Evaluates P and its shifts for the monomial ordering mon_ord, "
    "at table Tab, allowing the user to test P as a valid relation.\n"
    "If P fails, errors set to true allows the user to get the failing "
    "shifting monomials.\n"
    "See fromPolToTabRelation.";
local s;
    s:= map (u->`if`(fromPolToTabRelation (Tab,vars,P*u,p) <> 0,u,NULL),
             shift_mon);
    return evalb (s = []),`if`(errors,s,NULL);
end proc:

testGb:= proc (Tab,vars,mon_ord,G,shift_mon,p:=mychar,errors:=false)
description "Evaluates the polynomials of G and their shifts "
    "for monomial ordering mon_ord, at table Tab, allowing the user "
    "to test G as a valid "
    "Gröbner basis of relations.\n"
    "If G fails, errors set to true allows the user to get the failing "
    "shifting monomials.\n"
    "See testRelation.";
local B,b;
    B:= map (g->[testRelation (Tab,vars,mon_ord,g,shift_mon,p,errors)], G);
    for b in B do
        if (not b[1]) then return false,`if`(errors,map(m->m[2],B),NULL);
        end if:
    end do:
    return true,`if`(errors,[],NULL);
end proc:

from2SetsToMultiHankelMatrix:= proc (Tab,vars,rows_mon,cols_mon:=rows_mon)
description "Makes the multi-Hankel matrix of Tab for sets of terms cols_mon "
    "and rows_mon.";
    return Matrix (nops(rows_mon), nops(cols_mon),
                   (i,j)-> fromMonToTabElem (Tab,vars,cols_mon[j]*rows_mon[i]));
end proc:

mirrorPolynomial:= proc (P,vars,M)
description "Makes the mirror of polynomial P in variables vars.\n"
    "M is the monomial mapped to 1 and should be a multiple of the lcm of "
    "the support of P for the result to be a polynomial.";
    return expand (M * subs (map (v->v=1/v, vars),P));
end proc:

from1SetToTruncatedGeneratingSeries:= proc (Tab,vars,set_mon)
description "Makes the truncated generating series of Tab for set of terms set_mon.";
local t:
    return add (fromMonToTabElem (Tab,vars,t)*t, t in set_mon);
end proc:

from1SetToMirrorTruncatedGeneratingSeries:= proc (Tab,vars,set_mon)
description "Makes the miror of the truncated generating series of Tab for "
    "set of terms set_mon.";
    return mirrorPolynomial (from1SetToTruncatedGeneratingSeries (Tab,vars,set_mon),vars,
                             lcm(op(set_mon)));
end proc:

from2SetsToTruncatedGeneratingSeries:= proc (Tab,vars,set_mon1,set_mon2:=set_mon1)
description "Makes the truncated generating series of Tab for set of terms "
    "set_mon1 + set_mon2.";
local all_mon,u,t;
    all_mon:= [op({seq (seq (t*u, t in set_mon1), u in set_mon2)})];
    return from1SetToTruncatedGeneratingSeries (Tab,vars,all_mon);
end proc:

from2SetsToMirrorTruncatedGeneratingSeries:= proc (Tab,vars,set_mon1,set_mon2:=set_mon1)
description "Makes the miror of the truncated generating series of Tab for "
    "set of terms set_mon1 + set_mon2.";
    return mirrorPolynomial (from2SetsToTruncatedGeneratingSeries (Tab,vars,
                                                                   set_mon1,set_mon2),
                             vars,lcm(op(set_mon1))*lcm(op(set_mon2)));
end proc:

fromGbToRandomSequence:= proc (G,vars,mon_order,p:=mychar)
description "Returns a random table in characteristic p"
    "whose ideal of relations is <G>"
    "if <G> is Gorenstein or contains <G> otherwise.";
local Tab,lmG,remG,roll,i;
    lmG := map (g->Groebner:-LeadingMonomial (g,mon_order), G);
    remG:= [seq (redp (lmG[i]-G[i], p),i=1..nops(G))];
    roll:= `if`(p=0,rand(-2^16..2^16-1),rand(0..p-1));
    Tab := proc ()
    option remember;
    local idx,m,i;
        idx:= [_passed];
        m  := mul (vars[i]^idx[i],i=1..nops(idx));
        for i from 1 to nops (lmG) do
            if divide (m,lmG[i],'q') then
                return fromPolToTabRelation (thisproc,vars,q*remG[i],p);
            end if:
        end do:
        return roll ();
    end proc;
    return Tab;
end proc:

fromGbToExtendedTab:= proc (G,vars,mon_order,initTab,p:=mychar)
description "Returns a table extending initTab in characteristic p"
    "whose ideal of relations is <G>"
    "if <G> is Gorenstein or contains <G> otherwise.";
local Tab,lmG,remG,roll,i;
    lmG := map (g->Groebner:-LeadingMonomial (g,mon_order), G);
    remG:= [seq (redp (lmG[i]-G[i], p),i=1..nops(G))];
    Tab := proc ()
    option remember;
    local idx,m,i;
        idx:= [_passed];
        m  := mul (vars[i]^idx[i],i=1..nops(idx));
        for i from 1 to nops (lmG) do
            if divide (m,lmG[i],'q') then
                return fromPolToTabRelation (thisproc,vars,q*remG[i],p);
            end if:
        end do:
        return initTab (_passed);
    end proc;
    return Tab;
end proc:


#####################
# Table procedures  #
# Quasi-commutative #
#####################

fromQCMonToTabElem:= proc (Tab,vars,diffvars,m,p:=mychar)
description "Returns the table element of Tab corresponding to monomial "
    "m.\n"
    "For instance, if m=t1^j1*t2^j2*x1^i1*x2^i2, vars=[x1,x2,x3] "
    "and diffvars=[t1,t2,t3] "
    "it will return i1^j1*i2^j2*Tab(i1,i2,0).";
    local v,i:
    ASSERT (nops (vars) = nops (diffvars), "the numbers of variables must "
            "be the same");
    return redp (mul (degree (m,vars[i])^(degree (m,diffvars[i])),i=1..nops(vars))
                 *Tab(seq (degree (m,v), v in vars)),p);
end proc:

fromQCPolToTabRelation:= proc (Tab,vars,diffvars,P,p:=mychar)
description "Evaluates the polynomial P at the table Tab by "
    "replacing every monomial by its corresponding table element.\n"
    "See fromQCMonToTabElem.";
local c,m,size,i;
    c:= [coeffs (expand (P),[op(vars),op(diffvars)],'m')];
    m:= [m];
    size:= nops(c);
    return redp (add (c[i]*fromQCMonToTabElem (Tab,vars,diffvars,m[i]), i=1..size),p);
end proc:


testQCRelation:= proc (Tab,vars,diffvars,mon_ord,P,shift_mon,p:=mychar,errors:=false)
description "Evaluates P and its shifts for the monomial ordering mon_ord, "
    "at table Tab, allowing the user to test P as a valid relation.\n"
    "If P fails, errors set to true allows the user to get the failing "
    "shifting monomials.\n"
    "See fromPolToTabRelation.";
local s;
    s:= map (u->`if`(fromQCPolToTabRelation (Tab,vars,diffvars,P*u,p) <> 0,u,
                     NULL), shift_mon);
    return evalb (s = []),`if`(errors,s,NULL);
end proc:

testQCGb:= proc (Tab,vars,diffvars,mon_ord,G,shift_mon,p:=mychar,errors:=false)
description "Evaluates the polynomials of G and their shifts, up to max_shift "
    "for monomial ordering mon_ord, at table T, allowing the user to test G as a valid "
    "Gröbner basis of relations.\n"
    "If G fails, errors set to true allows the user to get the failing "
    "shifting monomials.\n"
    "See testRelation.";
local B,b;
    B:= map (g->[testQCRelation (Tab,vars,diffvars,mon_ord,g,shift_mon,p,errors)], G);
    for b in B do
        if (not b[1]) then return false,`if`(errors,map(m->m[2],B),NULL);
        end if:
    end do:
    return true,`if`(errors,[],NULL);
end proc:

from2QCSetsToMultiHankelMatrix:= proc (Tab,vars,diffvars,rows_mon,cols_mon:=rows_mon,
                                       p:=mychar)
description "Makes the multi-Hankel matrix of Tab for sets of terms cols_mon "
    "and rows_mon.";
    return Matrix (nops(rows_mon), nops(cols_mon),
                   (i,j)-> fromQCMonToTabElem (Tab,vars,diffvars,
                                               cols_mon[j]*rows_mon[i],p));
end proc:

precursiveCriterion:= proc (S,vars,diffvars,mon_ord)
local Sp, s, T, mon, U, u, i, deg, diffdeg,diffmon;
    Sp:= S;
    for s in S do
        T:= staircaseStabilization ([coeffs (s, vars)],diffvars,mon_ord);
        mon:= coeffs (s, diffvars);
        T:= {seq (diffmon*mon, diffmon in T)};
        if (not (T subset {op(S)})) then Sp:= remove (sp->sp=s,Sp);
        end if:
    end do:
    return Sp,nops(Sp);
end proc:


###################
# Main procedures #
###################

NaiveBerlekampMassey:= proc (Tab,vars:='x',max_mon:=vars^10,p:=mychar,full_answer:=false)
description "Compute the minimal relation satisfied by Tab for all the monomials "
    "up to max_mon in characteristic p.\n"
    "If full_answer is set to true, returns a triplet [m,r,s] where\n"
    "- m is a monomial;\n"
    "- r is a relation with leading monomial m;\n"
    "- s is the shift of the relation, that is up to which monomial "
    "it has been tested.\n"
    "Otherwise returns r only.";
local R,shift,F,fail,mon_ord,list_mon,discrepancy,cmp,mon;
    R:= 1; shift:= 0; F:= 0; fail:= vars^(-1);
    mon_ord := tdeg(vars);
    list_mon:= listMonLessOrEqualThan ([vars],mon_ord,max_mon);
    cmp     := orderComparison (mon_ord);
    userinfo (1,BerlekampMassey,`*** Tab `, Tab,
              `in variable`, vars, `up to`, max_mon,
              `in characteristic`, p, `***`);
    userinfo (2,BerlekampMassey,`***\t New potential relation`, 1, `***`);
    for mon in list_mon do
        userinfo (1,BerlekampMassey,mon);
        shift      := mon/LeadingMonomial (R,mon_ord);
        discrepancy:= redp (fromPolToTabRelation (Tab,vars,expand (R*shift),p),p);
        if (discrepancy <> 0) then
            userinfo (2,BerlekampMassey,`***\t FAILS with a shift`,shift,`***`);
            if (cmp (shift,fail)) then
                R:= redp (expand(R-discrepancy*fail*F/shift),p);
                userinfo (2,BerlekampMassey,`***\t\t but tweaked to`, R, `***`);
            else
                userinfo (2,BerlekampMassey,`***\t\t and is discarded ***`);
                R,F       := redp (expand(shift*R/fail-discrepancy*F),p),
                             redp (expand(R/discrepancy),p);
                fail,shift:= shift,fail;
                if (degree (shift,vars) < 0) then shift:= 0;
                end if:
                userinfo (2,BerlekampMassey,`***\t New potential relation`, R,
                          `***`);
            end if:
        else
            userinfo (2,BerlekampMassey,`***\t SUCCEEDS with a shift`,shift,`***`);
        end if:
    end do:
    if (full_answer) then
        return [LeadingMonomial (R,mon_ord),R,shift];
    else
        return R;
    end if:
end proc:

PolynomialBerlekampMassey:= proc (Tab,vars:='x',max_mon:=vars^10,p:=mychar,
                                  full_answer:=false)
description "Compute the minimal relation satisfied by Tab for all the monomials "
    "up to max_mon in characteristic p.\n"
    "If full_answer is set to true, returns a triplet [m,r,s] where\n"
    "- m is a monomial;\n"
    "- r is a relation with leading monomial m;\n"
    "- s is the shift of the relation, that is up to which monomial "
    "it has been tested.\n"
    "Otherwise returns r only.";
local R0,R1,R2,Q,l,mon_ord;
    userinfo (1,BerlekampMassey,`*** Tab `, Tab,
              `in variable`, vars, `up to`, max_mon,
              `in characteristic`, p, `***`);
    mon_ord:= tdeg(vars);
    R0:= [vars*max_mon,0];
    R1:= [from1SetToMirrorTruncatedGeneratingSeries
          (Tab,[vars],listMonLessOrEqualThan ([vars],mon_ord,max_mon)),1];
    userinfo (2,BerlekampMassey,`***\t New potential relation`, 1, `***`);
    while (degree (R1[1],vars) >= degree (R1[2],vars)) do
        userinfo (2,BerlekampMassey,`***\t FAILS with shift`,
                  max_mon/LeadingMonomial (R1[1],mon_ord),`***`);
        userinfo (2,BerlekampMassey,`***\t\t and is discarded ***`);
        if (p <> 0) then Q:= Quo (R0[1],R1[1],vars) mod p;
        else             Q:= quo (R0[1],R1[1],vars);
        end if:
        R0,R1:= R1,redp (expand (R0-Q*R1),p);
        R1   := redp (expand (R1/LeadingCoefficient (R1[2],mon_ord)),p);
        userinfo (2,BerlekampMassey,`***\t New potential relation`, R1[2], `***`);
    end do:
    userinfo (2,BerlekampMassey,`***\t SUCCEEDS! ***`);
    if (full_answer) then
        return [LeadingMonomial (R1[2],mon_ord),R1[2],
                max_mon/LeadingMonomial (R1[2],mon_ord)];
    else
        return R1[2];
    end if:
end proc:

MatrixBerlekampMassey:= proc (Tab,vars:='x',
                              cols_mon:=listMonLessOrEqualThan ([vars],tdeg(vars),
                                                                vars^10),
                              rows_mon:=cols_mon,p:=mychar,full_answer:=false)
description "Compute the minimal relation satisfied by Tab for all the monomials "
    "in cols_mon + rows_mon in characteristic p."
    "If full_answer is set to true, returns a triplet [m,r,s] where\n"
    "- m is a monomial;\n"
    "- r is a relation with leading monomial m;\n"
    "- s is the shift of the relation, that is up to which monomial "
    "it has been tested.\n"
    "Otherwise returns r only.";
local H,mon_ord,firstIndRows,usefulStaircase,fullStaircase,firstIndCols;
local subH,next_mon,R,rank;
    userinfo (1,BerlekampMassey,`*** Tab `, Tab,
              `in variable`, vars, `in`, nops(cols_mon), `times`, nops(rows_mon),
              `in characteristic`, p, `***`);
    mon_ord          := tdeg(vars);
    H                := from2SetsToMultiHankelMatrix (Tab,[vars],rows_mon,cols_mon);
    firstIndRows,rank:= firstIndependentRows (H,p,min(Dimension(H)));
    if (rank = 0) then
        if (full_answer) then
            return [1,1,rows_mon];
        else
            return 1;
        end if:
    end if:
    if (rows_mon = cols_mon) then
        firstIndCols:= firstIndRows;
    else
        firstIndCols,rank:= firstIndependentColumns (H,p,rank);
    end if:
    firstIndRows     := convert (firstIndRows,list);
    firstIndCols     := convert (firstIndCols,list);
    usefulStaircase:= map (m->cols_mon[m], firstIndCols);
    userinfo (2,BerlekampMassey,`Useful staircase`, usefulStaircase);
    fullStaircase  := staircaseStabilization (usefulStaircase,[vars],mon_ord);
    userinfo (2,BerlekampMassey,`Staircase`, fullStaircase);
    next_mon:= fullStaircase[-1]*vars;
    subH    := (<H[firstIndRows,firstIndCols]|
                from2SetsToMultiHankelMatrix (Tab,[vars],
                                              map (i->rows_mon[i],firstIndRows),
                                              [next_mon])>);
    if (p = 0) then
        R:= LinearSolve (subH);
    elif (p < 2^32) then
        subH:= LinearAlgebra[Modular][Mod] (p,subH,integer[8]);
        R   := LinearAlgebra[Modular][LinearSolve] (p,subH,1,inplace=false);
    else
        R   := LinearSolve (subH) mod p;
    end if:
    R       := redp (next_mon - DotProduct (R,usefulStaircase),p);
    userinfo (2,BerlekampMassey,`***\t New potential relation`, R, `***`);
    if (full_answer) then
        return [next_mon, R, rows_mon];
    else
        return R;
    end if:
end proc:

BerlekampMasseySakata:= proc (Tab,vars:=['x','y'],mon_ord:=tdeg(op(vars)),
                              max_mon:=vars[1]^10,p:=mychar,full_answer:=false)
description "Compute the minimal relations satisfied by Tab for all the monomials "
    "up to max_mon in characteristic p.\n"
    "If full_answer is set to true, returns a set of triplets [m,r,s] where\n"
    "- m is a monomial;\n"
    "- r is a relation with leading monomial m;\n"
    "- s is the shift of the relation, that is up to which monomial "
    "it has been tested.\n"
    "Otherwise returns a set of r only.";
local R,shift,F,fail,list_mon,toTest,discrepancy,cmp,mon,staircase,newStaircase,
    Gb,g,r,i,j,k;
    R:= [[1,1,0]]; F:=[];
    list_mon:= listMonLessOrEqualThan (vars,mon_ord,max_mon);
    cmp     := orderComparison (mon_ord);
    userinfo (1,BerlekampMasseySakata,`*** Tab `, Tab,
              `in variables`, vars, `up to`, max_mon,
              `in characteristic`, p, `***`);
    userinfo (2,BerlekampMasseySakata,`***\t New potential relation`, 1, `***`);
    for mon in list_mon do
        userinfo (1,BerlekampMasseySakata,mon);
        newStaircase:= F;
        toTest     := select (r->divide (mon,r[1]),R);
        shift      := map (r->mon/r[1],toTest);
        discrepancy:= [seq (redp (fromPolToTabRelation (Tab,vars,
                                                        expand (toTest[i][2]*shift[i]),
                                                        p),
                                  p),i=1..nops(toTest))];
        for i from 1 to nops (toTest) do
            if (discrepancy[i] <> 0) then
                newStaircase:= [op(newStaircase),[redp (toTest[i][2]/discrepancy[i],p),
                                                  shift[i]]];
            else
                R[i][3]:= shift[i];
            end if:
        end do:
        newStaircase:= staircaseMaximum (newStaircase,vars,mon_ord);
        Gb          := fromStaircaseToGb (map (s->s[2],newStaircase),vars,mon_ord);
        for i from 1 to nops(Gb) do
            g:= Gb[i];
            Gb[i]:= [g,g,0];
            for r in R do
                if (divide (g,r[1],'quogr')) then break;
                end if:
            end do:
            if (not divide (mon,g,'quomong')) then
                Gb[i][2]:= expand (quogr*r[2]);
                Gb[i][3]:= 0;
                j:= 0;
                while (cmp (list_mon[j+1]*g,mon)) do
                    j:= j+1;
                end do:
                if (j <> 0) then Gb[i][3]:= list_mon[j];
                end if:
            else
                for j from 1 to nops (F) do
                    if (divide (F[j][2],quomong,'quofailquomong')) then
                        Gb[i][2]:= redp (expand (quogr*r[2]
                                                 -fromPolToTabRelation (Tab,vars,
                                                                        r[2]*mon/r[1],
                                                                        p)*
                                                 quofailquomong*F[j][1]),p);
                        Gb[i][3]:= quomong;
                        break;
                    end if:
                end do:
                if (j = nops(F)+1) then
                    Gb[i][2]:= r[2];
                    Gb[i][3]:= quomong;
                end if:
            end if:
        end do:
        i:= 1; j:= 1; k:=1;
        while (i <= nops (toTest) and j <= nops (Gb)) do
            if (cmp (toTest[i][1],Gb[j][1])) then
                userinfo (2,BerlekampMasseySakata,`***\t`,toTest[i][2],`***`);
                userinfo (2,BerlekampMasseySakata,
                          `if`(discrepancy[i] = 0,
                               `***\t SUCCEEDS with a shift`,
                               `***\t FAILS with a shift`),
                          shift[i],`***`);
                if (toTest[i][1] = Gb[j][1]) then
                    if (toTest[i][2] <> Gb[j][2]) then
                        userinfo (2,BerlekampMasseySakata,
                                  `***\t\t but tweaked to`, Gb[j][2],`***`);
                    end if:
                    j:= j+1;
                else
                    userinfo (2,BerlekampMasseySakata,
                              `***\t\t and is discarded ***`);
                end if:
                i:= i+1;
                k:= k+1;
            else
                while (k <= nops (R) and not cmp (Gb[j][1],R[k][1])) do
                    k:= k+1;
                end do:
                if (Gb[j][1] <> R[k][1]) then
                    userinfo (2,BerlekampMasseySakata,`***\t New potential relation in`,
                              Gb[j][1],`\n`,Gb[j][2],`***`);
                end if:
                j:= j+1;
            end if:
        end do:
        while (i <= nops (toTest)) do
            userinfo (2,BerlekampMasseySakata,`***\t`,toTest[i][2],`***`);
            userinfo (2,BerlekampMasseySakata,
                      `if`(discrepancy[i] = 0,
                           `***\t SUCCEEDS with a shift`,
                           `***\t FAILS with a shift`),
                      shift[i],`***`);
            i:= i+1;
        end do:
        while (j <= nops (Gb)) do
            while (k <= nops (R) and not cmp (Gb[j][1],R[k][1])) do
                k:= k+1;
            end do:
            if (k > nops (R) or Gb[j][1] <> R[k][1]) then
                userinfo (2,BerlekampMasseySakata,`***\t New potential relation in`,
                          Gb[j][1],`\n`,Gb[j][2],`***`);
            end if:
            j:= j+1;
        end do:
        R:= Gb;
        F:= newStaircase;
    end do:
    if (full_answer) then
        return R;
    else
        return InterReduce (map (g->g[2],R),mon_ord,characteristic=p);
    end if:
end proc:

BerlekampMasseySakata_queries:= proc (Tab,vars,mon_ord,max_mon)
description "Compute the number of table queries performed by "
    "BerlekampMasseySakata to recover the minimal relations satisfied by Tab "
    "for all the monomials up to max_mon.";
local list_mon;
    list_mon:= listMonLessOrEqualThan (vars,mon_ord,max_mon);
    return nops (list_mon);
end proc:

BerlekampMasseySakata_basicop:= proc (Tab,vars,mon_ord,Gb_LM,staircase,max_mon)
description "Compute the number of basic operations performed by "
    "BerlekampMasseySakata to recover the minimal relations satisfied by Tab "
    "for all the monomials up to max_mon.";
local list_mon,cmp,basicop,current_staircase,current_Gb;
local new_staircase,new_Gb,mon,g,q,i;
    cmp     := orderComparison (mon_ord);
    list_mon:= listMonLessOrEqualThan (vars,mon_ord,max_mon);
    basicop := 0;
    current_staircase:= [];
    current_Gb       := [1];
    for mon in list_mon do
        new_Gb       := current_Gb;
        new_staircase:= current_staircase;
        for g in current_Gb do
            if (divide (mon,g,'q')) then
                if (cmp (q,g)) then
                    basicop:= basicop + 6*(1+nops(select (h->cmp (h,g),
                                                          current_staircase)));
                    basicop:= basicop + 2*nops(new_Gb) + 2*nops(new_staircase);
                    if (q = g and g in staircase) then
                        new_Gb:= mergeOrderedLists (remove (h->h=g,new_Gb),
                                                    [seq (g*vars[-i],i=1..nops(vars))],
                                                    mon_ord);
                        new_staircase:= mergeOrderedLists (new_staircase,
                                                           [g],mon_ord);
                    end if:
                end if:
            end if:
        end do:
        current_Gb       := new_Gb;
        current_staircase:= new_staircase;
    end do:
    return basicop;
end proc:

MatrixScalarFGLM:= proc (Tab,vars:=['x','y'],mon_ord:=tdeg(op(vars)),
                         cols_mon:=listMonLessOrEqualThan (vars,mon_ord,vars[1]^10),
                         rows_mon:=cols_mon,
                         lattice_mon:=cols_mon,p:=mychar,
                         full_answer:=false)
description "Compute the minimal relations satisfied by Tab for all the monomials "
    "in cols_mon + rows_mon in characteristic p."
    "If full_answer is set to true, returns a set of triplets [m,r,s] where\n"
    "- m is a monomial;\n"
    "- r is a relation with leading monomial m;\n"
    "- s is the shift of the relation, that is up to which monomial "
    "it has been tested.\n"
    "Otherwise returns a set of r only.";
local H,firstIndRows,usefulStaircase,fullStaircase,firstIndCols,
    subH,next_mon,R,rank,col_next_mon,u,i;
    userinfo (1,ScalarFGLM,`*** Tab `, Tab,
              `in variables`, vars, `in`, nops(cols_mon), `times`, nops(rows_mon),
              `in characteristic`, p, `***`);
    H                := from2SetsToMultiHankelMatrix (Tab,vars,rows_mon,cols_mon);
    firstIndRows,rank:= firstIndependentRows (H,p,min(Dimension(H)));
    if (rank = 0) then
        if (full_answer) then
            return [[1,1,rows_mon]];
        else
            return [1];
        end if:
    elif (rank = min(Dimension(H))) then return [];
    end if:
    if (rows_mon = cols_mon) then
        firstIndCols:= firstIndRows;
    else
        firstIndCols,rank:= firstIndependentColumns (H,p,rank);
    end if:
    firstIndRows     := convert (firstIndRows,list);
    firstIndCols     := convert (firstIndCols,list);
    usefulStaircase:= map (m->cols_mon[m], firstIndCols);
    userinfo (2,ScalarFGLM,`Useful staircase`, usefulStaircase);
    fullStaircase  := sparseStaircaseStabilization (usefulStaircase,cols_mon,
                                                    lattice_mon,mon_ord);
    userinfo (2,ScalarFGLM,`Staircase`, fullStaircase);
    next_mon    := [op({op(cols_mon)} minus {op(fullStaircase)})];
    if (next_mon = []) then return [];
    end if:
    next_mon    := InterReduce (next_mon,mon_ord);
    subH        := H[firstIndRows,firstIndCols];
    col_next_mon:= from2SetsToMultiHankelMatrix (Tab,vars,
                                                 map (i->rows_mon[i],firstIndRows),
                                                 next_mon);
    if (p = 0) then
        R:= LinearSolve (subH,col_next_mon);
    elif (p < 2^32) then
        subH        := LinearAlgebra[Modular][Mod] (p,subH,integer[8]);
        col_next_mon:= LinearAlgebra[Modular][Mod] (p,col_next_mon,integer[8]);
        R           := LinearAlgebra[Modular][LinearSolve] (p,<subH|col_next_mon>,
                                                            nops(next_mon),
                                                            inplace=false);
    else
        subH        := subH mod p;
        col_next_mon:= col_next_mon mod p;
        R           := LinearSolve (subH,col_next_mon) mod p;
    end if:
    R       := [seq (redp (next_mon[i] - DotProduct (Column (R,i),usefulStaircase),p),
                     i=1..nops(next_mon))];
    userinfo (2,ScalarFGLM,`***\t New potential relations`, R, `***`);
    if (full_answer) then
        return [seq ([next_mon[i], R[i], rows_mon],i=1..nops(next_mon))];
    else
        return R;
    end if:
end proc:

(*
PolynomialScalarFGLM:= proc (Tab,vars:=['x','y'],mon_ord:=tdeg(op(vars)),
                         cols_mon:=listMonLessOrEqualThan (vars,mon_ord,vars[1]^10),
                         rows_mon:=cols_mon,
                         p:=mychar,
                         full_answer:=false)
description "Compute the minimal relations satisfied by Tab for all the monomials "
    "in cols_mon + rows_mon in characteristic p."
    "If full_answer is set to true, returns a set of triplets [m,r,s] where\n"
    "- m is a monomial;\n"
    "- r is a relation with leading monomial m;\n"
    "- s is the shift of the relation, that is up to which monomial "
    "it has been tested.\n"
    "Otherwise returns a set of r only.";
local P,M,R,Rp,Rm,m,G,S,cmp;
    userinfo (1,ScalarFGLM,`*** Tab `, Tab,
              `in variables`, vars, `in`, nops(cols_mon), `times`, nops(rows_mon),
              `in characteristic`, p, `***`);
    cmp:= orderComparison (mon_ord);
    P  := from2SetsToMirrorTruncatedGeneratingSeries (Tab,vars,rows_mon,cols_mon);
    M  := lcm(op(rows_mon))*lcm(op(cols_mon));
    B  := map(v->[v^(1+degree(M,v)),0], vars);
    R  := [[P,1]];
    Rp := [];
    G  := [];
    S  := [];
    H  := [];
    while (op(R) <> []) do
        Rm:= R[1]; R:= subsop (1=NULL,R);
        m := LeadingMonomial(Rm,mon_ord);
        if (not cmp (Rm[1],cols_mon[-1]) or
            isPSFGLMValidRelation (Rm,cols_mon,rows_mon,mon_ord)) then
            G:= [op(G),[m,Rm,rows_mon]];
        else
            Rp:= [op(Rp),Rm];
            R := map (r-> nfLeftPart (r,vars,[op(B),Rm],mon_ord,p),R);
            S := [op(S),m];
            Hp:= remove (h->divide (m,h), H);
            Hp:= Groebner:-Basis ([op(Hp),seq (m*v, v in vars)],mon_ord);
            #TODO
        end if:
    end do:
    if (full_answer) then
        return [seq ([next_mon[i], R[i], rows_mon],i=1..nops(next_mon))];
    else
        return R;
    end if:
end proc:
*)

MatrixScalarFGLM_queries:= proc (Tab,vars,mon_ord,
                                 Gb_LM,
                                 staircase,
                                 cols_mon,
                                 rows_mon:=cols_mon)
description "Compute the number of table queries performed by "
    "ScalarFGLMMatrix to recover the minimal relations satisfied by Tab "
    "for all the monomials in cols_mon + rows_mon.";
local set_mon,t,u;
    userinfo (1,ScalarFGLM,`*** Tab `, Tab,
              `in variables`, vars, `in`, nops(cols_mon), `times`, nops(rows_mon),
              `in characteristic`, p, `***`);
    set_mon:= {seq (seq (t*u, u in cols_mon), t in rows_mon)};
    for t in {op(Gb_LM)} minus {op(cols_mon)} do
        set_mon:= set_mon union {seq (t*u, u in rows_mon)};
    end do;
    return nops (set_mon);
end proc:

MatrixScalarFGLM_basicop:= proc (Tab,vars,mon_ord,
                                 Gb_LM, staircase,
                                 cols_mon,
                                 rows_mon:=cols_mon)
description "Compute the number of basic operations performed by "
    "ScalarFGLMMatrix to recover the minimal relations satisfied by Tab "
    "for all the monomials in cols_mon + rows_mon.";
local m,n,d,r,basicop,i;
    userinfo (1,ScalarFGLM,`*** Tab `, Tab,
              `in variables`, vars, `in`, nops(cols_mon), `times`, nops(rows_mon),
              `in characteristic`, p, `***`);
    m:= max (nops (cols_mon), nops (rows_mon));
    n:= min (nops (cols_mon), nops (rows_mon));
    d:= nops (staircase);
    if (d = 0) then return m*n;
    end if:
    basicop:= 3*add ((m-i)*(n-i),i=1..d);
    r      := nops (Gb_LM);
    basicop:= basicop + d*(2*d-1)*r;
    return basicop;
end proc:

QCMatrixScalarFGLM:= proc (Tab,vars:=['x'],diffvars:=[`if`(nops(vars)=1,'t',
                                                           seq (cat('t',i),
                                                                i=1..nops(vars)))],
                           mon_ord:=tdeg(op(vars),op(diffvars)),
                           cols_mon:=listMonLessOrEqualThan ([op(vars),
                                                              op(diffvars)],
                                                             mon_ord,
                                                             vars[1]^10),
                           rows_mon:=listMonLessOrEqualThan (vars,mon_ord,
                                                             vars[1]^20),
                           p:=mychar,
                           full_answer:=false)
description "Compute the minimal relations satisfied by Tab for all the monomials "
    "in cols_mon + rows_mon in characteristic p."
    "If full_answer is set to true, returns a set of triplets [m,r,s] where\n"
    "- m is a monomial;\n"
    "- r is a relation with leading monomial m;\n"
    "- s is the shift of the relation, that is up to which monomial "
    "it has been tested.\n"
    "Otherwise returns a set of r only.";
local H,firstIndRows,usefulStaircase,fullStaircase,firstIndCols,
    subH,next_mon,R,rank,col_next_mon,rank2,S,i;
    userinfo (1,ScalarFGLM,`*** Tab `, Tab,
              `in variables`, vars, `and`, diffvars,
              `in`, nops(cols_mon),`times`,nops(rows_mon),
              `in characteristic`, p, `***`);
    H:= from2QCSetsToMultiHankelMatrix (Tab,vars,diffvars,
                                        rows_mon,cols_mon,p);
    userinfo (1,ScalarFGLM,`*** Matrix built`, `***`);
    #See also LinearAlgebra[Modular][RankProfile]
    firstIndCols,rank:= firstIndependentColumns (H,p,min(Dimension(H)));
    if (rank = 0) then
        if (full_answer) then
            return [[1,1,rows_mon]];
        else
            return [1];
        end if:
    elif (rank = min(Dimension(H))) then return [];
    end if:
    firstIndCols   := convert (firstIndCols,list);
    usefulStaircase:= map (m->cols_mon[m], firstIndCols);
    userinfo (2,ScalarFGLM,`Useful staircase`,
              cat(nops(usefulStaircase),"/",nops(cols_mon)));
    S:= map (i->cols_mon[i], firstIndCols);
    S,rank2      := precursiveCriterion (S,vars,diffvars,
                                         mon_ord);
    userinfo (2,ScalarFGLM,`Criterion -> removed columns:`, rank-rank2);
    fullStaircase:= sparseStaircaseStabilization (usefulStaircase,cols_mon,
                                                  cols_mon,mon_ord);
    userinfo (2,ScalarFGLM,`Staircase`,
              cat(nops(fullStaircase),"/",nops(cols_mon)));
    if (nops(fullStaircase) = min(Dimension(H))) then return [];
    end if:
    H           := from2QCSetsToMultiHankelMatrix (Tab,vars,diffvars,
                                                   rows_mon,S,p);
    firstIndRows,rank:= firstIndependentRows (H,p,rank);
    firstIndRows     := convert (firstIndRows,list);
    next_mon    := [op({op(cols_mon)} minus {op(fullStaircase)})];
    if (next_mon = []) then return [];
    end if:
    next_mon    := InterReduce (next_mon,mon_ord);
    subH        := H[firstIndRows,1..-1];
    col_next_mon:= from2QCSetsToMultiHankelMatrix (Tab,vars,diffvars,
                                                   map (i->rows_mon[i],firstIndRows),
                                                   next_mon,p);
    userinfo (1,ScalarFGLM,`Relations`);
    if (p = 0) then
        R:= LinearSolve (subH,col_next_mon);
    elif (p < 2^32) then
        subH        := LinearAlgebra[Modular][Mod] (p,subH,integer[8]);
        col_next_mon:= LinearAlgebra[Modular][Mod] (p,col_next_mon,integer[8]);
        R           := LinearAlgebra[Modular][LinearSolve] (p,<subH|col_next_mon>,
                                                            nops(next_mon),
                                                            inplace=false);
    else
        subH        := subH mod p;
        col_next_mon:= col_next_mon mod p;
        R           := LinearSolve (subH,col_next_mon) mod p;
    end if:
    R       := [seq (redp (next_mon[i] - DotProduct (Column (R,i),usefulStaircase),p),
                     i=1..nops(next_mon))];
    userinfo (2,ScalarFGLM,`***\t New potential relations`, next_mon, `***`);
    if (full_answer) then
        return [seq ([next_mon[i], R[i], rows_mon],i=1..nops(next_mon))];
    else
        return R;
    end if:
end proc:

orthProj:= proc (f,p,m,n,b,xvar,Tab,myprime:=0)
local g,xalph;
    g:= f;
    for xalph in b do
        g:= redp (g - fromPolToTabRelation (Tab,xvar,g*m[xalph],myprime)*p[xalph],myprime);
    end do:
    return g;
end proc:

nnext:= proc (n,xvar,b,d,c,s)
local nn,a, bb, i;
    nn:= [op({seq (seq (bb*xvar[i],i=1..n), bb in b)})];
    nn:= remove (i->(i in b),nn);
    nn:= select (i->(i in s),nn);
    nn:= remove (i->(i in d),nn);
    a := {op(b)} union {op(d)} union {op(s)};
    nn:= select (i->({op(expand (i*c))} subset (a)),nn);
    return nn;
end proc:

AGbb:= proc (Tab,xvar,a,mon_order,myprime)
local b,m,c,d,nn,s,t,l,p,q,k,n;
local xalph,palph,xgam,xgam_palph,malph,xgam_palphEval,boolgam,supp;
    n:= nops(xvar):
    b:=[];
    c:=[];
    d:= [];
    nn:=[1];
    s:= a;
    t:= a;
    p:= table ([]); m:= table ([]); q:= table ([]); k:= table ([]);
    userinfo (2,AGbb,`set of terms`,nn);
    while (nn <> []) do
        for xalph in [nn[1]] do
            userinfo (1,AGbb,`xalpha`,xalph);
            palph:= orthProj (xalph,p,m,n,b,xvar,Tab,myprime);
            userinfo (2,AGbb,`palpha`,palph);
            p[xalph]:= palph;
            boolgam:= false;
            for xgam in t do
                userinfo (1,AGbb,`\txgamma`,xgam);
                xgam_palph    := xgam*palph;
                xgam_palphEval:= redp (fromPolToTabRelation (Tab,xvar,xgam_palph,myprime),
                                       myprime);
                coeffs (xgam_palph,xvars,'supp'):
                supp:= {supp}:
                if (supp subset {op(a)} and
                    xgam_palphEval <> 0) then
                    boolgam:= true;
                    userinfo (2,AGbb,`\txgamma`,xgam,`!`);
                    break;
                end if:
            end do:
            if (boolgam) then
                if (myprime = 0) then
                    malph:= xgam/xgam_palphEval;
                else
                    malph:= xgam/xgam_palphEval mod myprime;
                end if:
                m[xalph]:= malph;
                userinfo (2,AGbb,`\tmalph`,malph);
                # q[xalph]:= orthProj (malph,q,p,n,b,xvar,Tab,myprime);
                # b:= [op({op(b),xalph})];
                b:= [op(b),xalph];
                userinfo (2,AGbb,`\tnew b`,b);
                s:= remove (i->i=xalph,s);
                userinfo (2,AGbb,`\tnew s`,s);
                # c:= [op({op(c),xgam})];
                c:= [op(c),xgam];
                userinfo (2,AGbb,`\tnew c`,c);
                t:= remove (i->i=xgam,t);
                userinfo (2,AGbb,`\tnew t`,t);
            else
                userinfo (2,AGbb,`\tno xgam !!`);
                k[xalph]:= palph;
                userinfo (1,AGbb,`relation`,xalph,`\n\t\t OK`);
                d:= [op({op(d),xalph})];
                userinfo (2,AGbb,`\tnew d`,d);
                s:= remove (i->i=xalph,s);
                userinfo (2,AGbb,`\tnew s`,s);
            end if:
        end do:
        nn:= nnext (n,xvar,b,d,c,s);
        nn:= sortListMon (nn,mon_order);
        userinfo (2,AGbb,`set of terms`,nn);
        if (xalpha = x^3) then return;
        end if:
    end do:
    return b,c,p,m,q,k,d;
end proc:



###################
# Main procedures #
#    Adaptive     #
###################

AdaptiveBerlekampMasseySakata:= proc (Tab,vars:=['x','y'],mon_ord:=tdeg(op(vars)),
                                      max_mon:=vars[1]^10,d:= 5,p:=mychar,
                                      full_answer:=false)
description "Compute the minimal relations satisfied by Tab for all the monomials "
    "up to max_mon in characteristic p.\n"
    "If full_answer is set to true, returns a set of triplets [m,r,s] where\n"
    "- m is a monomial;\n"
    "- r is a relation with leading monomial m;\n"
    "- s is the shift of the relation, that is up to which monomial "
    "it has been tested.\n"
    "Otherwise returns a set of r only."
    "Parameter d is an upper bound on the degree of the ideal "
    "to help skipping some testings.";
local R,shift,F,fail,list_mon,toTest,discrepancy,cmp,mon,staircase,newStaircase;
local Gb,g,r,i,j,k,fullStaircase,q,f;
    R:= [[1,1,0]]; F:=[]; fullStaircase:= [];
    list_mon:= listMonOfDegreeLessOrEqualThan (vars,mon_ord,
                                               degree(max_mon,vars));
    cmp     := orderComparison (mon_ord);
    userinfo (1,AdaptiveBerlekampMasseySakata,`*** Tab `, Tab,
              `in variables`, vars, `up to`, max_mon,
              `in characteristic`, p, `***`);
    userinfo (2,AdaptiveBerlekampMasseySakata,`***\t New potential relation`, 1, `***`);
    for mon in list_mon do
        userinfo (1,AdaptiveBerlekampMasseySakata,mon);
        newStaircase:= F;
        # toTest     := select (r->divide (mon,r[1],'q'),R);
        # userinfo (2,AdaptiveBerlekampMasseySakata,`To test`,toTest,`***`);
        # userinfo (2,AdaptiveBerlekampMasseySakata,`Failing`,F,`***`);
        # userinfo (2,AdaptiveBerlekampMasseySakata,`S`,fullStaircase,`***`);
        toTest     := select (r->divide (mon,r[1],'q')
                              and (`or`(seq (divide (f[2],q), f in F))
                                   or nops (mergeOrderedLists (fullStaircase,
                                                               staircaseStabilization
                                                               ([r[1],q],vars,mon_ord),
                                                               mon_ord)) <= d),
                              R);
        if (toTest = []) then
            userinfo (2,AdaptiveBerlekampMasseySakata,`***\tNothing to test\t***`);
            next;
        end if:
        shift      := map (r->mon/r[1],toTest);
        discrepancy:= [seq (redp (fromPolToTabRelation (Tab,vars,
                                                        expand (toTest[i][2]*shift[i]),
                                                        p),
                                  p),i=1..nops(toTest))];
        for i from 1 to nops (toTest) do
            if (discrepancy[i] <> 0) then
                newStaircase:= [op(newStaircase),[redp (toTest[i][2]/discrepancy[i],p),
                                                  shift[i]]];
            else
                R[i][3]:= shift[i];
            end if:
        end do:
        newStaircase:= staircaseMaximum (newStaircase,vars,mon_ord);
        Gb          := fromStaircaseToGb (map (s->s[2],newStaircase),vars,mon_ord);
        for i from 1 to nops(Gb) do
            g:= Gb[i];
            Gb[i]:= [g,g,0];
            for r in R do
                if (divide (g,r[1],'quogr')) then break;
                end if:
            end do:
            if (not divide (mon,g,'quomong')) then
                Gb[i][2]:= expand (quogr*r[2]);
                Gb[i][3]:= 0;
                j:= 0;
                while (cmp (list_mon[j+1]*g,mon)) do
                    j:= j+1;
                end do:
                if (j <> 0) then Gb[i][3]:= list_mon[j];
                end if:
            else
                for j from 1 to nops (F) do
                    if (divide (F[j][2],quomong,'quofailquomong')) then
                        Gb[i][2]:= redp (expand (quogr*r[2]
                                                 -fromPolToTabRelation (Tab,vars,
                                                                        r[2]*mon/r[1],
                                                                        p)*
                                                 quofailquomong*F[j][1]),p);
                        Gb[i][3]:= quomong;
                        break;
                    end if:
                end do:
                if (j = nops(F)+1) then
                    Gb[i][2]:= r[2];
                    Gb[i][3]:= quomong;
                end if:
            end if:
        end do:
        i:= 1; j:= 1; k:=1;
        while (i <= nops (toTest) and j <= nops (Gb)) do
            if (cmp (toTest[i][1],Gb[j][1])) then
                userinfo (2,AdaptiveBerlekampMasseySakata,`***\t`,toTest[i][2],`***`);
                userinfo (2,AdaptiveBerlekampMasseySakata,
                          `if`(discrepancy[i] = 0,
                               `***\t SUCCEEDS with a shift`,
                               `***\t FAILS with a shift`),
                          shift[i],`***`);
                if (toTest[i][1] = Gb[j][1]) then
                    if (toTest[i][2] <> Gb[j][2]) then
                        userinfo (2,BerlekampMasseySakata,`***\t\t but tweaked to`,
                                  Gb[j][2],`***`);
                    end if:
                    j:= j+1;
                else
                    userinfo (2,BerlekampMasseySakata,`***\t\t and is discarded ***`);
                end if:
                i:= i+1;
                k:= k+1;
            else
                while (k <= nops (R) and not cmp (Gb[j][1],R[k][1])) do
                    k:= k+1;
                end do:
                if (Gb[j][1] <> R[k][1]) then
                    userinfo (2,AdaptiveBerlekampMasseySakata,
                              `***\t New potential relation in`,Gb[j][1],`\n`,
                              Gb[j][2],`***`);
                end if:
                j:= j+1;
            end if:
        end do:
        while (i <= nops (toTest)) do
            userinfo (2,AdaptiveBerlekampMasseySakata,`***\t`,toTest[i][2],`***`);
            userinfo (2,AdaptiveBerlekampMasseySakata,
                      `if`(discrepancy[i] = 0,
                           `***\t SUCCEEDS with a shift`,
                           `***\t FAILS with a shift`),
                      shift[i],`***`);
            i:= i+1;
        end do:
        while (j <= nops (Gb)) do
            while (k <= nops (R) and not cmp (Gb[j][1],R[k][1])) do
                k:= k+1;
            end do:
            if (k > nops (R) or Gb[j][1] <> R[k][1]) then
                userinfo (2,AdaptiveBerlekampMasseySakata,
                          `***\t New potential relation in`,Gb[j][1],`\n`,
                          Gb[j][2],`***`);
            end if:
            j:= j+1;
        end do:
        fullStaircase:= staircaseStabilization (map (m->m[2],newStaircase),vars,mon_ord);
        R:= Gb;
        F:= newStaircase;
    end do:
    if (full_answer) then
        return R;
    else
        return map (g->g[2],R);
    end if:
end proc:

NaiveAdaptiveMatrixScalarFGLM:= proc (Tab,vars:=['x','y'],mon_ord:=tdeg(op(vars)),d:=10,
                                      p:=mychar,full_answer:=false)
description "Compute the minimal relations satisfied by Tab for all the monomials "
    "in characteristic p."
    "If full_answer is set to true, returns a set of triplets [m,r,s] where\n"
    "- m is a monomial;\n"
    "- r is a relation with leading monomial m;\n"
    "- s is the shift of the relation, that is up to which monomial "
    "it has been tested.\n"
    "Otherwise returns a set of r only.\n"
    "Parameter d is an upper bound on the degree of the ideal "
    "to prevent an infinite loop.";
local H,subH,rank,next_mon,staircase,R,col_next_mon,tau,r,i,g;
    userinfo (1,AdaptiveScalarFGLM,`*** Tab `, Tab,
              `in variables`, vars, `up to degree`, d,
              `in characteristic`, p, `***`);
    H        := Matrix (0,0);
    subH     := H;
    rank     := 0;
    next_mon := [1];
    staircase:= [];
    R        := [];
    while (next_mon <> []) do
        tau    := next_mon[1];
        userinfo (2,AdaptiveScalarFGLM,`***\t New monomial`, tau, `***`);
        col_next_mon:= from2SetsToMultiHankelMatrix (Tab,vars,staircase,[tau]);
        H           := <<subH,Transpose(col_next_mon)>|
                        <col_next_mon,fromMonToTabElem (Tab,vars,tau^2)>>;
        if (((p = 0 or p > 2^32) and Rank (H) > rank)
            or (p <> 0 and LinearAlgebra[Modular][Rank] (p,H) > rank)) then
            userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
            subH     := H;
            staircase:= [op(staircase),tau];
            next_mon := InterReduce ([op(next_mon[2..-1]),
                                      op(remove (l->`or` (seq (divide (l,g),
                                                               g in map (r->r[1],R))),
                                                 expand (tau*vars)))], mon_ord);
            rank:= rank+1;
            if (rank >= d) then
                userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
                col_next_mon:= from2SetsToMultiHankelMatrix (Tab,vars,
                                                             staircase,
                                                             next_mon);
                if (p = 0) then
                    col_next_mon:= LinearSolve (H,col_next_mon);
                elif (p < 2^32) then
                    H           := LinearAlgebra[Modular][Mod] (p,H,integer[8]);
                    col_next_mon:= LinearAlgebra[Modular][Mod] (p,col_next_mon,
                                                                integer[8]);
                    col_next_mon:= LinearAlgebra[Modular][LinearSolve] (p,
                                                                        <H|
                                                                         col_next_mon>,
                                                                        nops(next_mon),
                                                                        inplace=false);
                else
                    H           := H mod p;
                    col_next_mon:= col_next_mon mod p;
                    col_next_mon:= LinearSolve (subH,col_next_mon) mod p;
                end if:
                r:= nops (R);
                R:= [op(R),seq ([next_mon[i],
                                 redp (next_mon[i]
                                       - DotProduct (Column (col_next_mon,i),
                                                     staircase),p),
                                 staircase],
                                i=1..nops(next_mon))];
                map (r->userinfo (2,AdaptiveScalarFGLM,`***\t New potential relations`,
                          r[2],`***`),R[r+1..-1]);
                if (full_answer) then return R;
                else return map (r->r[2],R);
                end if:
            end if:
        else
            userinfo (2,AdaptiveScalarFGLM,`***\t\t NOT in the staircase! ***`);
            if (p = 0) then
                col_next_mon:= LinearSolve (subH,col_next_mon);
            elif (p < 2^32) then
                subH        := LinearAlgebra[Modular][Mod] (p,subH,integer[8]);
                col_next_mon:= LinearAlgebra[Modular][Mod] (p,col_next_mon,
                                                            integer[8]);
                col_next_mon:= LinearAlgebra[Modular][LinearSolve] (p,
                                                                    <subH|col_next_mon>,
                                                                    1,inplace=false);
            else
                subH        := subH mod p;
                col_next_mon:= col_next_mon mod p;
                col_next_mon:= LinearSolve (subH,col_next_mon) mod p;
            end if:
            R:= [op(R),[tau, redp (tau
                                   - DotProduct (col_next_mon,staircase),p),
                        [op(staircase),tau]]];
            # next_mon:= remove (l->divide (l,tau), next_mon);
            next_mon:= subsop(1=NULL,next_mon);
            userinfo (2,AdaptiveScalarFGLM,`***\t New potential relation`, R[-1][2],
                      `***`);
        end if:
    end do:
    if (full_answer) then return R;
    else return map (r->r[2],R);
    end if:
end proc:

AdaptiveMatrixScalarFGLM:= proc (Tab,vars:=['x','y'],mon_ord:=tdeg(op(vars)),d:=10,
                                 p:=mychar,full_answer:=false)
description "Compute the minimal relations satisfied by Tab for all the monomials "
    "in characteristic p."
    "If full_answer is set to true, returns a set of triplets [m,r,s] where\n"
    "- m is a monomial;\n"
    "- r is a relation with leading monomial m;\n"
    "- s is the shift of the relation, that is up to which monomial "
    "it has been tested.\n"
    "Otherwise returns a set of r only.\n"
    "Parameter d is an upper bound on the degree of the ideal "
    "to prevent an infinite loop.";
local L,U,rank,next_mon,staircase,R,tau,col_tau,row_tau,elem_tau2,r,i,g;
local col_next_mon,rtau;
    userinfo (1,AdaptiveScalarFGLM,`*** Tab `, Tab,
              `in variables`, vars, `up to degree`, d,
              `in characteristic`, p, `***`);
    elem_tau2:= fromMonToTabElem (Tab,vars,1);
    userinfo (2,AdaptiveScalarFGLM,`***\t New monomial`, 1, `***`);
    if (elem_tau2 = 0) then
        userinfo (2,AdaptiveScalarFGLM,`***\t New potential relation`, 1,
                  `***`);
        if (full_answer) then return [[1,1,1]];
        else return [1];
        end if:
    end if:
    userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
    L        := Matrix (1,1,1);
    U        := Matrix (1,1,elem_tau2);
    rank     := 1;
    next_mon := sortListMon (vars,mon_ord);
    staircase:= [1];
    if (d = 1) then
        userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
        R:= map (v->[v,redp (v-fromMonToTabElem (Tab,vars,v)/elem_tau2,p),1],next_mon);
        if (full_answer) then return R;
        else return map (r->r[2],R);
        end if:
    end if:
    R        := [];
    while (next_mon <> []) do
        tau    := next_mon[1];
        userinfo (2,AdaptiveScalarFGLM,`***\t New monomial`, tau, `***`);
        col_tau:= from2SetsToMultiHankelMatrix (Tab,vars,staircase,[tau]);
        if (p = 0) then
            col_tau, row_tau:= LinearSolve (L,col_tau),
                               Transpose (LinearSolve (Transpose (U),col_tau));
        elif (p < 2^32) then
            L              := LinearAlgebra[Modular][Mod] (p,L,integer[8]);
            U              := LinearAlgebra[Modular][Mod] (p,U,integer[8]);
            col_tau         := LinearAlgebra[Modular][Mod] (p,col_tau,integer[8]);
            col_tau, row_tau:= LinearAlgebra[Modular][LinearSolve] (p,<L|col_tau>,
                                                                    1,inplace=false),
                               Transpose (LinearAlgebra[Modular][LinearSolve]
                                          (p,<Transpose(U)|col_tau>,
                                           1,inplace=false));
        else
            L               := L mod p;
            U               := U mod p;
            col_tau, row_tau:= col_tau mod p, row_tau mod p;
            col_tau, row_tau:= LinearSolve (L,col_tau) mod p,
                               Transpose (LinearSolve (Transpose (U),col_tau)) mod p;
        end if:
        elem_tau2:= fromMonToTabElem (Tab,vars,tau^2);
        elem_tau2:= redp (-row_tau.col_tau+elem_tau2,p);
        if (type (elem_tau2,Matrix)) then elem_tau2:= elem_tau2[1,1];
        end if:
        # redp (<<L,row_tau>|<Matrix (rank,1),1>>.
        #       <<U,Matrix (1,rank)>|<col_tau,elem_tau2>>,p);
        # redp (from2SetsToMultiHankelMatrix (Tab,vars,[op(staircase),tau]),p);
        # redp (<<L,row_tau>|<Matrix (rank,1),1>>
        #       .<<U,Matrix (1,rank)>|<col_tau,elem_tau2>>
        #       - from2SetsToMultiHankelMatrix (Tab,vars,[op(staircase),tau]),p);
        if (elem_tau2 <> 0) then
            userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
            L        := <<L,row_tau>|<Matrix (rank,1),1>>;
            U        := <<U,Matrix (1,rank)>|<col_tau,elem_tau2>>;
            staircase:= [op(staircase),tau];
            next_mon := InterReduce ([op(next_mon[2..-1]),
                                      op(remove (l->`or` (seq (divide (l,g),
                                                               g in map (r->r[1],R))),
                                                 expand (tau*vars)))], mon_ord);
            rank:= rank+1;
            if (rank >= d) then
                userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
                col_next_mon:= from2SetsToMultiHankelMatrix (Tab,vars,
                                                             staircase,
                                                             next_mon);
                if (p = 0) then
                    col_next_mon:= LinearSolve (U,
                                                LinearSolve (L,col_next_mon));
                elif (p < 2^32) then
                    L:= LinearAlgebra[Modular][Mod] (p,L,integer[8]);
                    U:= LinearAlgebra[Modular][Mod] (p,U,integer[8]);
                    col_next_mon:= LinearAlgebra[Modular][Mod] (p,col_next_mon,
                                                                integer[8]);
                    col_next_mon:= LinearAlgebra[Modular][LinearSolve] (p,
                                                                        <L|
                                                                         col_next_mon>,
                                                                        nops(next_mon),
                                                                        inplace=false);
                    col_next_mon:= LinearAlgebra[Modular][LinearSolve] (p,
                                                                        <U|
                                                                         col_next_mon>,
                                                                        nops(next_mon),
                                                                        inplace=false);
                else
                    L           := L mod p;
                    U           := U mod p;
                    col_next_mon:= col_next_mon mod p;
                    col_next_mon:= LinearSolve (L,col_next_mon) mod p;
                    col_next_mon:= LinearSolve (U,col_next_mon) mod p;
                end if:
                r:= nops (R);
                R:= [op(R),seq ([next_mon[i],
                                 redp (next_mon[i]
                                       - DotProduct (Column (col_next_mon,i),
                                                     staircase),p),
                                 staircase],
                                i=1..nops(next_mon))];
                map (r->userinfo (2,AdaptiveScalarFGLM,`***\t New potential relations`,
                          r[2],`***`),R[r+1..-1]);
                if (full_answer) then return R;
                else return map (r->r[2],R);
                end if:
            end if:
        else
            userinfo (2,AdaptiveScalarFGLM,`***\t\t NOT in the staircase! ***`);
            if (p = 0) then
                col_tau:= LinearSolve (U,col_tau);
            elif (p < 2^32) then
                U      := LinearAlgebra[Modular][Mod] (p,U,integer[8]);
                col_tau:= LinearAlgebra[Modular][Mod] (p,col_tau,integer[8]);
                col_tau:= LinearAlgebra[Modular][LinearSolve] (p,<U|col_tau>,
                                                                    1,inplace=false);
            else
                U      := U mod p;
                col_tau:= col_tau mod p;
                col_tau:= LinearSolve (U,col_tau) mod p;
            end if:
            rtau:= redp (tau - DotProduct (col_tau,staircase),p);
            R   := [op(R),[tau, rtau, [op(staircase),tau]]];
            # next_mon:= remove (l->divide (l,tau), next_mon);
            next_mon:= subsop(1=NULL,next_mon);
            userinfo (2,AdaptiveScalarFGLM,`***\t New potential relation`, rtau,
                      `***`);
        end if:
    end do:
    if (full_answer) then return R;
    else return map (r->r[2],R);
    end if:
end proc:

AdaptiveMatrixScalarFGLM_queries:= proc (Tab,vars,mon_ord,
                                         Gb_LM,S)
description "Compute the number of table queries performed by "
    "AdaptiveScalarFGLM to recover the minimal relation satisfied by Tab.";
local queried_mon,next_mon,staircase,rank,R,col_next_mon,tau,r,s,n,g;
    userinfo (1,AdaptiveScalarFGLM,`*** Tab `, Tab,
              `in variables`, vars, `up to degree`, nops(S));
    userinfo (2,AdaptiveScalarFGLM,`***\t New monomial`, 1, `***`);
    if (Gb_LM = [1] or S = []) then
        userinfo (2,AdaptiveScalarFGLM,`***\t\t NOT in the staircase! ***`);
        return 1;
    end if:
    userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
    if (S = [1]) then
        userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
        return 1+nops(vars);
    end if:
    queried_mon:= [1];
    rank     := 1;
    next_mon := sortListMon (vars,mon_ord);
    staircase:= [1];
    R        := [];
    while (next_mon <> []) do
        tau    := next_mon[1];
        userinfo (2,AdaptiveScalarFGLM,`***\t New monomial`, tau, `***`);
        queried_mon:= mergeOrderedLists (queried_mon,expand ([op(staircase),tau]*tau),
                                         mon_ord);
        if (tau in S) then
            userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
            staircase:= [op(staircase),tau];
            next_mon := InterReduce ([op(next_mon[2..-1]),
                                      op(remove (l->`or` (seq (divide (l,g),
                                                               g in R)),
                                                 expand (tau*vars)))], mon_ord);
            rank:= rank+1;
            if (rank >= nops(S)) then
                userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
                userinfo (2,AdaptiveScalarFGLM,`***\t New monomials`,
                          next_mon,`***`);
                queried_mon:= mergeOrderedLists (queried_mon,
                                                 sortListMon([op({seq
                                                                  (seq (s*n,
                                                                        n in next_mon),
                                                                   s in staircase)})],
                                                             mon_ord), mon_ord);
                return nops(queried_mon);
            end if:
        else
            userinfo (2,AdaptiveScalarFGLM,`***\t\t NOT in the staircase! ***`);
            # next_mon:= remove (l->divide (l,tau), next_mon);
            next_mon:= subsop(1=NULL,next_mon);
            R:= [op(R),tau];
            # userinfo (2,AdaptiveScalarFGLM,
            #           `***\t New potential relation with leading monomial`, R[-1],
            #           `***`);
        end if:
    end do:
    return nops(queried_mon);
end proc:

AdaptiveMatrixScalarFGLM_basicop:= proc (Tab,vars,mon_ord,
                                         Gb_LM,S)
description "Compute the number of basic operations performed by "
    "AdaptiveScalarFGLM to recover the minimal relation satisfied by Tab.";
local basicop,rank,next_mon,staircase,R,tau,g;
    userinfo (1,AdaptiveScalarFGLM,`*** Tab `, Tab,
              `in variables`, vars, `up to degree`, nops(d));
    userinfo (2,AdaptiveScalarFGLM,`***\t New monomial`, 1, `***`);
    if (Gb_LM = [1] or S = []) then
        userinfo (2,AdaptiveScalarFGLM,`***\t\t NOT in the staircase! ***`);
        return 1;
    end if:
    userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
    if (S = [1]) then
        userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
        return 1+2*nops(vars);
    end if:
    basicop  := 1;
    rank     := 1;
    next_mon := sortListMon (vars,mon_ord);
    staircase:= [1];
    R        := [];
    while (next_mon <> []) do
        tau    := next_mon[1];
        userinfo (2,AdaptiveScalarFGLM,`***\t New monomial`, tau, `***`);
        basicop:= basicop + (rank+1)*(rank+4) + 1;
        if (tau in S) then
            userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
            staircase:= [op(staircase),tau];
            next_mon := InterReduce ([op(next_mon[2..-1]),
                                      op(remove (l->`or` (seq (divide (l,g),
                                                               g in map (r->r[1],R))),
                                                 expand (tau*vars)))], mon_ord);
            rank:= rank+1;
            if (rank >= nops(S)) then
                userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
                basicop:= basicop + rank*(rank+1)*nops (next_mon);
                return basicop;
            end if:
        else
            userinfo (2,AdaptiveScalarFGLM,`***\t\t NOT in the staircase! ***`);
            basicop:= basicop + rank*(rank+1)/2;
            # next_mon:= remove (l->divide (l,tau), next_mon);
            next_mon:= subsop(1=NULL,next_mon);
            R:= [op(R),tau];
        end if:
    end do:
    return basicop;
end proc:

NaiveAdaptivePolynomialScalarFGLM:= proc (Tab,vars:=['x','y'],mon_ord:=tdeg(op(vars)),
                                          d:=10,p:=mychar,full_answer:=false)
description "Compute the minimal relations satisfied by Tab for all the monomials "
    "in characteristic p."
    "If full_answer is set to true, returns a set of triplets [m,r,s] where\n"
    "- m is a monomial;\n"
    "- r is a relation with leading monomial m;\n"
    "- s is the shift of the relation, that is up to which monomial "
    "it has been tested.\n"
    "Otherwise returns a set of r only.\n"
    "Parameter d is an upper bound on the degree of the ideal "
    "to prevent an infinite loop.";
local PolS,BS,RS,next_mon,staircase,G,lcmS,lcmS2,TwoS,tau,lcmStau,lcmStau2,extra_mon;
local PolStau,BStau,RStau,rtau,rank,g,s,i;
    userinfo (1,AdaptiveScalarFGLM,`*** Tab `, Tab,
              `in variables`, vars, `up to degree`, d,
              `in characteristic`, p, `***`);
    extra_mon:= fromMonToTabElem (Tab,vars,1);
    userinfo (2,AdaptiveScalarFGLM,`***\t New monomial`, 1, `***`);
    if (extra_mon = 0) then
        userinfo (2,AdaptiveScalarFGLM,`***\t New potential relation`, 1,
                  `***`);
        if (full_answer) then return [[1,1,1]];
        else return [1];
        end if:
    end if:
    userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
    BS       := map (v->[v,0], vars);
    RS       := [[extra_mon,1]];
    next_mon := sortListMon (vars,mon_ord);
    staircase:= [1];
    rank     := 1;
    PolS     := extra_mon;
    if (d = 1) then
        userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
        G:= map (v->[v,redp (v-fromMonToTabElem (Tab,vars,v)/elem_tau2,p),1],next_mon);
        if (full_answer) then return G;
        else return map (g->g[2],G);
        end if:
    end if:
    G    := [];
    lcmS := 1;
    lcmS2:= 1;
    TwoS := [1];
    while (next_mon <> []) do
        tau     := next_mon[1];
        userinfo (2,AdaptiveScalarFGLM,`***\t New monomial`, tau, `***`);
        lcmStau  := lcm (tau,lcmS);
        lcmStau2 := lcmStau^2;
        extra_mon:= remove (m->m in TwoS, map(n->n*tau,[op(staircase),tau]));
        PolStau  := redp (mirrorPolynomial (from1SetToTruncatedGeneratingSeries
                                            (Tab,vars,extra_mon),
                                            vars,lcmStau2),p);
        BStau    := map (v->[v^(degree (lcmStau2,v)+1),0], vars);
        RStau    := map (r->nfLeftPart([expand (lcmStau2/lcmS2 * r[1] +
                                                redp(PolStau * r[2],p)),r[2]],vars,
                                       BStau,mon_ord,p),RS);
        PolStau  := RStau[1][1];
        rtau     := nfLeftPart ([expand (PolStau*tau),tau],vars,BStau,mon_ord,p);
        rtau     := nfLeftHigherPart (rtau,vars,RStau,lcmStau2,
                                      map (g->g[1],G),mon_ord,p);
        if (LeadingMonomial (leftHigherPart (rtau[1],vars,lcmStau2,map(g->g[1],G),
                                             mon_ord,p),mon_ord) = lcmStau2/tau) then
            userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
            userinfo (2,AdaptiveScalarFGLM,`***\t\t Leading Monomials higher part`,
                      LeadingMonomial (leftHigherPart (rtau[1],vars,lcmStau2,
                                                       map(g->g[1],G), mon_ord,p),mon_ord),
                      `and truly`,
                      LeadingMonomial (rtau[1],mon_ord), `***`);
            staircase:= [op(staircase),tau];
            next_mon := InterReduce ([op(next_mon[2..-1]),
                                      op(remove (l->`or` (seq (divide (l,g),
                                                               g in map (r->r[1],G))),
                                                 expand (tau*vars)))], mon_ord);
            rank:= rank+1;
            RS       := [op(RStau),rtau];
            BS       := BStau;
            TwoS     := mergeOrderedLists (TwoS,extra_mon,mon_ord);
            lcmS     := lcmStau;
            lcmS2    := lcmStau2;
            PolS     := PolStau;
            if (rank >= d) then
                userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
                lcmStau  := lcm (op(next_mon),lcmS);
                lcmStau2 := lcmS*lcmStau;
                extra_mon:= sortListMon ([op(remove (m->m in TwoS,
                                                     {seq (seq (s*tau,s in staircase),
                                                           tau in next_mon)}))],mon_ord);
                PolStau  := redp (mirrorPolynomial
                                  (from1SetToTruncatedGeneratingSeries
                                   (Tab,vars,extra_mon),vars,lcmStau2),p);
                BStau    := map (v->[v^(degree (lcmStau2,v)+1),0], vars);
                RStau    := map (r->nfLeftPart([expand (lcmStau/lcmS * r[1] +
                                                        redp(PolStau * r[2],p)),
                                                r[2]],vars,BStau,mon_ord,p),RS);
                PolStau  := RStau[1][1];
                next_mon := map (n->nfLeftPart ([expand (PolStau*n),n],
                                                vars,BStau,mon_ord,p),next_mon);
                next_mon := map (n->[n[2],op(nfLeftHigherPart (n,vars,RStau,lcmStau2,
                                                               map (g->g[1],G),
                                                               mon_ord,p))],next_mon);
                for g in next_mon do
                    userinfo (2,AdaptiveScalarFGLM,`***\t New potential relation`,
                              g[3],`***`);
                end do:
                G:= [op(G),seq ([next_mon[i][1],next_mon[i][3],staircase],
                                i=1..nops(next_mon))];
                if (full_answer) then return G;
                else return map (g->g[2],G);
                end if:
            end if:
        else
            userinfo (2,AdaptiveScalarFGLM,`***\t\t NOT in the staircase! ***`);
            userinfo (2,AdaptiveScalarFGLM,`***\t\t Leading Monomials higher part`,
                      LeadingMonomial (leftHigherPart (rtau[1],vars,lcmStau2,
                                                       map(g->g[1],G), mon_ord,p),mon_ord),
                      `and truly`,
                      LeadingMonomial (rtau[1],mon_ord), `***`);
            G:= [op(G),[tau,rtau[2],[op(staircase),tau]]];
            # next_mon:= remove (l->divide (l,tau), next_mon);
            next_mon:= subsop(1=NULL,next_mon);
            userinfo (2,AdaptiveScalarFGLM,`***\t New potential relation`, rtau[2],
                      `***`);
        end if:
    end do:
    if (full_answer) then return G;
    else return map (g->g[2],G);
    end if:
end proc:

AdaptivePolynomialScalarFGLM:= proc (Tab,vars:=['x','y'],mon_ord:=tdeg(op(vars)),
                                     d:=10,p:=mychar,full_answer:=false)
description "Compute the minimal relations satisfied by Tab for all the monomials "
    "in characteristic p."
    "If full_answer is set to true, returns a set of triplets [m,r,s] where\n"
    "- m is a monomial;\n"
    "- r is a relation with leading monomial m;\n"
    "- s is the shift of the relation, that is up to which monomial "
    "it has been tested.\n"
    "Otherwise returns a set of r only.\n"
    "Parameter d is an upper bound on the degree of the ideal "
    "to prevent an infinite loop.";
local PolS,BS,RS,next_mon,staircase,GS,GStau,lcmS,lcmS2,TwoS,tau,lcmStau,lcmStau2;
local extra_mon,PolStau,BStau,RStau,rtau,rank,g,i,initiated,T,n,lm,s;
    userinfo (1,AdaptiveScalarFGLM,`*** Tab `, Tab,
              `in variables`, vars, `up to degree`, d,
              `in characteristic`, p, `***`);
    extra_mon:= fromMonToTabElem (Tab,vars,1);
    userinfo (2,AdaptiveScalarFGLM,`***\t New monomial`, 1, `***`);
    if (PolS = 0) then
        userinfo (2,AdaptiveScalarFGLM,`***\t New potential relation`, 1,
                  `***`);
        if (full_answer) then return [[1,1,1]];
        else return [1];
        end if:
    end if:
    userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
    BS       := map (v->[v,0], vars);
    RS       := [[extra_mon,1]];
    next_mon := sortListMon (vars,mon_ord);
    staircase:= [1];
    T        := table ([1=1]);
    rank     := 1;
    PolS     := extra_mon;
    if (d = 1) then
        userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
        GS:= map (v->v=[v,redp (v-fromMonToTabElem (Tab,vars,v)/elem_tau2,p),1],
                  next_mon);
        if (full_answer) then return GS;
        else return map (g->g[2],GS);
        end if:
    end if:
    GS   := [];
    lcmS := 1;
    lcmS2:= 1;
    TwoS := [1];
    while (next_mon <> []) do
        tau      := next_mon[1];
        userinfo (2,AdaptiveScalarFGLM,`***\t New monomial`, tau, `***`);
        lcmStau  := lcm (tau,lcmS);
        lcmStau2 := lcmStau^2;
        extra_mon:= remove (m->m in TwoS, map(n->n*tau,[op(staircase),tau]));
        PolStau  := redp (mirrorPolynomial (from1SetToTruncatedGeneratingSeries
                                            (Tab,vars,extra_mon),
                                            vars,lcmStau2),p);
        BStau    := map (v->[v^(degree (lcmStau2,v)+1),0], vars);
        RStau    := map (r->nfLeftPart([expand (lcmStau2/lcmS2 * r[1] +
                                                redp(PolStau * r[2],p)),r[2]],vars,BStau,
                                       mon_ord,p),RS);
        GStau    := map (g->[g[1],
                             op(nfLeftPart([expand (lcmStau2/lcmS2 * g[2] +
                                                    redp(PolStau * g[3],p)),g[3]],
                                           vars,BStau,mon_ord,p)),
                             g[4]],GS);
        PolStau  := RStau[1][1];
        initiated:= false;
        for i from 1 to nops (vars) do
            if (degree (tau,vars[i]) >= 2) then
                rtau:= nfLeftHigherPart (RStau[T[tau/vars[i]^2]],vars,
                                         [RStau[T[tau/vars[i]]],op(BStau)],
                                         lcmStau2,map (g->g[1],GStau),
                                         mon_ord,p);
                rtau:= nfLeftPart(rtau,vars,BStau,mon_ord,p);
                initiated:= true;
                lm       := tau/vars[i]^2;
                break;
            end if:
        end do:
        if (not initiated) then
            for i from 1 to nops (vars) do
                if (degree (tau,vars[i]) = 1) then
                    rtau:= nfLeftPart (vars[i]*RStau[T[tau/vars[i]]],vars,
                                       BStau,mon_ord,p);
                    initiated:= true;
                    lm       := tau/vars[i];
                    break;
                end if:
            end do:
        end if:
        if (not initiated) then error "BUG";
        end if:
        rtau:= nfRightPart (rtau,vars,GStau,mon_ord,p);
        rtau:= nfLeftPart (rtau,vars,BStau,mon_ord,p);
        rtau:= nfLeftHigherPart (rtau,vars,RStau,lcmStau2,
                                 map (g->g[1],GStau),mon_ord,p);
        if (LeadingMonomial (leftHigherPart (rtau[1],vars,lcmStau2,map(g->g[1],GStau),
                                             mon_ord,p),mon_ord) = lcmStau2/tau) then
            userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
            userinfo (2,AdaptiveScalarFGLM,`***\t\t Leading Monomials higher part`,
                      LeadingMonomial (leftHigherPart (rtau[1],vars,lcmStau2,
                                                       map(g->g[1],GStau), mon_ord,p),mon_ord),
                      `and truly`,
                      LeadingMonomial (rtau[1],mon_ord), `***`);
            staircase:= [op(staircase),tau];
            next_mon := InterReduce ([op(next_mon[2..-1]),
                                      op(remove (l->`or` (seq (divide (l,g),
                                                               g in map (r->r[1],GS))),
                                                 expand (tau*vars)))], mon_ord);
            rank     := rank+1;
            RS       := [op(RStau),rtau];
            BS       := BStau;
            GS       := GStau;
            TwoS     := mergeOrderedLists (TwoS,extra_mon,mon_ord);
            lcmS     := lcmStau;
            lcmS2    := lcmStau2;
            PolS     := PolStau;
            T[tau]   := rank;
            if (rank >= d) then
                userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
                userinfo (2,AdaptiveScalarFGLM,`***\t rank ***`,rank,`* d * `, d);
                lcmStau  := lcm (op(next_mon),lcmS);
                lcmStau2 := lcmS*lcmStau;
                extra_mon:= sortListMon ([op(remove (m->m in TwoS,
                                                     {seq (seq (s*tau, s in staircase),
                                                           tau in next_mon)}))],mon_ord);
                PolStau := redp (mirrorPolynomial
                                 (from1SetToTruncatedGeneratingSeries
                                  (Tab,vars,extra_mon),vars,lcmStau2),p);
                BStau   := map (v->[v^(degree (lcmStau2,v)+1),0], vars);
                RStau   := map (r->nfLeftPart([expand (lcmStau/lcmS * r[1] +
                                                       redp(PolStau * r[2],p)),
                                               r[2]],vars,BStau,mon_ord,p),RS);
                GStau    := map (g->[g[1],
                                     op(nfLeftPart([expand (lcmStau/lcmS * g[2] +
                                                            redp(PolStau * g[3],p)),
                                                    g[3]],vars,BStau,mon_ord,p)),
                                     g[4]],GS);
                PolStau:= RStau[1][1];
                for n from 1 to nops(next_mon) do
                    tau      := next_mon[n];
                    initiated:= false;
                    for i from 1 to nops (vars) do
                        if (degree (tau,vars[i]) >= 2) then
                            rtau:= nfLeftPart (RStau[T[tau/vars[i]^2]],vars,
                                               [RStau[T[tau/vars[i]]],op(BStau)],
                                               mon_ord,p);
                            initiated:= true;
                            lm       := tau/vars[i]^2;
                            break;
                        end if:
                    end do:
                    if (not initiated) then
                        for i from 1 to nops (vars) do
                            if (degree (tau,vars[i]) = 1) then
                                rtau:= nfLeftPart (vars[i]*RStau[T[tau/vars[i]]],vars,
                                                   BStau,mon_ord,p);
                                initiated:= true;
                                lm       := tau/vars[i];
                                break;
                            end if:
                        end do:
                    end if:
                    rtau       := nfRightPart (rtau,vars,GStau,mon_ord,p);
                    rtau       := nfLeftPart (rtau,vars,BStau,mon_ord,p);
                    # rtau       := nfLeftHigherPart (rtau,vars,RStau[T[lm]+1..-1],
                    #                                 lcmStau2,
                    #                                 map (g->g[1],GStau),mon_ord,p);
                    rtau       := nfLeftHigherPart (rtau,vars,RStau,lcmStau2,
                                                    map (g->g[1],GStau),mon_ord,p);
                    userinfo (2,AdaptiveScalarFGLM,`***\t New potential relation`,
                              rtau[2],`***`);
                    next_mon[n]:= [tau,rtau[2],staircase];
                end do:
                GStau:= [op(map (g->subsop(2=NULL,g),GStau)),
                         op(next_mon)];
                if (full_answer) then return GStau;
                else return map (g->g[2],GStau);
                end if:
            end if:
        else
            userinfo (2,AdaptiveScalarFGLM,`***\t\t NOT in the staircase! ***`);
            userinfo (2,AdaptiveScalarFGLM,`***\t\t Leading Monomials higher part`,
                      LeadingMonomial (leftHigherPart (rtau[1],vars,lcmStau2,
                                                       map(g->g[1],GS), mon_ord,p),mon_ord),
                      `and truly`,
                      LeadingMonomial (rtau[1],mon_ord), `***`);
            # GS:= [op(GS),[tau,rtau[1],rtau[2],[op(staircase),tau]]];
            GS:= [op(GS),[tau,Groebner:-NormalForm(rtau[2]*
                                                   from1SetToMirrorTruncatedGeneratingSeries
                                                   (Tab,vars,TwoS),map(g->g[1],BS),mon_ord,
                                                   characteristic=p),
                          rtau[2],[op(staircase),tau]]];
            # next_mon:= remove (l->divide (l,tau), next_mon);
            next_mon:= subsop(1=NULL,next_mon);
            userinfo (2,AdaptiveScalarFGLM,`***\t New potential relation`, rtau[2],
                      `***`);
        end if:
    end do:
    if (full_answer) then return map (g->[g[1],g[3],g[4]],GS);
    else return map (g->g[3],GS);
    end if:
end proc:

ScalarNumerator:= proc (Tab,vars:=['x'],d:=10,P:=0,p:=mychar)
description "Returns the Scalar Numerator Omega (Tab,P), "
    "See Hyun et al. 2020, where m=1.";
    local i:
    return quo (redp (expand (P*add (Tab(d-1-i)*vars^i,i=0..d-1)),p),vars^d,vars);
end proc:

BlockParametrization:= proc (Tab,vars:=['x','y'],lf:=vars[-1],
                             extravar:=vars[-1],d:=10,p:=mychar)
description "Computes the parametrizations of the radical of the ideal "
    "of relations of Tab, assuming it is in shape position for "
    "LEX (vars[1]>...>vars[-1]>T) where T is a new variable satisfying "
    "T=lf, a linear form in vars[1],...,vars[-1]. See Hyun et al. 2020.";
local P,Q,n,C,k,lfpow,lfTab,st,i;
    st:= time ();
    n:= nops (vars);
    lfpow:= [seq (redp (expand (lf^i),p),i=0..2*d-1)];
    lfTab:= proc (i)
    option remember;
        return fromPolToTabRelation (Tab,vars,lfpow[i+1],p);
    end proc:
    SeqVector (lfTab,2*d);
    # printf("Sequence computed\n"):
    P:= PolynomialBerlekampMassey (lfTab,extravar,extravar^(2*d-1),p,false);
    # printf("Elimination polynomial computed\n"):
    if (p = 0) then Q:= sqrfree (P);
    else Q:= Sqrfree (P) mod p;
    end if:
    Q:= redp (expand (Q[1]*mul (Q[2][i][1],i=1..nops(Q[2]))),p);
    C:= table ([]);
    C[1]:= ScalarNumerator (lfTab,extravar,d,P,p);
    # print ("C1=",C[1]);
    if (p = 0) then gcdex (C[1],Q,extravar,'invC1');
    else Gcdex (C[1],Q,extravar,'invC1') mod p;
    end if:
    # print ("invC1=",invC1);
    for k from 1 to n do
        lfTab:= proc (i)
        option remember;
            return fromPolToTabRelation (Tab,vars,lfpow[i+1]*vars[k],p);
        end proc:
        C[vars[k]]:= redp (expand (ScalarNumerator
                                  (lfTab,extravar,d,P,p)*invC1),p);
        if (p = 0) then C[vars[k]]:= rem (C[vars[k]],Q,extravar);
        else C[vars[k]]:= Rem (C[vars[k]],Q,extravar) mod p;
        end if:
    end do:
    # printf ("elapsed time: %f\n",time()-st);
    return [Q,seq (vars[n-k]-C[vars[n-k]] mod p,k=0..n-1),P];
end proc:

NaiveQCAdaptiveMatrixScalarFGLM:= proc (Tab,vars:=['x'],diffvars:=['t'],
                                        mon_ord:=tdeg(op(vars),op(diffvars)),
                                        d:=10, p:=mychar,full_answer:=false)
description "Compute the minimal relations with polynomial coefficients"
    "satisfied by Tab for all the monomials "
    "in characteristic p."
    "If full_answer is set to true, returns a set of triplets [m,r,s] where\n"
    "- m is a monomial;\n"
    "- r is a relation with leading monomial m;\n"
    "- s is the shift of the relation, that is up to which monomial "
    "it has been tested.\n"
    "Otherwise returns a set of r only.\n"
    "Parameter d is an upper bound on the degree of the ideal "
    "to prevent an infinite loop.";
local H,subH,rank,next_mon,next_shift,staircase,shifts,R;
local row_next_shift,col_next_mon,tau,upsilon,r,i,g;
    userinfo (1,AdaptiveScalarFGLM,`*** Tab `, Tab,
              `in variables`, diffvars, `and`, vars, `up to degree`, d,
              `in characteristic`, p, `***`);
    H            := Matrix (0,0);
    subH         := H;
    rank         := 0;
    next_mon     := [1];
    next_shift   := [1];
    staircase    := [];
    shifts       := [];
    R            := [];
    while ((not indets (next_mon) subset {op(diffvars)} and next_shift <> [])
           or next_mon = [1]) do
        tau    := next_mon[1];
        upsilon:= next_shift[1];
        userinfo (2,AdaptiveScalarFGLM,`***\t New monomials`, tau,
                  `and`,upsilon,`***`);
        col_next_mon  := from2QCSetsToMultiHankelMatrix (Tab,vars,diffvars,
                                                         shifts,[tau],p);
        row_next_shift:= from2QCSetsToMultiHankelMatrix (Tab,vars,diffvars,
                                                         [upsilon],staircase,p);
        H           := <<subH,row_next_shift>|
                        <col_next_mon,fromQCMonToTabElem (Tab,vars,diffvars,
                                                          tau*upsilon,p)>>;
        if (((p = 0 or p > 2^32) and Rank (H) > rank)
            or (p <> 0 and LinearAlgebra[Modular][Rank] (p,H) > rank)) then
            userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
            subH     := H;
            staircase:= [op(staircase),tau];
            shifts   := [op(shifts),upsilon];
            next_mon := InterReduce ([op(next_mon[2..-1]),
                                      op(remove (l->`or` (seq (divide (l,g),
                                                               g in map (r->r[1],R))),
                                                 [op(expand (tau*vars)),
                                                  op(expand (tau*diffvars))]))],
                                     mon_ord);
            next_shift:= InterReduce ([op(next_shift[2..-1]),
                                       op(remove (l->`or` (seq (divide (l,g),
                                                                g in map (r->r[1],R))),
                                                  expand (upsilon*vars)))],
                                      mon_ord);
            rank:= rank+1;
            if (rank >= d) then
                userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
                col_next_mon:= from2QCSetsToMultiHankelMatrix (Tab,vars,diffvars,
                                                               shifts,
                                                               next_mon,p);
                if (p = 0) then
                    col_next_mon:= LinearSolve (H,col_next_mon);
                elif (p < 2^32) then
                    H           := LinearAlgebra[Modular][Mod] (p,H,integer[8]);
                    col_next_mon:= LinearAlgebra[Modular][Mod] (p,col_next_mon,
                                                                integer[8]);
                    col_next_mon:= LinearAlgebra[Modular][LinearSolve] (p,
                                                                        <H|
                                                                         col_next_mon>,
                                                                        nops(next_mon),
                                                                        inplace=false);
                else
                    H           := H mod p;
                    col_next_mon:= col_next_mon mod p;
                    col_next_mon:= LinearSolve (subH,col_next_mon) mod p;
                end if:
                r:= nops (R);
                R:= [op(R),seq ([next_mon[i],
                                 redp (next_mon[i]
                                       - DotProduct (Column (col_next_mon,i),
                                                     staircase),p),
                                 shifts],
                                i=1..nops(next_mon))];
                map (r->userinfo (2,AdaptiveScalarFGLM,`***\t New potential relations`,
                          r[2],`***`),R[r+1..-1]);
                if (full_answer) then return R;
                else return map (r->r[2],R);
                end if:
            end if:
        else
            userinfo (2,AdaptiveScalarFGLM,`***\t\t NOT in the staircase! ***`);
            if (p = 0) then
                col_next_mon:= LinearSolve (subH,col_next_mon);
            elif (p < 2^32) then
                subH        := LinearAlgebra[Modular][Mod] (p,subH,integer[8]);
                col_next_mon:= LinearAlgebra[Modular][Mod] (p,col_next_mon,
                                                            integer[8]);
                col_next_mon:= LinearAlgebra[Modular][LinearSolve] (p,
                                                                    <subH|col_next_mon>,
                                                                    1,inplace=false);
            else
                subH        := subH mod p;
                col_next_mon:= col_next_mon mod p;
                col_next_mon:= LinearSolve (subH,col_next_mon) mod p;
            end if:
            R:= [op(R),[tau, redp (tau
                                   - DotProduct (col_next_mon,staircase),p),
                        [op(shifts),tau]]];
            next_mon  := subsop(1=NULL,next_mon);
            next_shift:= subsop(1=NULL,next_shift);
            userinfo (2,AdaptiveScalarFGLM,`***\t New potential relation`, R[-1][2],
                      `***`);
        end if:
    end do:
    if (next_shift = []) then
        userinfo (1,AdaptiveScalarFGLM,`***\t Not enough rows`,`***`);
    end if:
    if (full_answer) then return R;
    else return map (r->r[2],R);
    end if:
end proc:

QCAdaptiveMatrixScalarFGLM:= proc (Tab,vars:=['x'],diffvars:=['t'],
                                   mon_ord:=tdeg(op(vars),op(diffvars)),
                                   d:=10, p:=mychar,full_answer:=false)
description "Compute the minimal relations with polynomial coefficients"
    "satisfied by Tab for all the monomials "
    "in characteristic p."
    "If full_answer is set to true, returns a set of triplets [m,r,s] where\n"
    "- m is a monomial;\n"
    "- r is a relation with leading monomial m;\n"
    "- s is the shift of the relation, that is up to which monomial "
    "it has been tested.\n"
    "Otherwise returns a set of r only.\n"
    "Parameter d is an upper bound on the degree of the ideal "
    "to prevent an infinite loop.";
local L,U,rank,next_mon,next_shift,staircase,shifts,R,tau,upsilon;
local col_tau,row_upsilon,elem_tauupsilon,r,col_next_mon,rtau,i,g;
    userinfo (1,AdaptiveScalarFGLM,`*** Tab `, Tab,
              `in variables`, diffvars, `and`, vars, `up to degree`, d,
              `in characteristic`, p, `***`);
    elem_tauupsilon:= fromMonToTabElem (Tab,vars,1);
    userinfo (2,AdaptiveScalarFGLM,`***\t New monomial`, 1, `***`);
    if (elem_tauupsilon = 0) then
        userinfo (2,AdaptiveScalarFGLM,`***\t New potential relation`, 1,
                  `***`);
        if (full_answer) then return [[1,1,1]];
        else return [1];
        end if:
    end if:
    userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
    L         := Matrix (1,1,1);
    U         := Matrix (1,1,elem_tauupsilon);
    rank      := 1;
    next_mon  := sortListMon ([op(vars),op(diffvars)],mon_ord);
    next_shift:= sortListMon (vars,mon_ord);
    staircase := [1];
    shifts    := [1];
    if (d = 1) then
        userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
        R:= map (v->[v,redp (v-fromMonToTabElem (Tab,vars,v)/elem_tauupsilon,p),
                     1],next_mon);
        if (full_answer) then return R;
        else return map (r->r[2],R);
        end if:
    end if:
    R        := [];
    while (not indets (next_mon) subset {op(diffvars)} and next_shift <> []) do
        tau    := next_mon[1];
        upsilon:= next_shift[1];
        userinfo (2,AdaptiveScalarFGLM,`***\t New monomials`, tau,
                  `and`, upsilon,`***`);
        col_tau    := from2QCSetsToMultiHankelMatrix (Tab,vars,diffvars,
                                                      shifts,[tau],p);
        row_upsilon:= from2QCSetsToMultiHankelMatrix (Tab,vars,diffvars,
                                                      staircase,[upsilon],p);
        if (p = 0) then
            col_tau    := LinearSolve (L,col_tau);
            row_upsilon:= Transpose (LinearSolve (Transpose (U),row_upsilon));
        elif (p < 2^32) then
            L          := LinearAlgebra[Modular][Mod] (p,L,integer[8]);
            U          := LinearAlgebra[Modular][Mod] (p,U,integer[8]);
            col_tau    := LinearAlgebra[Modular][Mod] (p,col_tau,integer[8]);
            row_upsilon:= LinearAlgebra[Modular][Mod] (p,row_upsilon,integer[8]);
            col_tau    := LinearAlgebra[Modular][LinearSolve] (p,<L|col_tau>,
                                                               1,inplace=false);
            row_upsilon:= Transpose (LinearAlgebra[Modular][LinearSolve]
                                     (p,<Transpose(U)|row_upsilon>,
                                      1,inplace=false));
        else
            L          := L mod p;
            U          := U mod p;
            col_tau    := LinearSolve (L,col_tau mod p) mod p;
            row_upsilon:= Transpose (LinearSolve (Transpose (U),row_upsilon mod p)
                                     mod p);
        end if:
        elem_tauupsilon:= fromQCMonToTabElem (Tab,vars,diffvars,tau*upsilon);
        elem_tauupsilon:= redp (-row_upsilon.col_tau+elem_tauupsilon,p);
        if (type (elem_tauupsilon,Matrix)) then elem_tauupsilon:= elem_tauupsilon[1,1];
        end if:
        if (elem_tauupsilon <> 0) then
            userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
            L        := <<L,row_upsilon>|<Matrix (rank,1),1>>;
            U        := <<U,Matrix (1,rank)>|<col_tau,elem_tauupsilon>>;
            staircase:= [op(staircase),tau];
            shifts   := [op(shifts),upsilon];
            next_mon := InterReduce ([op(next_mon[2..-1]),
                                      op(remove (l->`or` (seq (divide (l,g),
                                                               g in map (r->r[1],R))),
                                                 [op(expand (tau*vars)),
                                                  op(expand (tau*diffvars))]))],
                                     mon_ord);
            next_shift:= InterReduce ([op(next_shift[2..-1]),
                                       op(remove (l->`or` (seq (divide (l,g),
                                                                g in map (r->r[1],R))),
                                                  expand (upsilon*vars)))],
                                      mon_ord);
            rank:= rank+1;
            if (rank >= d) then
                userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
                col_next_mon:= from2QCSetsToMultiHankelMatrix (Tab,vars,diffvars,
                                                               shifts,
                                                               next_mon,p);
                if (p = 0) then
                    col_next_mon:= LinearSolve (U,
                                                LinearSolve (L,col_next_mon));
                elif (p < 2^32) then
                    L:= LinearAlgebra[Modular][Mod] (p,L,integer[8]);
                    U:= LinearAlgebra[Modular][Mod] (p,U,integer[8]);
                    col_next_mon:= LinearAlgebra[Modular][Mod] (p,col_next_mon,
                                                                integer[8]);
                    col_next_mon:= LinearAlgebra[Modular][LinearSolve] (p,
                                                                        <L|
                                                                         col_next_mon>,
                                                                        nops(next_mon),
                                                                        inplace=false);
                    col_next_mon:= LinearAlgebra[Modular][LinearSolve] (p,
                                                                        <U|
                                                                         col_next_mon>,
                                                                        nops(next_mon),
                                                                        inplace=false);
                else
                    L           := L mod p;
                    U           := U mod p;
                    col_next_mon:= col_next_mon mod p;
                    col_next_mon:= LinearSolve (L,col_next_mon) mod p;
                    col_next_mon:= LinearSolve (U,col_next_mon) mod p;
                end if:
                r:= nops (R);
                R:= [op(R),seq ([next_mon[i],
                                 redp (next_mon[i]
                                       - DotProduct (Column (col_next_mon,i),
                                                     staircase),p),
                                 shifts],
                                i=1..nops(next_mon))];
                map (r->userinfo (2,AdaptiveScalarFGLM,`***\t New potential relations`,
                          r[2],`***`),R[r+1..-1]);
                if (full_answer) then return R;
                else return map (r->r[2],R);
                end if:
            end if:
        else
            userinfo (2,AdaptiveScalarFGLM,`***\t\t NOT in the staircase! ***`);
            if (p = 0) then
                col_tau:= LinearSolve (U,col_tau);
            elif (p < 2^32) then
                U      := LinearAlgebra[Modular][Mod] (p,U,integer[8]);
                col_tau:= LinearAlgebra[Modular][Mod] (p,col_tau,integer[8]);
                col_tau:= LinearAlgebra[Modular][LinearSolve] (p,<U|col_tau>,
                                                               1,inplace=false);
            else
                U      := U mod p;
                col_tau:= col_tau mod p;
                col_tau:= LinearSolve (U,col_tau) mod p;
            end if:
            rtau:= redp (tau - DotProduct (col_tau,staircase),p);
            R   := [op(R),[tau, rtau, [op(shifts),upsilon]]];
            # next_mon:= remove (l->divide (l,tau), next_mon);
            next_mon:= subsop(1=NULL,next_mon);
            next_shift:= subsop(1=NULL,next_shift);
            userinfo (2,AdaptiveScalarFGLM,`***\t New potential relation`, rtau,
                      `***`);
        end if:
    end do:
    if (next_shift = []) then
        userinfo (1,AdaptiveScalarFGLM,`***\t Not enough rows`,`***`);
    end if:
    if (full_answer) then return R;
    else return map (r->r[2],R);
    end if:
end proc:


QCAdaptiveMatrixScalarFGLM:= proc (Tab,vars:=['x'],diffvars:=['t'],
                                   mon_ord:=tdeg(op(vars),op(diffvars)),
                                   d:=10, p:=mychar,full_answer:=false,
                                   cone_gen:=vars,
                                   to_exclude:=[])
description "Compute the minimal relations with polynomial coefficients"
    "satisfied by Tab for all the monomials "
    "in characteristic p."
    "If full_answer is set to true, returns a set of triplets [m,r,s] where\n"
    "- m is a monomial;\n"
    "- r is a relation with leading monomial m;\n"
    "- s is the shift of the relation, that is up to which monomial "
    "it has been tested.\n"
    "Otherwise returns a set of r only.\n"
    "Parameter d is an upper bound on the degree of the ideal "
    "to prevent an infinite loop.";
local L,U,rank,next_mon,next_shift,staircase,shifts,R,tau,upsilon;
local col_tau,row_upsilon,elem_tauupsilon,r,col_next_mon,rtau;
local excluded_cols,excluded_rows,i,g;
    userinfo (1,AdaptiveScalarFGLM,`*** Tab `, Tab,
              `in variables`, diffvars, `and`, vars,
              `and cone generators`,cone_gen,
              `up to degree`, d,
              `in characteristic`, p, `***`);
    elem_tauupsilon:= fromMonToTabElem (Tab,vars,1);
    userinfo (2,AdaptiveScalarFGLM,`***\t New monomial`, 1, `***`);
    if (elem_tauupsilon = 0) then
        userinfo (2,AdaptiveScalarFGLM,`***\t New potential relation`, 1,
                  `***`);
        if (full_answer) then return [[1,1,1]];
        else return [1];
        end if:
    end if:
    userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
    L         := Matrix (1,1,1);
    U         := Matrix (1,1,elem_tauupsilon);
    rank      := 1;
    next_mon  := sortListMon ([op(cone_gen),op(diffvars)],mon_ord);
    next_shift:= sortListMon (cone_gen,mon_ord);
    staircase := [1];
    shifts    := [1];
    excluded_cols:= to_exclude;
    excluded_rows:= select (m->indets (m) subset {op(vars)},to_exclude);
    if (d = 1) then
        userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
        R:= map (v->[v,redp (v-fromMonToTabElem (Tab,vars,v)/elem_tauupsilon,p),
                     1],next_mon);
        if (full_answer) then return R;
        else return map (r->r[2],R);
        end if:
    end if:
    R        := [];
    while (not indets (next_mon) subset {op(diffvars)} and next_shift <> []) do
        tau    := next_mon[1];
        upsilon:= next_shift[1];
        userinfo (2,AdaptiveScalarFGLM,`***\t New monomials`, tau,
                  `and`, upsilon,`***`);
        col_tau    := from2QCSetsToMultiHankelMatrix (Tab,vars,diffvars,
                                                      shifts,[tau],p);
        row_upsilon:= from2QCSetsToMultiHankelMatrix (Tab,vars,diffvars,
                                                      staircase,[upsilon],p);
        if (p = 0) then
            col_tau    := LinearSolve (L,col_tau);
            row_upsilon:= Transpose (LinearSolve (Transpose (U),row_upsilon));
        elif (p < 2^32) then
            L          := LinearAlgebra[Modular][Mod] (p,L,integer[8]);
            U          := LinearAlgebra[Modular][Mod] (p,U,integer[8]);
            col_tau    := LinearAlgebra[Modular][Mod] (p,col_tau,integer[8]);
            row_upsilon:= LinearAlgebra[Modular][Mod] (p,row_upsilon,integer[8]);
            col_tau    := LinearAlgebra[Modular][LinearSolve] (p,<L|col_tau>,
                                                               1,inplace=false);
            row_upsilon:= Transpose (LinearAlgebra[Modular][LinearSolve]
                                     (p,<Transpose(U)|row_upsilon>,
                                      1,inplace=false));
        else
            L          := L mod p;
            U          := U mod p;
            col_tau    := LinearSolve (L,col_tau mod p) mod p;
            row_upsilon:= Transpose (LinearSolve (Transpose (U),row_upsilon mod p)
                                     mod p);
        end if:
        elem_tauupsilon:= fromQCMonToTabElem (Tab,vars,diffvars,tau*upsilon);
        elem_tauupsilon:= redp (-row_upsilon.col_tau+elem_tauupsilon,p);
        if (type (elem_tauupsilon,Matrix)) then elem_tauupsilon:= elem_tauupsilon[1,1];
        end if:
        if (elem_tauupsilon <> 0) then
            userinfo (2,AdaptiveScalarFGLM,`***\t\t in the staircase! ***`);
            L        := <<L,row_upsilon>|<Matrix (rank,1),1>>;
            U        := <<U,Matrix (1,rank)>|<col_tau,elem_tauupsilon>>;
            staircase:= [op(staircase),tau];
            shifts   := [op(shifts),upsilon];
            next_mon := MakeUnique (
                sortListMon ([op(next_mon[2..-1]),
                              op(remove (l->`or` (seq (divide (l,g,'q') and
                                                       isInTheCone (q,vars,cone_gen),
                                                       g in excluded_cols)),
                                         [op(expand (tau*cone_gen)),
                                          op(expand (tau*diffvars))]))],
                             mon_ord),1);
            next_shift:= MakeUnique (
                sortListMon ([op(next_shift[2..-1]),
                              op(remove (l->`or` (seq (divide (l,g,'q') and
                                                       isInTheCone (q,vars,cone_gen),
                                                       g in excluded_rows)),
                                         expand (upsilon*cone_gen)))],
                             mon_ord),1);
            rank:= rank+1;
            if (rank >= d) then
                userinfo (2,AdaptiveScalarFGLM,`***\t Early termination ***`);
                col_next_mon:= from2QCSetsToMultiHankelMatrix (Tab,vars,diffvars,
                                                               shifts,
                                                               next_mon,p);
                if (p = 0) then
                    col_next_mon:= LinearSolve (U,
                                                LinearSolve (L,col_next_mon));
                elif (p < 2^32) then
                    L:= LinearAlgebra[Modular][Mod] (p,L,integer[8]);
                    U:= LinearAlgebra[Modular][Mod] (p,U,integer[8]);
                    col_next_mon:= LinearAlgebra[Modular][Mod] (p,col_next_mon,
                                                                integer[8]);
                    col_next_mon:= LinearAlgebra[Modular][LinearSolve] (p,
                                                                        <L|
                                                                         col_next_mon>,
                                                                        nops(next_mon),
                                                                        inplace=false);
                    col_next_mon:= LinearAlgebra[Modular][LinearSolve] (p,
                                                                        <U|
                                                                         col_next_mon>,
                                                                        nops(next_mon),
                                                                        inplace=false);
                else
                    L           := L mod p;
                    U           := U mod p;
                    col_next_mon:= col_next_mon mod p;
                    col_next_mon:= LinearSolve (L,col_next_mon) mod p;
                    col_next_mon:= LinearSolve (U,col_next_mon) mod p;
                end if:
                r:= nops (R);
                R:= [op(R),seq ([next_mon[i],
                                 redp (next_mon[i]
                                       - DotProduct (Column (col_next_mon,i),
                                                     staircase),p),
                                 shifts],
                                i=1..nops(next_mon))];
                map (r->userinfo (2,AdaptiveScalarFGLM,`***\t New potential relations`,
                          r[2],`***`),R[r+1..-1]);
                if (full_answer) then return R;
                else return map (r->r[2],R);
                end if:
            end if:
        else
            userinfo (2,AdaptiveScalarFGLM,`***\t\t NOT in the staircase! ***`);
            if (p = 0) then
                col_tau:= LinearSolve (U,col_tau);
            elif (p < 2^32) then
                U      := LinearAlgebra[Modular][Mod] (p,U,integer[8]);
                col_tau:= LinearAlgebra[Modular][Mod] (p,col_tau,integer[8]);
                col_tau:= LinearAlgebra[Modular][LinearSolve] (p,<U|col_tau>,
                                                               1,inplace=false);
            else
                U      := U mod p;
                col_tau:= col_tau mod p;
                col_tau:= LinearSolve (U,col_tau) mod p;
            end if:
            rtau:= redp (tau - DotProduct (col_tau,staircase),p);
            R   := [op(R),[tau, rtau, [op(shifts),upsilon]]];
            next_mon  := remove (l->divide (l,tau,'q') and isInTheCone (q,vars,
                                                                        cone_gen),
                                 next_mon);
            excluded_cols:= [op(excluded_cols),tau];
            if (indets (tau) subset {op(vars)}) then
                excluded_rows:= [op(excluded_rows), tau];
            end if:
            # next_shift:= remove (l->divide (l,upsilon,'q') and isInTheCone (q,vars,
            #                                                                 cone_gen),
            #                      next_shift);
            userinfo (2,AdaptiveScalarFGLM,`***\t New potential relation`, rtau,
                      `***`);
        end if:
    end do:
    if (next_shift = []) then
        userinfo (1,AdaptiveScalarFGLM,`***\t Not enough rows`,`***`);
    end if:
    if (full_answer) then return R;
    else return map (r->r[2],R);
    end if:
end proc:

ScalarNumerator2:= proc (UTMsV,P,w,i,ai,Ls)
description "Returns the Scalar Numerator Omega (Tab,P), "
    "See Hyun et al. 2020, where m=1.";

end proc:

BlockParametrization:= proc (Tab,vars:=['x','y'],lf:=vars[-1],
                             extravar:=vars[-1],d:=10,p:=mychar)
description "Computes the parametrizations of the radical of the ideal "
    "of relations of Tab, assuming it is in shape position for "
    "LEX (vars[1]>...>vars[-1]>T) where T is a new variable satisfying "
    "T=lf, a linear form in vars[1],...,vars[-1]. See Hyun et al. 2020.";
local P,Q,n,C,k,lfpow,lfTab,st,i;
    st:= time ();
    n:= nops (vars);
    lfpow:= [seq (redp (expand (lf^i),p),i=0..2*d-1)];
    lfTab:= proc (i)
    option remember;
        return fromPolToTabRelation (Tab,vars,lfpow[i+1],p);
    end proc:
    SeqVector (lfTab,2*d);
    # printf("Sequence computed\n"):
    P:= PolynomialBerlekampMassey (lfTab,extravar,extravar^(2*d-1),p,false);
    # printf("Elimination polynomial computed\n"):
    if (p = 0) then Q:= sqrfree (P);
    else Q:= Sqrfree (P) mod p;
    end if:
    Q:= redp (expand (Q[1]*mul (Q[2][i][1],i=1..nops(Q[2]))),p);
    C:= table ([]);
    C[1]:= ScalarNumerator (lfTab,extravar,d,P,p);
    # print ("C1=",C[1]);
    if (p = 0) then gcdex (C[1],Q,extravar,'invC1');
    else Gcdex (C[1],Q,extravar,'invC1') mod p;
    end if:
    # print ("invC1=",invC1);
    for k from 1 to n do
        lfTab:= proc (i)
        option remember;
            return fromPolToTabRelation (Tab,vars,lfpow[i+1]*vars[k],p);
        end proc:
        C[vars[k]]:= redp (expand (ScalarNumerator
                                  (lfTab,extravar,d,P,p)*invC1),p);
        if (p = 0) then C[vars[k]]:= rem (C[vars[k]],Q,extravar);
        else C[vars[k]]:= Rem (C[vars[k]],Q,extravar) mod p;
        end if:
    end do:
    # printf ("elapsed time: %f\n",time()-st);
    return [Q,seq (vars[n-k]-C[vars[n-k]] mod p,k=0..n-1),P];
end proc: