#
# Author: Felipe Cucker 1991
#
constants:=constants, HASH;

# sth computes the Sturm-Habicht sequence of two polynomials w.r.t. a 
# particular variable. Their inputs are:
#	D - a univariate polynomial ring (over an arbitrary ring)
#	p,q - two elements of D

sth:= proc(D,p,q)
    local suc,n1,n2,lp,lq,r,R,S0,C0,S1,C1,d,NS,Sr,j,O,n;
    if not hasCategory(D,UnivariatePolynomial) then 
	ERROR(`1st argument must be a PolynomialRing`)
    fi;
    O:=D[CoefficientRing];
    if D[`=`](q,D[0]) then RETURN(p)
    fi;
    n1:=D[Degree](p);
    n2:=D[Degree](q);
    lp:=D[Lcoeff](p);
    lq:=D[Lcoeff](q);
    suc:=p,q;
    if n1>n2 then n:=n1-1;
         d:=n-n2;
         if d<>0 then
             C1:= O[`^`](O[`*`](lp,lq),d);
             S1:=parity(D,(d*(d+1)/2),D[`.`](C1,q));
             C0:=O[`^`](lp,d);
             S0:=parity(D,(d+2)*(d+3)/2),D[`.`](C0,D[PseudoRem](p,q));
             if D[`<>`](S0,D[0]) then
                  suc:=suc,S1,S0
             else suc:=suc,S1
             fi;
          else S1:=q;
             S0:= D[`-`](D[PseudoRem](p,q));
             if D[`<>`](S0,D[0]) then suc:= suc,S0
             fi;
          fi;
          j:=n2-1;
    else
        n:=n2;
        S1:=q;
        j:=n-1;
        if n1=n2 then 
            S0:=D[`-`](D[PseudoRem](D[`.`](lq,p),q))
        else
            S0:=D[`-`](D[`.`](O[`^`](lq,2),p))
        fi;
        if D[`<>`](S0,D[0]) then suc:=suc,S0
        fi;
    fi;
    R:=D[Lcoeff](S1);
    if D[`=`](S0,D[0]) then r:=-1
            else r:=D[Degree](S0);
    fi;
    while r>=0 do
      if r>0 then
         if r=j then 
             NS:=D[Div](D[PseudoRem](S1,S0),
		D[Constant](O[`^`](O[`-`](R),j-r+2)));
             NS:=parity(D,(j-r+2)*(j-r+3)/2,NS);
             if D[`<>`](NS,D[0]) then suc:=suc,NS
             fi;
             S1:=S0;
             S0:=NS;
             R:=D[Lcoeff](S1);
             j:=r-1;
          else 
             Sr:=D[Div](D[PseudoRem](S1,S0),D[Constant](O[`^`](R,j-r+2)));
             C1:=O[`^`](D[Lcoeff](S0),j-r);
             S1:=D[Div](D[`.`](C1,S0),D[Constant](O[`^`](R,j-r)));
             S1:=parity(D,(j-r)*(j-r+1)/2,S1);
             S0:=parity(D,(j-r+2)*(j-r+3)/2,Sr);
             if D[`<>`](S0,D[0]) then suc:=suc,S1,S0
                 else suc:=suc,S1
             fi;
             R:=D[Lcoeff](S1);
             j:= r-1;
         fi;
       else 
         if j>0 then
           C1:=O[`^`](D[Lcoeff](S0),j);
           S1:=D[Div](D[`.`](C1,S0),D[Constant](O[`^`](R,j)));
           S1:=parity(D,j*(j+1)/2,S1);
           suc:=suc,S1;
         fi;
         S0:=D[0];
      fi;
      if D[`<>`](S0,D[0]) then r:=D[Degree](S0)
             else r:=-1
      fi;
    od;
    suc;
end:


# parity is a private function that computes a*(-1)^n in an arbitrary ring D.

parity:=proc(D,n,a);
    if irem(n,2)=0 then a 
    else D[`-`](a)
    fi;
end:
 

# ss computes the difference of sign variations in -infty. and +infty of the 
# Sturm-Habicth sequence of two polynomials over an ordered ring.

ss:= proc(D,p,q)
     local s;
     s:= sth(D,p,q);
     difvar(D,s);
end:


# nroots computes the number of roots of a polynomial over an ordered ring in 
# the real closure of that ring.

nroots:= proc(D,p)
    local s,O;
    if not hasCategory(D,UnivariatePolynomial) then 
	ERROR(`1st argument must be a PolynomialRing`)
    fi;
    O:=D[CoefficientRing];
    if not hasCategory(O,OrderedDomain) then 
	ERROR(`1st argument must be over an OrderedDomain`)
    fi;
    if D[Degree](p)<1 then
        if D[`=`](p,D[0]) then RETURN(p, `is the zero polynomial`)
             else RETURN(0)
        fi;
    fi;
    s:=sth(D,p,D[Diff](p));
    difvar(D,s);
end:


# difvar computes the difference of sign variations of a sequence of polynomials 
# in -infty. and +infty. Its arguments are:
#	D - a univariate polynomial ring over an ordered domain
#	s - a sequence of polynomials in D

difvar:=proc(D)
    local suc1,suc2,i,s,O;
    if not hasCategory(D,UnivariatePolynomial) then 
	ERROR(`1st argument must be a PolynomialRing`)
    fi;
    O:=D[CoefficientRing];
    if not hasCategory(O,OrderedDomain) then 
	ERROR(`1st argument must be over an OrderedDomain`)
    fi;
    suc1:=NULL;
    suc2:=NULL;
    for i from 2 to nargs do
        if D[`<>`](args[i],D[0]) then
            s:=O[Sign](D[Lcoeff](args[i]));
            suc2:=suc2,s;
            s:=s*(1-2*(irem(D[Degree](args[i]),2)));
            suc1:=suc1,s;
        fi;
    od;
    var(suc1)-var(suc2);
end:

var:= proc()
    local i,k,v;
    k:=0;
    v:=0;
    if nargs>1 then
    for i from 1 to nargs do
         if args[i]>0 then
              if k=0 then k:=1;
                   elif k<0 then v:= v+1; k:=1;
              fi;
         elif args[i]<0 then
              if k=0 then k:=-1;
                   elif k>0 then v:=v+1; k:=-1;
              fi;
         fi;
    od;
    fi;
    v;
end:


# kron(A,B) returns the kroneker product of the matrixes A and B.

kron:=proc(A,B)
    local i,j,k,l,n,d,C;
    d:=linalg[rowdim](A);
    n:=linalg[rowdim](B);
    C:=array(1..n*d,1..n*d);
    for i from 1 to n do
        for j from 1 to n do
            for k from 1 to d do
 		for l from 1 to d do
		    C[d*(i-1)+k,d*(j-1)+l]:=A[k,l]*B[i,j];
		od;
            od;
        od;
    od;
    C;
end:


# clean(A,n,L) returns a submatrix of A having the same rank, obtained by 
# eliminating some rows. A second parameter must be given to the procedure 
# with the number of the first rows we indeed know to be lin. indep. In the 
# third parameter (which must be unevaluated) the list of conserved rows is 
# returned.

clean:=proc(A,n,L)
    local C,E,l,i,j,F,k,s,rank,r;
    r:=linalg[coldim](A);
    if n>r then ERROR(`wrong size of arguments`)
    fi;
    C:=array(1..r,1..r);
    E:=array(1..r,1..r);
    l:=[];
    for i from 1 to n do
         for j from 1 to r do
              C[i,j]:=A[i,j];
         od;
    l:=[op(l),i];
    od;
    F:=linalg[gausselim](linalg[submatrix](C,1..n,1..r));
    for i from 1 to n do
         for j from 1 to r do
              E[i,j]:=F[i,j];
         od;
    od;
    k:=n;
    for i from n+1 while k<r do
         for j from 1 to r do
         E[k+1,j]:=A[i,j];
         od;
         linalg[gausselim](linalg[submatrix](E,1..k+1,1..r),'rank');
         if rank=k+1 then k:=k+1;
                        l:=[op(l),i];
                        for s from 1 to r do
                             C[k,s]:=A[i,s];
                        od;
         fi;
    od;
    L:=l;
    C;
end:


# reduce(D,p) takes out from p, a polynomial over an ordered domain, all the 
# coefficients that vanish for the ordering in that domain. Its arguments are:
#	D - a univariate polynomial ring over an integer domain
#	p - an element of D.

reduce:=proc(D,p)
    local O,np,d,C,vanish;
    if not hasCategory(D,UnivariatePolynomial) then 
        ERROR(`1st argument must be a polynomial ring`)
    fi;
    if not D[Type](p) then
        ERROR(`2nd argument must be a polynomial`)
    fi;
    O:=D[CoefficientRing];
    np:=p;
    C:=D[Lcoeff](np);
    vanish:=O[Zero](C);
    d:=D[Degree](np);
    while (vanish and d>0) do
        np:=D[Reductum](np);
        C:=D[Lcoeff](np);
        vanish:=O[Zero](C);
        d:=D[Degree](np);
    od;
    if vanish then np:=D[Reductum](np) 
    fi;
    np;
end:


# squarefree(D,p) squarefrees the polynomial p belonging to the ring of 
# univariated polynomials over an ordered domain. 

squarefree:=proc(D,p)
    local suc,R,n,M,O;
    O:=D[CoefficientRing];
    R:=reduce(D,p);
    suc:=[sth(D,R,D[Diff](R))];
    n:=nops(suc);
    M:=reduce(D,op(n,suc));
    while D[`=`](M,D[0]) do
        n:=n-1;
        M:=reduce(D,op(n,suc));
    od;
    if D[Degree](M)>0 then
        M:=D[Primpart](M);
        R:=D[`.`](O[`^`](D[Lcoeff](R),D[Degree](R)-D[Degree](M)+1),R);
        R:=D[Div](R,M);
    fi;
    R:=D[Primpart](R);
end:
    
# SIadd adds the information of a polynomial to a Tarski Data Type. 
# They parameters are:
#	D - a univariate polynomial ring over an ordered domain
#	T - a Tarski Data Type
# 	Q - an element of D
#	m - (optional) the number of multiple roots of the main polynomial in T.

SIadd:=proc(D,T,Q,m)
    local O,P,P1,suc,M,d,V,E,s,r,ns,ns2,i,NE,NE1,NE2,NE3,NV,NR,
          Mat,A,NC,l,l1,Nr,EF,CF,B,L,VF,j,RF,mult,NA,R2; 
    if not hasCategory(D,UnivariatePolynomial) then 
        ERROR(`1st argument must be a polynomial ring`)
    fi;
    if not D[Type](Q) then
        ERROR(`3rd argument must be a polynomial`)
    fi; 
    O:=D[CoefficientRing];
    if not hasCategory(O,OrderedDomain) then 
	ERROR(`1st argument must be over an OrderedDomain`)
    fi;
    P:=op(1,T);
    if not D[Type](P) then
        ERROR(`incorrect 2nd argument`)
    fi;
    if  D[Degree](Q)<1 then
	    E:=op(3,T);
	    r:=nops(E);
            if D[`=`](Q,D[0]) then s:=0 else s:=O[Sign](D[Coeff](Q,0))
	    fi;
	    NE:=[];
	    for i from 1 to r do
		NE:=[op(NE),[op(op(i,E)),s]]
	    od;
	    RETURN([op(1,T),op(2,T),NE,op(4,T),op(5,T),op(6,T)])
    fi;
    if op(3,T)=[] then RETURN(T)
    fi;
    P1:=D[Diff](P);
    suc:=sth(D,P,D[`*`](P1,Q));
    M:=op(nops([suc]),[suc]);
    d:=D[Degree](M);
    V:=op(5,T);
    E:=op(3,T);
    s:=V[1];
    r:=linalg[vectdim](V);
    ns:=difvar(D,suc);
    if s=ns then NE:=[];
	 for i from 1 to r do
	      NE:=[op(NE),[op(op(i,E)),1]]
	 od;
	 RETURN([P,op(2,T),NE,op(4,T),op(5,T),op(6,T)]);
    elif s=-ns then NE:=[];
	 for i from 1 to r do
	      NE:=[op(NE),[op(op(i,E)),-1]]
 	 od;
	 RETURN([P,op(2,T),NE,op(4,T),op(5,T),op(6,T)]);
    fi;
    NR:=op(2,T);
    if d=0 then 
        NE1:=[];
        NE2:=[];
        NV:=array(1..2*r);
        for i from 1 to r do
            NE1:=[op(NE1),[op(op(i,E)),1]];
            NE2:=[op(NE2),[op(op(i,E)),-1]];
            NV[i]:=V[i];
       	    NR:=[op(NR),D[`*`](op(i,NR),Q)];
        od;
        NE:=[op(NE1),op(NE2)];
        NV[r+1]:=ns;
        for i from 2 to r do
	    NV[r+i]:=ss(D,P,D[`*`](P1,op(i+r,NR)));
        od;
        Mat:=array(1..2,1..2);
        Mat[1,1]:=1;
        Mat[1,2]:=1;
        Mat[2,1]:=1;
        Mat[2,2]:=-1;
    else 
	if nargs=4 then mult:=args[4]
	    else mult:=nroots(D,D[Gcd](P,P1)) 
	fi;
      	if mult=0 then
	    ns2:=s - nroots(D,M)
	else
	    ns2:=ss(D,P,D[`*`](P1,D[`^`](Q,2)))
        fi;
        if ns2=0 then NE:=[];
	    for i from 1 to r do
	       NE:=[op(NE),[op(op(i,E)),0]]
	    od;
	    RETURN([P,op(2,T),NE,op(4,T),op(5,T),op(6,T)]);
        elif ns=ns2 then 
            NE1:=[];
            NE2:=[];
            NV:=array(1..2*r);
            for i from 1 to r do
                NE1:=[op(NE1),[op(op(i,E)),0]];
                NE2:=[op(NE2),[op(op(i,E)),1]];
                NV[i]:=V[i];
                NR:=[op(NR),D[`*`](op(i,NR),Q)];
            od;
            NE:=[op(NE1),op(NE2)];
            NV[r+1]:=ns;
            for i from 2 to r do
	        NV[r+i]:=ss(D,P,D[`*`](P1,op(i+r,NR)));
            od;
	    Mat:=array(1..2,1..2);
            Mat[1,1]:=1;
            Mat[1,2]:=1;
            Mat[2,1]:=0;
            Mat[2,2]:=1;
        elif ns=-ns2 then
            NE1:=[];
            NE2:=[];
            NV:=array(1..2*r);
            for i from 1 to r do
                NE1:=[op(NE1),[op(op(i,E)),0]];
                NE2:=[op(NE2),[op(op(i,E)),-1]];
                NV[i]:=V[i];
                NR:=[op(NR),D[`*`](op(i,NR),Q)];
            od;
            NE:=[op(NE1),op(NE2)];
            NV[r+1]:=ns;
            for i from 2 to r do
	        NV[r+i]:=ss(D,P,D[`*`](P1,op(i+r,NR)));
            od;
	    Mat:=array(1..2,1..2);
            Mat[1,1]:=1;
            Mat[1,2]:=1;
            Mat[2,1]:=0;
            Mat[2,2]:=-1;
        else NE1:=[];
            NE2:=[];
	    NE3:=[];
            R2:=[];
            NV:=array(1..3*r);
            for i from 1 to r do
                NE1:=[op(NE1),[op(op(i,E)),0]];
                NE2:=[op(NE2),[op(op(i,E)),1]];
                NE3:=[op(NE3),[op(op(i,E)),-1]];
                NV[i]:=V[i];
	        R2:=[op(R2),D[`*`](op(i,NR),D[`^`](Q,2))];
                NR:=[op(NR),D[`*`](op(i,NR),Q)];
            od;
            NE:=[op(NE1),op(NE2),op(NE3)];
            NR:=[op(NR),op(R2)];
            NV[r+1]:=ns;
	    NV[2*r+1]:=ns2;
            for i from 2 to r do
	        NV[r+i]:=ss(D,P,D[`*`](P1,op(i+r,NR)));
	        NV[2*r+i]:=ss(D,P,D[`*`](P1,op(i+2*r,NR)));
            od;
            Mat:=array(1..3,1..3);
            Mat[1,1]:=1;
            Mat[1,2]:=1;
            Mat[1,3]:=1;
            Mat[2,1]:=0;
            Mat[2,2]:=1;
            Mat[2,3]:=-1;
            Mat[3,1]:=0;
            Mat[3,2]:=1;
            Mat[3,3]:=1;
        fi;
    fi;
    A:=kron(op(6,T),Mat);
    NC:=linalg[linsolve](A,NV);
    l:=[];
    l1:=[];
    Nr:=0;
    EF:=[];
    CF:=[];
    for i from 1 to linalg[vectdim](NV) do
        l1:=[op(l1),i];
        if NC[i]<>0 then Nr:=Nr+1;
            EF:=[op(EF),op(i,(NE))];
            CF:=[op(CF),NC[i]];
            l:=[op(l),i];
        fi;
    od;
    if Nr<linalg[vectdim](NV) then
        B:=linalg[submatrix](A,l1,l);
        NA:=clean(B,r,'L');
    else NA:=A;
        L:=l1;
    fi;
    VF:=array(1..Nr);
    RF:=[];
    for i from 1 to Nr do
        j:=op(i,L);
        VF[i]:=NV[j];
        RF:=[op(RF),op(j,NR)];
    od;
    [P,RF,EF,CF,VF,NA];
end:


# SI(D,P,Q1,...,Qk) returns a TDT with the signs of the Qi's on the roots 
# of P. All P,Q1,...,Qk belongs to D, a univariate polynomial ring over an 
# orderd domain.

SI:=proc(D)
    local O,k,P,suc,r,M,A,V,C,NT,i;
    k:=nargs;
    P:=args[2];
    if not D[Type](P) then
         ERROR(`arguments must be polynomials`)
    fi;
    if D[Degree](P)=0 then RETURN(SIadd(D,P,D[1]))
    fi;
    O:=D[CoefficientRing];
    if not hasCategory(O,OrderedDomain) then 
	ERROR(`1st argument must be over an OrderedDomain`)
    fi;    
    suc:=sth(D,P,D[Diff](P));
    r:=difvar(D,suc);
    if r=0 then RETURN([])
    fi;
    M:=op(nops([suc]),[suc]);
    if D[Degree](M)>0 then 
        M:=D[Primpart](M);
        P:=D[`.`](O[`^`](D[Lcoeff](P),D[Degree](P)-D[Degree](M)+1),P);
        P:=D[Div](P,M);
    fi;
    A:=array(1..1,1..1);
    A[1,1]:=1;
    V:=array(1..1);
    V[1]:=r;
    C:=array(1..1);
    C[1]:=r;
    NT:=[P,[1],[[]],C,V,A];
    if k=1 then RETURN(NT)
    else for i from 3 to k do
             NT:=SIadd(NT,args[i],0)
         od;
    fi;
    NT;
end:


# sqfran(D,R,nder) computes the Thom code of the roots of a polynomial R
# jumping to the first derivative once the roots are separated in groups 
# of at most two roots. Its parameters are:
#	D - a univariate polynomial ring over an ordered domain
#	R - an element of D
#	nder - an unevaluated name where the number of used derivatives is
#	       returned

sqfran:=proc(D,R,nder)
    local O,suc,r,A,V,C,M,NT,d,E,NE,P,i,der;
    if not hasCategory(D,UnivariatePolynomial) then 
	ERROR(`1st argument must be a polynomial ring`)
    fi;
    O:=D[CoefficientRing];
    if not hasCategory(O,OrderedDomain) then 
	ERROR(`1st argument must be over an OrderedDomain`)
    fi;
    if not D[Type](R) then
        ERROR(`2nd argument must be a polynomial`)
    fi;
    P:=R;
    A:=array(1..1,1..1);
    A[1,1]:=1;
    V:=array(1..1);
    V[1]:=0;
    C:=array(1..1);
    C[1]:=0; 
    if D[Degree](P)=0 then 
        RETURN([P,[1],[],C,V,A])
    fi;
    suc:=sth(D,P,D[Diff](P));
    r:=difvar(D,suc);
    V[1]:=r;
    C[1]:=r;
    if r=0 then RETURN([P,[1],[],C,V,A])
        elif r=1 then 
            if nargs>2 then nder:=0;
            fi;
            RETURN([P,[1],[[]],C,V,A])
    fi;
    M:=op(nops([suc]),[suc]);
    if D[Degree](M)>0 then
        M:=D[Primpart](M);
        P:=D[`.`](O[`^`](D[Lcoeff](P),D[Degree](P)-D[Degree](M)+1),P);
        P:=D[Div](P,M);
    fi;
    P:=D[Primpart](P);
    NT:=[P,[1],[[]],C,V,A];
    d:=D[Degree](P)-1;
    while r>1 do 
        if (r=2 and d>1) then E:=op(3,NT);
            NE:=[];
            for i from 1 to nops(E) do
                NE:=[op(NE),[op(op(i,E)),HASH]]
            od;
            NT:=[op(1,NT),op(2,NT),NE,op(4,NT),op(5,NT),op(6,NT)];
            NT:=SIadd(D,NT,D[Diff](P),0);
        else
            der:=P;
            for i from 1 to d do
                der:=D[Diff](der)
            od;
            der:=D[Div](der, D[Coerce](d!));
            der:=D[Primpart](der);
            NT:=SIadd(D,NT,der,0)
        fi;
        d:=d-1;
        r:=maximo(op(4,NT));
    od;
    if nargs>2 then nder:=D[Degree](P)-d-1
    fi;
    NT;
end:

maximo:=proc(v)
    local m,i,n;
    n:=linalg[vectdim](v);
    m:=v[1];
    if n>1 then
         for i from 2 to n do
              if v[i]>m then m:=v[i]
              fi;
         od;
    fi;
    m;
end:


# ransi(D,T,Q1,...,Qk) adds the information of the polynomials Q1,...,Qk 
# to the TDT T, without changing the order of the elements coded by T.

ransi:=proc(D)
    local n,T,i;
    if not hasCategory(D,UnivariatePolynomial) then 
	ERROR(`1st argument must be a polynomial ring`)
    fi;
    O:=D[CoefficientRing];
    if not hasCategory(O,OrderedDomain) then 
	ERROR(`1st argument must be over an OrderedDomain`)
    fi;    
    n:=nargs;
    if n<2 then ERROR(`incorrect number of arguments`)
    elif (n=2 or op(3,args[2])=[]) then RETURN(args[2])
    fi;
    T:=mark(args[2]);
    for i from 3 to n do
        T:=SIadd(D,T,args[i])
    od;
    T:=unmark(T);
end:


mark:=proc(T)
    local E,NE,i;
    E:=op(3,T);
    NE:=[];
    for i from 1 to nops(E) do
    	NE:=[op(NE),[i,op(op(i,E))]]
    od;
    [op(1,T),op(2,T),NE,op(4,T),op(5,T),op(6,T)];
end:

unmark:=proc(T)
    local E,A,i,l,NE,j,k,found,rt,NA,NC,C,h,s;
    E:=op(3,T);
    s:=nops(E);
    NE:=[];
    A:=op(6,T);
    C:=op(4,T);
    NC:=[];
    NA:=array(1..s,1..s);
    for l from 1 to s do
        i:=s-l+1;
    	found:=false;
	j:=s+1;
	while not found do 
	    j:=j-1;
	    if op(1,op(j,E))=i then 
   		found:=true;
		rt:=[];
		for k from 2 to nops(op(j,E)) do
		    rt:=[op(rt),op(k,op(j,E))]
		od;
		NE:=[rt,op(NE)];
		for h from 1 to s do
		    NA[h,i]:=A[h,j]
		od;
		NC:=[op(NC),op(j,C)];
	    fi;
	od;
    od;
    [op(1,T),op(2,T),NE,NC,op(5,T),NA];
end: 

# minor(D,L1,L2,c,h) compares two Thom codes L1 and L2 and returns true if 
# L1 is smaller than L2 or the two roots are consecutive (and its relative 
# position depends on the other coded roots). Its inputs are:
#	D - a univariate polynomial ring over an ordered domain
#	L1,L2 - two Thom codes
#	c - the leading coefficient of the polynomial coding L1 and L2
#	h - (optional) an unevaluated name that returns true when some more 
#	   information is needed to decide the sorting.

minor:=proc(D,L1,L2,c,h)
    local O,i,j,M;
    O:=D[CoefficientRing];
    j:=0;
    if nargs =5 then h:=false
    fi;
    for i from 1 while j=0 do
        if (op(i,L1)=HASH or (i=nops(L1) and op(i,L1)=op(i,L2))) 
            then j:=1;
            if nargs =5 then h:=true
            fi;
            M:=true
        elif op(i,L1)<op(i,L2) then j:=1;
            if i>1 then 
                if op(i-1,L1)>0 then M:=true
                      else M:=false
                fi;
            else 
                if O[Sign](c)=1 then M:=true
                     else M:=false
                fi;
            fi;
        elif op(i,L1)>op(i,L2) then j:=1;
            if i>1 then 
                if op(i-1,L1)<0 then M:=true
                      else M:=false
                fi;
            else 
                if O[Sign](c)=-1 then M:=true
                     else M:=false
                fi;
            fi;
        fi;
    od;
    M;
end:


# ransort(D,T) sorts the roots coded by the TDT T, whose main polynomial 
# belongs to D.

ransort:=proc(D,T)
    local O,v,i,j,U,s1,n,k,E,P,c,s;
    O:=D[CoefficientRing];
    E:=op(3,T);
    P:=op(1,T);
    c:=D[Lcoeff](P);
    n:=nops(E);
    if n<2 then RETURN(E)
    fi;
    v:=array(1..n);
    for i from 1 to n do
         v[i]:=op(i,E);
    od;
    for i from 1 to n-1 do
        for j from i+1 to n do
            if minor(D,v[j],v[i],c) then U:=v[j];
                                     v[j]:=v[i];
                                     v[i]:=U;
            fi;
        od;
    od;
    i:=0;
    while i<n-1 do
         k:=0;
         for j from 1 to min(nops(v[n-i]),nops(v[n-i-1])) while k=0 do
              if op(j,v[n-i])<>op(j,v[n-i-1]) then
                   k:=1;
                   i:=i+1;
                                              else
                   if op(j,v[n-i])=HASH then
                        k:=1;
                        s1:=(irem(i,2)*2-1)*op(j+1,v[n-i]);
                        if s1*O[Sign](c)>0 then U:=v[n-i];
                                 v[n-i]:=v[n-i-1];
                                 v[n-i-1]:=U;
                                 i:=i+2
                                else i:=i+1
                        fi;
                   fi;
              fi;
         od;
    od;
    s:=NULL;
    for i from 1 to n do
         s:=s,v[i];
    od;
    s:=[s];
end:


# position(rt,T) returns the position of the Thom code rt in T, a TDT. It 
# returns 0 if rt is not in T.

position:=proc(rt,T)
    local E,l,r,i,reached,rt1,l1,coincide,k;
    E:=op(3,T);
    l:=nops(rt);
    r:=nops(E);
    i:=1;
    reached:=false;
    while (not reached and i<=r) do
	rt1:=op(i,E);
	l1:=nops(rt1);
	if l=l1 then
	    k:=1;
	    coincide:=true;
	    while (coincide and k<=l) do 
		if op(k,rt)<>op(k,rt1) then 
		    coincide:=false
		fi;
		k:=k+1;
	    od;
	    if coincide then reached:=true
	    fi;
	fi;
	i:=i+1;
    od;
    if (not reached) then i:=1;
    fi;
    i-1;
end:

save `IF.m`;
quit
