with(Groebner);

##### minus für Listen #########
local `minus`:= overload([

     proc(L::list, E::{list,set})
     option overload;
     local R:= table(sparse), e;

          for e in E do R[e]:= R[e]+1 end do;
          remove(e-> if R[e] > 0 then R[e]:= R[e]-1; true else false end if, L)
     end proc,

     :-`minus`
]):

########### help function: sort by signature, then by lt(v) ##########
FHelp:=proc(pp,qq) ### e.g. pp=[[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar]
global Vars, tord, ord;
if ord = 1 then ## pot-lift
	if pp[1][3]<qq[1][3] then
		return true;
	elif qq[1][3]<pp[1][3] then
		return false;
	fi;
	if not TestOrder(pp[1][2],qq[1][2], tord) then
		return false;
	elif not TestOrder(qq[1][2],pp[1][2], tord) then
		return true;
	fi; ### so signatures are equal
	if not TestOrder(pp[2][2],qq[2][2], tord) then
		return false;
	elif not TestOrder(qq[2][2],pp[2][2], tord) then
		return true;
	fi;
elif ord = 2 then ## top lift
	if not TestOrder(pp[1][2],qq[1][2], tord) then
		return false;
	elif not TestOrder(qq[1][2],pp[1][2], tord) then
		return true;
	fi; ### so signatures might only differ in position
	if pp[1][3]<qq[1][3] then
		return true;
	elif qq[1][3]<pp[1][3] then
		return false;
	fi;### now they are the same
	if not TestOrder(pp[2][2],qq[2][2], tord) then
		return false;
	elif not TestOrder(qq[2][2],pp[2][2], tord) then
		return true;
	fi;
fi;
return true; ### for syzygies
end:
	
### find the index for sorting an element that might not be contained in a SORTED G ##
FindIndex:=proc(pp,GG)### e.g. pp=[[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar]
local i, F1,F2,index,g;
global Vars, tord, ord;

F1:=FHelp(pp,GG[1]):
if F1 then
	return 0;
fi;
F2:=FHelp(pp,GG[2]):
for i from 2 to nops(GG)-1 do
	if not F1 and F2 then
		return i;
	else
		F1:=F2:
		F2:=FHelp(pp,GG[i+1]):
	fi;
od:
return nops(GG);
end:

### Class for a term p (Term has coeff 1) ####
Class:=proc(p)
global Vars;
local Ind, i;
if p=0 then
	return 0;
else
	Ind:=indets(p):
	if Ind={} then
		return nops(Vars);
	fi;
	for i from 1 to nops(Vars) do
		if Vars[i] in Ind then
			return i;
		fi;
	od:
fi;
end:


########### Does q divide p wrt Pommaret-division? #########################
InvDiv:=proc(p,q) ### p,q are terms
local cp, ep, cq, eq, i;
global Vars;
cp:=Class(p):
cq:=Class(q):
if cp > cq then
	return false;
elif cq=0 then ### not used
	return false;
elif cp=0 then ## not used
        return true;
else  #### so cp <= cq and both p,q are non-zero
	ep := [seq(degree(p, Vars[i]), i=1..nops(Vars))]:
	eq := [seq(degree(q, Vars[i]), i=1..nops(Vars))]:
    if ep[cq] < eq[cq] then
   		return false;
   	else
        for i from cq+1 to nops(Vars) do 
         	if ep[i]<>eq[i] then
          		return false;
          	fi;
       	od:
    fi;
	return true;
fi;
end:



##### Does G contain a inv. divisor for p wrt Pommaret division? If yes return a list of all divisors #########
InvDivG:=proc(p,GG) ####### p is a term! Elements of GG, and qq1 are of the form [[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar]
local i, pqi;
global Vars;
pqi:=[]:
for i from 1 to nops(GG) do
    if GG[i][2][3]<>0 then 
    	if InvDiv(p,GG[i][2][2]) then
    		pqi:=[op(pqi), GG[i] ]:
    	fi;
    fi;
od:
return pqi; ### can return []
end:



#### Leading term for POT-Lift with ei<ej for i<j ####
PotLM:=proc(up) ### up is a vector of polynomials
local i,k, LM, LC;
global tord, Vars;
if up=[seq(0,k=1..nops(up))] then 
	return [];
fi;
for i from 0 to nops(up)-1 do ### comment on paper
	if up[nops(up)-i]<>0 then  
		LC,LM:=LeadingTerm(collect(up[nops(up)-i],Vars,distributed), tord):
		return [LC,LM,nops(up)-i];
	fi;
od:
end:


#### Leading term for TOP-Lift with ei<ej for i<j ####
TopLM:=proc(up) ### up is a vector of polynomials
local i,j,k,LM, LC, c,t;
global tord,Vars;
if up=[seq(0,k=1..nops(up))] then 
	return [];
fi;
LC,LM:=LeadingTerm(collect(up[1],Vars,distributed),tord):
j:=1:
for i from 2 to nops(up) do
	if up[i]<>0 then
		c,t:=LeadingTerm(collect(up[i],Vars,distributed),tord):
		if not TestOrder(t,LM,tord) then
			j:=i:
			LC,LM:=c,t:
		elif  t=LM and j<i	then
			j:=i:
			LC,LM:=c,t:
		fi;
	fi;
od:
return [LC, LM, j];
end:


####### is pp inv covered by qq? ###############
Covered:=proc(pp,qq,Syzbool)     #### e.g. pp=[[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar]
local LTup, LTuq, LTvp, LTvq;
global tord, Vars;
if pp[1][3]<>qq[1][3] then
	return false;
fi;
if Syzbool then
	if not InvDiv(pp[1][2],qq[1][2]) then
		return false;
	elif pp[2][3]<>0 and qq[2][3]<>0 then
		return not TestOrder(pp[2][2],qq[2][2]*pp[1][2]/qq[1][2],tord);
	elif pp[2][3]=0 and qq[2][3]<>0 then
		return false;
	else ## so inv covered by a syzygy
		return true;
	fi;
else 
	if not divide(pp[1][2],qq[1][2]) then
		return false;
	elif pp[2][3]<>0 and qq[2][3]<>0 then
		return not TestOrder(pp[2][2],qq[2][2]*pp[1][2]/qq[1][2],tord);
	elif pp[2][3]=0 and qq[2][3]<>0 then
		return false;
	else ## so inv covered by a syzygy
		return true;
	fi;
fi;

end:
	

###### is p inv covered by GG? ######## 
CoveredGG:=proc(pp,GG,Syzbool) #### Elements in GG, and pp are of the form [[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar]
local i;
global tord, Vars;
if GG=[] then
	return false;
fi;
for i from 1 to nops(GG) do
	if Covered(pp, GG[i],Syzbool) then
		return true;
	fi;
od:
return false;
end:


##### is pp inv super reducible by qq? #####
InvIrreg:=proc(pp,qq,Syzbool)     #### e.g. pp=[[LCup,LTup,i,up], [LCvp,LTvp,vp], JVar] 
global Vars;
if pp[1][3]<>qq[1][3] then
	return false;
fi;
if qq[1][1]=0 then
	return false:
fi;
if Syzbool then
	if not InvDiv(pp[1][2],qq[1][2]) then
		return false;
	fi;
else
	if not divide(pp[1][2],qq[1][2]) then
		return false;
	fi;
fi;
if qq[2][3]=0 then
    return true;
fi;
##### so qq[2][3]<>0
if pp[2][3]=0 then
   return false;
fi;         

  ## so here both are <>0	
if not InvDiv(pp[2][2],qq[2][2]) then
	return false;
else
	if pp[1][1]*pp[1][2]/(qq[1][1]*qq[1][2])<>pp[2][1]*pp[2][2]/(qq[2][1]*qq[2][2]) then
		return false;
	else
		return true;
	fi;
fi;
end:

#### is pp inv super reducible by GG\{pp}? #######
InvIrregGG:=proc(pp,GG,Syzbool) #### Elements in GG, and pp are of the form [[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar]
local i;
global Vars;
if GG=[] or pp=[] then   
	return [false];
fi;
if member(pp, GG, 'l') then
	if member(pp, GG minus [GG[l]], 'k') then ## l \leq k
		return [true, k+1]; ## k+1 because we took out an element.
	fi;
fi;
for i from 1 to nops(GG) do
	if [pp[1],pp[2]]<>[GG[i][1],GG[i][2]] and InvIrreg(pp, GG[i],Syzbool) then ### with minus function
		return [true,i];
	fi;
od:
return [false];
end:

															
#### is pp inv super reducible by GG where pp not to be said an element in GG
InvIrregG:=proc(pp,GG,Syzbool) #### Elements in GG, and pp are of the form [[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar]
local i;
global Vars;
if GG=[] or pp=[] then   
	return [false];
fi;
for i from 1 to nops(GG) do
	if InvIrreg(pp, GG[i],Syzbool) then 
		return [true,i];
	fi;
od:
return [false];
end:





##### is pp inv. reg reduzible by qq? (POT-Lift ord=1 or TOP-Lift ord=2) #####
InvReg:=proc(pp,qq)     #### e.g. pp=[[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar]
global tord, Vars, ord;
if [pp[1],pp[2]]=[qq[1],qq[2]] then ### actually covered by the rest of the code, but its a fast check at the beginning
	return [false,[]];
fi;
if ord=1 then
	if pp[1][3]<qq[1][3] then 
		return [false,[]];
	fi;
fi;

if pp[2][3]=0 or qq[2][3]=0 then 
    return [false,[]];
fi;
if not InvDiv(pp[2][2],qq[2][2]) then
	return [false,[]];
fi;

if ord=1 then
	if pp[1][3]>qq[1][3] then
		return [true,pp[1][1]];
	fi;
fi; ### positions are equal for ord=1

if TestOrder(qq[1][2]*pp[2][2]/qq[2][2],pp[1][2],tord) then
	if qq[1][2]*pp[2][2]/qq[2][2]=pp[1][2] then
		if ord=2 then
			if pp[1][3]>qq[1][3] then	
				return [true,pp[1][1]];
			elif pp[1][3]<qq[1][3] then
				return [false,[]];
			fi; ### now positions are equal for ord=2, too
		fi;
		if pp[1][1]/qq[1][1]=pp[2][1]/qq[2][1] then
			return [false,[]];
		else 
			return [true,pp[1][1]-pp[2][1]*qq[1][1]/qq[2][1]];
		fi;
	fi;
	return [true,pp[1][1]];
else
	return [false,[]];
fi;
end:


#### is pp inv reg red. by GG? If yes, return the first element with that property ###
InvRegG:=proc(pp,GG, index) ##### elements in GG, and pp are of the form [[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar] ######
local i,booltupel;
global tord, Vars, ord;
for i from 1 to index do
	if not i> nops(GG) then
		booltupel:=InvReg(pp,GG[i]):
		if booltupel[1] then
			return [GG[i],booltupel[2]];
		fi;
	fi;
od:
return []; ## if no reg red. is possible
end:

########### perform a inv regular red of pp by qq. ATTENTION: IT IS ASSUMED THAT THE REDUCTION IS POSSIBLE!! ############
InvRegRed:=proc(pp,qq,cc) #### e.g. pp=[[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar]
local u,v,c,t,cls,var,i;
global tord, Vars, ord, SigOnly;
if ord=1 then
	if not SigOnly then
		u:=[seq(collect(pp[1][4][i] - pp[2][1]*pp[2][2]/(qq[2][1]*qq[2][2])*qq[1][4][i], Vars, distributed),i=1..nops(pp[1][4]))]: ### up-lm(vp)/lm(vq) *uq
		v:=collect(pp[2][3]-pp[2][1]*pp[2][2]/(qq[2][1]*qq[2][2])*qq[2][3], Vars, distributed): ### vp - lm(vp)/lm(vq) *vq
		c,t:=LeadingTerm(v,tord):
		return [[op(PotLM(u)),u],[c,t,v], JVars(t)]; ### actually one does not need POTLM
	else
		#u:=[seq(collect(pp[1][4][i] - pp[2][1]*pp[2][2]/(qq[2][1]*qq[2][2])*qq[1][4][i], Vars, distributed),i=1..nops(pp[1][4]))]: ### up-lm(vp)/lm(vq) *uq
		v:=collect(pp[2][3]-pp[2][1]*pp[2][2]/(qq[2][1]*qq[2][2])*qq[2][3], Vars, distributed): ### vp - lm(vp)/lm(vq) *vq
		c,t:=LeadingTerm(v,tord):
		return [[cc,pp[1][2],pp[1][3],pp[1][4]],[c,t,v], JVars(t)];
	fi;
elif ord=2 then
	if not SigOnly then
		u:=[seq(collect(pp[1][4][i] - pp[2][1]*pp[2][2]/(qq[2][1]*qq[2][2])*qq[1][4][i], Vars, distributed),i=1..nops(pp[1][4]))]: ### up-lm(vp)/lm(vq) *uq
		v:=collect(pp[2][3]-pp[2][1]*pp[2][2]/(qq[2][1]*qq[2][2])*qq[2][3], Vars, distributed): ### vp - lm(vp)/lm(vq) *vq
		c,t:=LeadingTerm(v,tord):
		return [[op(TopLM(u)),u],[c,t,v], JVars(t)];  ### actually one does not need TOPLM
		#return [[cc,pp[1][2],pp[1][3],u],[c,t,v], JVars(t)]
	else
		#u:=[seq(collect(pp[1][4][i] - pp[2][1]*pp[2][2]/(qq[2][1]*qq[2][2])*qq[1][4][i], Vars, distributed),i=1..nops(pp[1][4]))]: ### up-lm(vp)/lm(vq) *uq
		v:=collect(pp[2][3]-pp[2][1]*pp[2][2]/(qq[2][1]*qq[2][2])*qq[2][3], Vars, distributed): ### vp - lm(vp)/lm(vq) *vq
		c,t:=LeadingTerm(v,tord):
		return [[cc,pp[1][2],pp[1][3],pp[1][4]],[c,t,v], JVars(t)];
	fi;
fi;
end:

###### calculate the inv reg normalform of pp by GG

InvRegNF:=proc(pp,GG, index) ###### elements in GG, and pp are of the form [[LCup,LTup,i,up], [LCvp,LTvp, vp],JVar] 
local qq,p,helptupel;
global tord, Vars, ord, SigOnly;
if pp[2][3]=0 then
	return pp;
fi;
#print(index);
helptupel:=InvRegG(pp,GG,index):
if helptupel <>[] then
	qq:=helptupel[1]:
fi;
p:=copy(pp):
while helptupel<>[] do  ## so reduction is possible
	p:=InvRegRed(p,qq,helptupel[2]):
	helptupel:=InvRegG(p,GG, index):
	if helptupel<>[] then
		qq:= helptupel[1]:
	fi;
od:
return p;
end:



#####################################
#########################################
################ FOR TESTBASIS ONLY ######################
#### 
Inv:=proc(pp,qq)     #### e.g. pp=[[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar]
global tord, Vars, ord;
if not InvDiv(pp[2][2],qq[2][2]) then
	return false;
fi;
return true;
end:


##########is pp inv reg red. by GG? If yes, return the first element with that property ###
InvG:=proc(pp, GG) ##### elements in GG, and pp are of the form [[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar] ######
local i;
global tord, Vars, ord;
for i from 1 to nops(GG) do
	if Inv(pp,GG[i]) then
		return GG[i];
	fi;
od:
return []; ## if no reg red. is possible
end:

InvRed:=proc(pp,qq) ### reduce pp by qq. It is assumed that this is possible!
local v,c,t,p;
global tord, Vars;
p:=copy(pp):
v:=collect(pp[2][3]-pp[2][1]*pp[2][2]/(qq[2][1]*qq[2][2])*qq[2][3], Vars, distributed): ### vp - lm(vp)/lm(vq) *vq
if v=0 then
	return 0;
fi;
c,t:=LeadingTerm(v,tord):
p:=subsop(2=[c,t,v],p):
return p;
end:


InvNF:=proc(pp,GG) ###### elements in GG, and pp are of the form [[LCup,LTup,i,up], [LCvp,LTvp, vp],JVar] 
local qq,p;
global tord, Vars, ord;
if pp[2][3]=0 then
	return 0;
fi;
qq:=InvG(pp,GG):
p:=copy(pp):
while qq<>[] do  ## so reduction is possible
	p:=InvRed(p,qq): 
	if p=0 then
		return 0;
	fi;
	qq:=InvG(p,GG):  		
od:
return p;
end:


##############
#########################
####################


######## returns the maximal degree of v-part of GG3
maxdeg:=proc(GG3)
local i, deg;
deg:=degree(GG3[1][2][2]):
for i from 2 to nops(GG3) do
	if degree(GG3[i][2][2])>deg then
		deg:=degree(GG3[i][2][2]):
	fi;
od:
return deg;
end:
	

############ Fill the set BG given by the inv. reg. autoreduced set gg, which are induced by q. For q=[] every element in gg is looked at) ) ################
InvBG:=proc(q,GG,HH,Syzbool,degbound)    ###########  q=[[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar];  
local bg, dg, bgg,pp,BG,JP3,JP, i,j,c,k,t, l,s,m, cls,var, qq,invdivg, g, div, n,p, vari;   ### pp is never a syzygy
global Vars, ord, SigOnly;
pp:=copy(q):
if GG=[] then
	return [];
fi;
if SigOnly then
	dg:=degbound-1:
else
	dg:=degbound:
fi;
JP:=[]:
bg:=[]:
BG:=[]:
JP3:=[]:
vari:=[]:
if pp=[] then
	i:=1:
	while not i > nops(GG) do
		if not degree(GG[i][2][2])>dg+3 then 
			invdivg:=InvDivG(GG[i][2][2], GG minus [ GG[i] ]): ## list of all p where p[2][2] is a inv divisor of GG[i][2][2]. same formate as pp (see above)
			if invdivg <>[] then
				for m from 1 to nops(invdivg) do
					div:=invdivg[m]:
					c:=GG[i][2][2]/div[2][2]:
					bgg:= [        [ div[1][1],collect(c*div[1][2],Vars,distributed), div[1][3], 
						[ seq( collect(c*div[1][4][j], Vars, distributed), j=1..nops(div[1][4]) )]] , 
								   [div[2][1], GG[i][2][2], collect(c*div[2][3],Vars,distributed)], JVars(GG[i][2][2])  ] :
					if Syzbool then
						if InvDiv(c*div[1][2],div[1][2]) then
							if InvReg(bgg,GG[i])[1] then
								if not CoveredGG(bgg, [op(HH),op(GG)], Syzbool) then
									BG:=[op(BG),bgg]:
								fi;
							fi;
						fi;
					else 
						if (InvReg(bgg,GG[i])[1] and not CoveredGG(bgg, [op(HH),op(GG)],Syzbool)) then
							BG:=[op(BG),bgg]:
						fi;
					fi;
				od:
			fi;
		fi;
		i:=i+1:
	od:
else 
	invdivg:=InvDivG(pp[2][2], GG): 
	if nops(invdivg) >0 then  
		for t from 1 to nops (invdivg) do
			if [pp[1],pp[2]]<> [invdivg[t][1],invdivg[t][2]] and not degree(pp[2][2])>dg+3 then ## pp  in invdivg shall not be considered
				div:= invdivg[t]: 
				#print(hiiiiier);
				c:=pp[2][2]/div[2][2]:
				bgg:= [        [ div[1][1],collect(c*div[1][2],Vars,distributed), div[1][3], 
					[ seq( collect(c*div[1][4][j], Vars, distributed), j=1..nops(div[1][4]) )]] , 
		                	   [div[2][1], pp[2][2], collect(c*div[2][3],Vars,distributed)], JVars(pp[2][2])  ] :
				if Syzbool then
					if InvDiv(c*div[1][2],div[1][2]) then
						if InvReg(bgg,pp)[1] then
							if not CoveredGG(bgg, [op(HH),op(GG)], Syzbool) then
								BG:=[op(BG),bgg]:
							fi;
						fi;
					fi;
				else 
					if (InvReg(bgg,pp)[1] and not CoveredGG(bgg, [op(HH),op(GG)],Syzbool)) then
						BG:=[op(BG),bgg]:
					fi;
				fi;
			fi;
		od:	
	fi;	
	for l from 1 to nops(GG) do
		if [GG[l][1],GG[l][2]]<>[pp[1],pp[2]] then
			if InvDiv(GG[l][2][2], pp[2][2]) and not degree(GG[l][2][2])>dg+3 then
				div:=pp:
				p:= GG[l]:
				c:=p[2][2]/div[2][2]:
				bgg:=[        [ div[1][1],collect(c*div[1][2],Vars,distributed), div[1][3], 
					[ seq( collect(c*div[1][4][j], Vars, distributed), j=1..nops(div[1][4]) )]] , 
		                	   [div[2][1], p[2][2], collect(c*div[2][3],Vars,distributed)], JVars(p[2][2])  ]:
				if Syzbool then
					if (InvDiv(c*div[1][2],div[1][2]) and  InvReg(bgg,p)[1] and not CoveredGG(bgg, [op(HH),op(GG)], Syzbool)) then
						BG:=[op(BG),bgg]:
					fi;
				else 
					if (InvReg(bgg,p)[1] and not CoveredGG(bgg, [op(HH),op(GG)],Syzbool)) then
						BG:=[op(BG),bgg]:
					fi;
				fi;
			fi;
		fi;
	od:
fi;
return BG;
end:


###### perform a inv reg head autoreduction of GG , collect all obtained syzygies #########

InvRegAuto:=proc(GG) ######## Elements in GG are of the form [[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar], INEFFICIENT CODE!!
local pp, rr,t, gg, HH,s, l, JP,jp,Q, i;           
global tord, Vars, SigOnly;								
gg:=copy(GG):
if not nops(gg) >1 then
	return [gg,[],[]];
fi;
t:=1:
HH:=[]:
while not t >  nops(gg) do ## combine it with collecting sygygies
	pp:=gg[t]:
	#debug(InvRegNF);
   	rr:= InvRegNF(pp,gg,t):
	if rr<>pp then
        if rr[2][3]<>0 then
			gg:=sort(subsop(t=rr,gg),FHelp):
		else
			HH:=[op(HH),rr]:
			gg:=subsop(t=NULL,gg):
        fi;
	    t:=1:
	else
		t:=t+1:	
	fi;
od:
return [gg, HH];

end:

	
###### Create GG out of G=[[u1,v1],...,[un,vn]] #######
SetGG:=proc(G)
global tord, Vars, ord;
local GG, i,j,c,t,cls,var;

GG:=[]:
if ord=1 then
	if not nops(G)<1 then
		c,t:=LeadingTerm(G[1][2],tord):
		var:=JVars(t):
		GG:=[ [[op(PotLM(G[1][1])),G[1][1]],[c,t,G[1][2]] ,var ] ]:
	else
		return [];
	fi;
	for i from 2 to nops(G) do
		c,t:=LeadingTerm(G[i][2],tord):
		var:=JVars(t):
		GG:=[op(GG),    [[op(PotLM(G[i][1])),G[i][1]],[c,t,G[i][2]],var ]  ]:
	od:
elif ord=2 then
	if not nops(G)<1 then
		c,t:=LeadingTerm(G[1][2],tord):
		var:=JVars(t):
		GG:=[ [[op(TopLM(G[1][1])),G[1][1]],[c,t,G[1][2]] ,var ] ]:
	else
		return [];
	fi;
	for i from 2 to nops(G) do
		c,t:=LeadingTerm(G[i][2],tord):
		var:=JVars(t):
		GG:=[op(GG),    [[op(TopLM(G[i][1])),G[i][1]],[c,t,G[i][2]],var ]  ]:
	od:
fi;
return collect(GG,Vars,distributed);
end:

### calculate JVar for some term p ############
JVars:=proc(p) #### p is a term
local i,cls,jvar;
global Vars;
if p=0 then	#### should never be the case...
	return []:
fi;
jvar:=[]:
cls:=Class(p):
if cls<nops(Vars) then
	for i from cls+1 to nops(Vars) do
		jvar:=[op(jvar), Vars[i]]:
	od:
fi;
return jvar;
end:
		

#### calculate the smallest involutive J-Pair of pp  wrt x1<x2<... if Vars=[x1,x2,...] #############
InvJP:=proc(pp,q,Syzbool) ### pp:=[[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar] 
local c, i, j, JP, jpp;
global Vars;
if pp=[] or pp[2][3]=0 or degree(pp[2][2])>q+2 then	
	return [];
fi;
if nops(pp[3])=0 then
	return [];
else
	JP:=[]:
	jpp:=[ [ pp[1][1], collect(pp[3][1]* pp[1][2],Vars,distributed), pp[1][3], [seq(collect(pp[3][1]*pp[1][4][j],Vars,distributed), j=1..nops(pp[1][4]))] ]   ,[pp[2][1], seq(collect(pp[3][1]* pp[2][j],Vars,distributed),j=2..3)], JVars(pp[2][2]) ]:
fi;
return [jpp];
end:

####### Obtain the correct format for the input out of F=[v1,...,vn] ##########
G0:=proc(F)
local i, j, g0, e;
global tord,Vars, ord;
if F=[] then
	return [];
fi;
e:=[seq(0, j=1..nops(F))]:
g0:=[]:
for i from 1 to nops(F) do	
	g0:=[op(g0),[subsop(i=1,e), F[i]] ]:
od:
return g0;
end:

###############################


###### Transform a syzygy into the right formate  [[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar] #######
SyzHH:=proc(H) ######### Elements in H are lists of polynomials
local h, i,j,k,e,h0,t;
global tord, Vars;
h:=[]:
for i from 1 to nops(H) do
	h:=[op(h), [H[i],0] ]:
od:
return SetGG(h);
end:
	
#### perform an inv. super reduction of pp  by qq where both are syzygies- ATTENTION: IT IS ASSUMED THAT THE REDUCTION IS POSSIBLE #####
InvIrregRed:=proc(pp,qq) #### e.g. pp=[[LCup,LTup,i,up], [0,1, 0], JVar]
local u,i,j;
global tord, Vars, ord;
u:=[seq(collect(pp[1][4][i] - pp[1][1]*pp[1][2]/(qq[1][1]*qq[1][2])*qq[1][4][i], Vars, distributed),i=1..nops(pp[1][4]))]:
if u=[seq(0, j=1..nops(u))] then	
	return [];
fi;
if ord=1 then
	return [[op(PotLM(u)),u],[0,1,0], []];
elif ord=2 then
	return [[op(TopLM(u)),u],[0,1,0], []];
fi;
end:

### calculate the inv. normalform of the syzygy pp wrt to the set HH of all syzygies obtained so far #####
InvIrregNF:=proc(pp,HH,Syzbool) ##### elements in HH, and pp are of the form [[LCup,LTup,i,up], [0,1,0], JVar] ######
local qq,p;
global tord, Vars, ord;
qq:=InvIrregGG(pp,HH,Syzbool): ### qq =[false] or [true, Position of element that reduces pp]
p:=copy(pp):
while qq[1] do
	p:=InvIrregRed(p,HH[qq[2]]):
	qq:=InvIrregGG(p,HH,Syzbool):
od:
return p;
end:

#### Autoreduziere Menge HH von Syzygien 

InvIrregAuto:=proc(HH,Syzbool)		###BAUSTELLE
local pp, rr, k,l, m,t, hh,s, Num;
global tord, Vars, ord;
hh:=copy(HH):
if not nops(hh) >1 then
	return hh;
fi;
t:=1:
Num:= nops(hh):
if Num = 1 then
	return hh;
fi;
while not t >  Num do
	pp:=hh[t]:
	rr:=InvIrregNF(pp,hh,Syzbool):
	if rr<>pp then
		if rr<>[] then
			hh:=subsop(t=rr,hh):
       	else
			hh:= subsop(t=NULL,hh):
			Num:= Num-1:
		fi;	
	t:=1:
	else
		t:=t+1:	
	fi;
od:
return hh;
end:
	

### transform GG/ HH to G /H ###

SetG:=proc(GG)
local G, i;
G:=[]:
for i from 1 to nops(GG) do
	G:=[op(G), [GG[i][1][4],GG[i][2][3]] ]:
od:
return G:
end:

### is pp pseudo reducible by GG?
PseudoRed:=proc(pp,GG,Syzbool)
local i;
global tord, Wars, ord;
for i from 1 to nops(GG) do
	if [pp[1],pp[2]]<>[GG[i][1],GG[i][2]] then
		if pp[1][3]<>GG[i][3] then
			return false;
		fi;
		if GG[i][1]=0 then
			return false:
		fi;
		if Syzbool then
			if not InvDiv(pp[1][2],GG[i][1][2]) then
				return false;
			fi;
		else
			if not divide(pp[1][2],GG[i][1][2]) then
				return false;
			fi;
		fi;
		if GG[i][2][3]=0 then
			return false;
		fi;      

		  ## so here both are <>0	
		if not divide(pp[2][2],GG[i][2][2]) then
			return false;
		else
			if pp[1][2]/(GG[i][1][2])<>pp[2][2]/(GG[i][2][2]) then
				return false;
			else
				return true;
			fi;
		fi;	
	fi;
od:
return false;	
end:

############################ TEST THE OUTPUT

TestBasis:=proc(GH,Syzbool,q) ###elements in GH=GG\cup HH are of the form [[LCup,LTup,i,up], [LCvp,LTvp, vp], JVar]
local GG, HH, JPG,JPH, i,j,m,rr,pp, k,l,hh,qq,gg,cls,a; ### inefficient reduction algorithm
global tord, Vars, ord, SigOnly;
GG:=GH[1]:
JPG:=[]:
HH:=GH[2]:
JPH:=[]:
if Syzbool and not SigOnly then
	for k from 1 to nops(HH) do	
		HH:=subsop(k=subsop(2=subsop(1=HH[k][1][1], HH[k][2]),HH[k]),HH):###HH[k][2][1]:=HH[k][1][1]:
		HH:=subsop(k=subsop(2=subsop(2=HH[k][1][2], HH[k][2]),HH[k]),HH):#HH[k][2][2]:=HH[k][1][2]:
		HH:=subsop(k=subsop(2=subsop(3=HH[k][1][4][HH[k][1][3]], HH[k][2]),HH[k]),HH):#HH[k][2][3]:=HH[k][1][4][HH[k][1][3]]:		
		cls:=Class(HH[k][2][2]):
		if cls<nops(Vars) then
			for i from cls+1 to nops(Vars) do
				HH:=subsop(k=subsop(3=i,HH[k]),HH):				#HH[k][3]:=[i]:
				JPH:=[op(JPH),op(InvJP(HH[k],q,Syzbool))]:
				HH:=subsop(k=subsop(3=[],HH[k]),HH):	
			od:
		fi;
		HH:=subsop(k=subsop(2=subsop(1=0, HH[k][2]),HH[k]),HH):###HH[k][2][1]:=0
		HH:=subsop(k=subsop(2=subsop(2=1, HH[k][2]),HH[k]),HH):#HH[k][2][2]:=1
		HH:=subsop(k=subsop(2=subsop(3=0, HH[k][2]),HH[k]),HH):#HH[k][2][3]:=0	
	od:
	for l from 1 to nops(JPH) do
		JPH:=subsop(l=subsop(2=subsop(1=0, JPH[l][2]),JPH[l]),JPH):###JPH[l][2][1]:=0
		JPH:=subsop(l=subsop(2=subsop(2=1, JPH[l][2]),JPH[l]),JPH):#JPH[l][2][2]:=1
		JPH:=subsop(l=subsop(2=subsop(3=0, JPH[l][2]),JPH[l]),JPH):#JPH[l][2][3]:=0	
		for j from 1 to nops(JPH) do
			hh:=JPH[1]:
			while not member(hh,HH) and hh<>[] do
				hh:=InvIrregNF(hh,HH,Syzbool):
				if hh<>[] then
					return print("Output not a weak Pommaret basis of Syz(F)");
				fi;
			od:
		od:
	od:
	print("True for syzygy module"); 
fi;
for i from 1 to nops(GG) do
	if nops(GG[i][3])>0 then
		for m from 1 to nops(GG[i][3]) do
			gg:=[[0,1,1,[seq(0, j=1..nops(GG[i][1][4]))]],GG[i][2],[GG[i][3][m]]]:
			JPG:=[op(JPG), op(InvJP(gg,q,Syzbool))]:
		od:
	fi;
od:
for j from 1 to nops(JPG) do
	rr:=InvNF(JPG[j],GG):
	if rr<>0 then
		return print("Output not a Pommaret basis of I");
	fi;
od:
print("Output contains a Pommaret-basis for the ideal");
end:

DegFSig:=proc(GG) ####find the element in GG with min DegF wrt sygnature and return the degree and first elem with that degree
local degf, deg, i,pp,qq;
global ord, tord, Vars;
pp:=GG[1]:
degf:=degree(pp[1][2]):
for i from 2 to nops(GG) do
	 qq:=GG[i]:
	 deg:=degree(qq[1][2]):
	 if deg<degf then
		degf:=deg:
		pp:=qq:
	fi;
od:
return [pp, degf];
end:

DegFV:=proc(GG) ####find the element in GG with min DegF wrt v-part and return the degree and first elem with that degree
local degf, deg, i,pp,qq;
global ord, tord, Vars;
pp:=GG[1]:
degf:=degree(pp[2][2]):
for i from 2 to nops(GG) do
	 qq:=GG[i]:
	 deg:=degree(qq[2][2]):
	 if deg<degf then
		degf:=deg:
		pp:=qq:
	fi;
od:
if degf<0 then
	degf:=-1:
fi;
return [pp, degf];
end:

findtrans:=proc(rr,JPP,BGG,GG, HH,string)
local T,j,i, cls, cls2, continue,k,q, trans, bigdiv;
global tord, Vars, ord;
T:=[rr,op(JPP),op(BGG)]:
trans:=[]:
if string="Syz" then
	for j from 1 to nops(T) do ######## 
		cls:=Class(T[j][1][2]):
		for i from 1 to nops(HH) do
			if divide(T[j][1][2], HH[i][1][2],'q') and not InvIrregGG(T[j][1][2], HH,true)[1] then ### if it is ordinary divisible but not inv.
				cls2:= Class(HH[i][1][2]): 				#####then take non-mult vars of HH[i][1][2] (which corresponds to some indets(q)[k]) 
				continue:=true:
				k:=1:
				while continue and not k>nops(indets(q)) do
					if Class(indets(q)[k])>cls2 and not member([Vars[cls2],indets(q)[k]],trans) then
						trans:=[op(trans),[Vars[cls2],indets(q)[k]]]:
#						continue:=false:				# taking only the first/smallest obstruction
					fi;
					k:=k+1:
				od:
			fi;
		od:
	od:
	return trans;
elif string="Ideal" then 
	for j from 1 to nops(T) do
		cls:=Class(T[j][2][2]):
		for i from 1 to nops(GG) do
			bigdiv:=BigDiv(T[j],GG):
			if bigdiv<>[] then
				q:=T[j][2][2]/bigdiv:
				if InvDivG(T[j][2][2], GG)=[] then ### if it is ordinary divisible but not inv. 
					##cls2:= Class(GG[i][2][2]): 				#####then take non-mult vars of GG[i][1][2] (which corresponds to some indets(q)[k]) 
					cls2:=Class(bigdiv):
					continue:=true:
					k:=1:
					while continue and not k>nops(indets(q)) do
						if Class(indets(q)[k])>cls2 and not member([Vars[cls2],indets(q)[k]],trans) then
							trans:=[op(trans),[Vars[cls2],indets(q)[k]]]:
#							continue:=false:				# taking only the first obstruction
						fi;
						k:=k+1:
					od:
				fi;
			fi;
		od:
	od:
	return trans;
fi;
end:


IndexOfSafety:=proc(gg,HH,q,Syzbool) ## only counting superfluous inv. J-Pairs!
local jpp, i, j, GG, minjp, Out,Q, ind;
jpp:=[]:

GG:=copy(gg):
for i from 1 to nops(GG) do
	if nops(GG[i][3])>0 then
		for j from 1 to nops(GG[i][3]) do
			jpp:=[op(jpp),op(InvJP(GG[i],q,Syzbool))]:
			GG:=subsop(i=subsop(3=subsop(1=NULL,GG[i][3]),GG[i]),GG); ######GG[i][3]:=subsop(1=NULL, GG[i][3]):
		od:
	fi;
	#for j from 1 to nops(jpp) do
	#	jpp:=subsop(j=subsop(3=[],jpp[j]),jpp):
	#od:
od:
L:=sort(InvBG([],gg,HH,Syzbool,q+2),FHelp):

minjp:=1:
Out:=false:
for j from 1 to nops(jpp) do ## ALLE KRITERIEN HIER HIN!!
	if not Out then
		if not(CoveredGG(jpp[j],sort([op(HH), op(GG)],FHelp),Syzbool)) then  
			#if jpp[j][3]=[] then
				if not(InvIrregG(jpp[j],sort([op(HH), op(GG)],FHelp),Syzbool)[1]) then ## Vorher im "Not" : jpp[j][3]=1 and 
					minjp:= j;
					Out:=true:
				fi;
			#else
			#	return j;
			#fi;
		fi;
	fi;
od:
if nops(jpp)+nops(L)>0 then
	Q:=sort([op(jpp),op(L)],FHelp):
fi;


ind:=FindIndex(jpp[minjp],Q):

if ind<nops(Q) then
	return [ind,minjp];
else 
	return [[]]; ### if all are superfluous
fi;
end:

################ fill JP for the POT lift with all J-pairs of the next element in a sorted GG
FilljpPot:=proc(GG,q,Syzbool)
local i,j, JPP,gg;
gg:=copy(GG):
JPP:=[]:
for i from 1 to nops(GG) do
	if nops(GG[i][3])>0 then
		for j from 1 to nops(GG[i][3]) do
			JPP:=[op(JPP),op(InvJP(gg[i],q,Syzbool))]:
			gg:=subsop(i=subsop(3=subsop(1=NULL,gg[i][3]),gg[i]),gg); ######gg[i][3]:=subsop(1=NULL, gg[i][3]):
		od:
	fi;
	if JPP<>[] then
		return [JPP,gg];
	fi;
od:
return [JPP,gg];
end:


BigDiv:=proc(pp,GG)
local i,j,g,p;
global Vars, tord;
p:=[]:
for i from 1 to nops(GG) do
	if divide(pp[2][2],GG[i][2][2]) then	
		p:=[op(p),GG[i][2][2]]:
	fi;
od:
if nops(p)=1 then
	return p[1];
fi;
if p=[] then
	return [];
fi;
g:=p[1]:
for j from 2 to nops(p) do
	if Class(g[2][2])<Class(p[j][2][2]) then
	#if Class(g)<Class(p[j]) then
		g:=p[j]:
	fi;
od:
return g;
end:
	




###############################################
######## FINAL ALGORITHM ################ 
StrPBas:=proc(F,H,Syzbool,q,RestartBool) #### + POT-Lift (ord=1) or Top-Lift (ord=2) of tord with ei<ej for i<j
local G2,G1,JG,H2,m,maxdg,superfl, GH,vari, GG, BGG, JPP, HH,G,H0,Tran,safe,transk,s, firstbytes, secondbytes, secondtime,das,help,minsigtupel, minsig,j,k, minJPtupel, minjp, minBGtupel, minbg,min, rr, l1,l2,l11, l22, hilfe, ende, balance,CountRed, CountCrit, i, bg,JP,g, index;
global tord, vars, ord, firsttime,CountTran, SigOnly,maxsafe, Superfl;

if not RestartBool then
	firsttime, firstbytes := kernelopts(cputime, bytesused);
	GG:=SetGG(G0(F)):
	GH:=InvRegAuto(collect(GG,Vars,distributed)):			
	GG:=sort(GH[1],FHelp):
	HH:=[op(GH[2]),op(SyzHH(H))]:
	maxsafe:=0:
	Superfl:=0:
else
	#print(F);
	#GH:=InvRegAuto(collect(F,Vars,distributed)):			## Here F and H are already in right formate
	#if nops(GH[1])>1 then
	#	GG:=sort(GH[1],FHelp):
	#	print("");
	#else
	#	GG:=GH[1]:
	#fi;
	#print(GG);
	#HH:=[op(GH[2]),op(H)]:
	#HH:=InvIrregAuto([op(GH[2]),op(H)],Syzbool):
	GG:=sort(F,FHelp):
	HH:=H:
fi;
CountRed:=0:
CountCrit:=0:	


#HH:=InvIrregAuto(sort([op(GH[2]),op(SyzHH(H))],FHelp),Syzbool):

JPP:=[]:
BGG:=[]:
if maxsafe<>[] then
	BGG:=InvBG([],GG,HH,Syzbool,q):
	#if ord=2 then
		for i from 1 to nops(GG) do
			if nops(GG[i][3])>0 then
				for m from 1 to nops(GG[i][3]) do 
					JPP:=[op(JPP),op(InvJP(GG[i],q,Syzbool))]:
					GG:=subsop(i=subsop(3=subsop(1=NULL,GG[i][3]),GG[i]),GG); ######GG[i][3]:=subsop(1=NULL, GG[i][3]):
				od:
			fi;
		od:
	#	if ord=1 then
	#		JG:=FilljpPot(GG,q,Syzbool):
	#		JPP:=[op(JPP),op(JG[1])]:
	#		GG:=JG[2]:
	#	fi;
	#fi;
	JPP:=sort(JPP,FHelp):
	BGG:=sort(BGG,FHelp):
	for i from 1 to Superfl -1 do		### here we are removing the elements before index maxsafe(=indexofsafety)
	#print(maxsafe), print(nops(JPP)+nops(BGG))	;
#		if nops(JPP)>0 and nops(BGG)>0 then
#			superfl:=sort([JPP[1],BGG[1]],FHelp):
#			if member(superfl,JPP) then
				JPP:=subsop(1=NULL,JPP):
				CountCrit:=CountCrit+1:
#			else
#				BGG:=subsop(1=NULL,BGG):
#			fi;
#		fi;
#		if nops(JPP)=0 and nops(BGG)>0 then
#			BGG:=subsop(1=NULL,BGG):
#		elif nops(JPP)>0 and nops(BGG)=0 then
#			JPP:=subsop(1=NULL,JPP):
#		fi;			
	od:
fi;
l1:=nops(JPP):
l2:=nops(BGG):
minjp:=[]:
minbg:=[]:
ende:=1:
balance:=false:             
l11:=true:
l22:=true:
while l1 + l2 <> 0 do  
	if l1<>0 and l11 then    ### find candiate for the next reduction from JP. dont enter if u already have a candidate (l11=false)
		#minJPtupel:=MinSig(JPP,false):
        #	minjp:=minJPtupel[1]:
		minjp:=JPP[1]:
	elif l1=0 then    
		minjp:=[]:
	fi;
    
	if l2<>0 and l22 then       ### find candiate for the next reduction from BG. dont enter if u already have a candidate (l22=false)
	#	minBGtupel:=MinSig(BGG,false):
    #	minbg:=minBGtupel[1]:
		minbg:=BGG[1]:
	elif l2=0 then
		minbg:=[]:
	fi;
	if minjp=[] and minbg<>[] then  
            min:=minbg:
           # BGG:=subsop(minBGtupel[2]= NULL, BGG):
			BGG:=subsop(1= NULL, BGG):
            l2:=l2-1:
            hilfe:=2:
            l11:=true: ### since we have no candidate for minJP
            l22:=true:
    elif minjp<>[] and minbg=[] then
       	    min:=minjp:
		  	#JPP:=subsop(minJPtupel[2]= NULL, JPP):
			JPP:=subsop(1= NULL, JPP):
            l1:= l1-1:
        	hilfe:=1:
        	l22:=true: ### since we have no candidate for minBG
            l11:= true:
    fi;
	if minjp<>[] and minbg<>[] then             #### choose the one with smaller signature       
		if ord=1 then
			if minbg[1][3]<minjp[1][3] or (minbg[1][3]= minjp[1][3] and not TestOrder(minjp[1][2], minbg[1][2], tord) ) then
				min:= minbg:
				#BGG:=subsop(minBGtupel[2]= NULL, BGG):
				BGG:=subsop(1= NULL, BGG):
				l2:=l2-1:
				hilfe:= 2:
				#l11:=false:
				l22:=true:
			else
				min:= minjp:
				#JPP:=subsop(minJPtupel[2]= NULL, JPP):
				JPP:=subsop(1= NULL, JPP):
				l1:=l1-1:
				hilfe:=1:
				if balance then ### if no J-Pair was reduced, we cant be sure, that we follow the strategy with the smallest signature
					#l22:=false:   #### so we are searching again for minBG after new elements have added to BG. 
				else
					l22:=true:
				fi;
				l11:=true:
			fi;
		fi;
		if ord=2 then
			if not TestOrder(minjp[1][2], minbg[1][2], tord) or (minjp[1][2]= minbg[1][2] and minbg[1][3]<minjp[1][3]) then
				min:= minbg:
				#BGG:=subsop(minBGtupel[2]= NULL, BGG):
				BGG:=subsop(1= NULL, BGG):
				l2:=l2-1:
				hilfe:= 2:
				#l11:=false:
				l22:=true:
			else
				min:= minjp:
				#JPP:=subsop(minJPtupel[2]= NULL, JPP):
				JPP:=subsop(1= NULL, JPP):
				l1:=l1-1:
				hilfe:=1:
				if balance then ### if no J-Pair was reduced, we cant be sure, that we follow the strategy with the smallest signature
					#l22:=false:   #### so we are searching again for minBG after new elements have added to BG. 
				else
					l22:=true:
				fi;
				l11:=true:
			fi;
		fi;
	fi;
	if min<>[] then
	#	print(min[1][2],min[1][3]);
		if not(CoveredGG(min,[op(HH), op(BGG), op(GG), op(JPP)],Syzbool) or PseudoRed(min,BGG)) then      #### if its covered or pseudo reducible it can be discarded
		   if not(hilfe=2 and InvIrregGG(min, BGG minus [min],Syzbool)[1] )then                   #### if its from BG and inv. irreg. reducible by BG\{minjp} it can be discarded
				if not(hilfe=1 and InvIrregGG(min , [op(HH), op(JPP), op(BGG), op(GG)] minus [min],Syzbool)[1])then ### if its from JP and its inv super reducible it can be discarded
					if hilfe=1 and not balance then	
						balance:= true:
					fi;
					index:= FindIndex(min,GG);
                    rr:=InvRegNF(min,GG, index):  
				#	print(InvDivG(rr[2][2], GG));
					if rr[2]<>min[2] then
						CountRed:=CountRed +1:
					fi;
			        if rr[2][3]=0 then
						HH:=[op(HH), rr]:
						#HH:=InvIrregAuto([op(HH), rr],Syzbool):
						if Syzbool and not SigOnly and DegFSig([rr,op(JPP),op(BGG)])[2]>q+1 then
							print("Syzygy module not quasi stable. Searching for index of safety");
							CountTran:=CountTran+1:
							#j:=1:
							#while not j > nops(GG) do
							#	if InvDivG(GG[j][2][2],GG minus [GG[j]])<>[] then
							#		GG:=subsop(j=NULL,GG):
							#		j:=j-1:
							#	fi;
							#	j:=j+1:
							#od:
							printf("%-1s %1s %1s: %3a\n", Num,of,syzygies,nops(HH));
							printf("%-1s %1s %1s %1s: %3a\n", Num,of,ideal, elements,nops(GG));
							printf("%-1s %1s %1s %1s %1s: %3g\n", Num,of,regular, normal, forms,CountRed);
							printf("%-1s %1s %1s %1s: %3a\n", Num,of,discarded, elements,CountCrit);
							Tran:=findtrans(rr,JPP,BGG,GG, HH,"Syz");
							safe:=[]:
			###				for k from 1 to 1 do #### saving time -- just take the first obstruction
							for k from 1 to nops(Tran) do
								G2:=subs(Tran[k][1]=Tran[k][1]+Tran[k][2],GG):
								#GH2:=InvRegAuto(subs(Tran[k][1]=Tran[k][1]+Tran[k][2],GG)):
								#G2:=GH2[1]:
								#H2:=GH2[2]:
								G1:=[]:
								for i from 1 to nops(G2) do
									G1:=[op(G1),[collect(G2[i][1][4],Vars,distributed),collect(G2[i][2][3],Vars,distributed)]]:
								od:
								GH2:=InvRegAuto(SetGG(G1)):
								G1:=GH2[1]:
								H0:=[]:
								if HH<>[] then
									H2:=subs(Tran[k][1]=Tran[k][1]+Tran[k][2],HH):
									#H2:=InvIrregAuto([op(H2),op(subs(Tran[k][1]=Tran[k][1]+Tran[k][2],HH))],Syzbool):
									for i from 1 to nops(H2) do
										H0:=[op(H0), collect(H2[i][1][4], Vars, distributed)]:
									od:
									H0:=InvIrregAuto([op(SyzHH(H0)),op(GH2[2])], Syzbool):
									#H0:=[op(SyzHH(H0)),op(GH2[2])]:
								fi;
								safe:=[op(safe), [IndexOfSafety(G1,H0,q,Syzbool),k]]:
								#print(safe);
							od:
							maxsafe:=safe[1][1][1]:
							Superfl:=safe[1][1][2]:
							transk:=safe[1][2]:
							if maxsafe=[] then
								print("Performing transformation:"),print(Tran[transk][1]),print(Tran[transk][1]+Tran[transk][2]);
								print("Index of safety:"),print(maxsafe);
								G2:=subs(Tran[transk][1]=Tran[transk][1]+Tran[transk][2],GG):
							#	GH:=InvRegAuto(subs(Tran[transk][1]=Tran[transk][1]+Tran[transk][2],GG)):
							#	G2:=GH[1]:
							#	H2:=GH2[2]:
								H2:=subs(Tran[transk][1]=Tran[transk][1]+Tran[transk][2],HH):
								#H2:=InvIrregAuto([op(H2),op(subs(Tran[k][1]=Tran[k][1]+Tran[k][2],HH))],Syzbool):
								G1:=[]:
								for i from 1 to nops(GG) do
									G1:=[op(G1),[collect(G2[i][1][4],Vars,distributed),collect(G2[i][2][3],Vars,distributed)]]:
								od:
								H0:=[]:
								for i from 1 to nops(HH) do
									H0:=[op(H0), collect(H2[i][1][4], Vars, distributed)]:
								od:
								GH2:=InvRegAuto(SetGG(G1)):
								G1:=GH2[1]:
								H0:=InvIrregAuto([op(SyzHH(H0)),op(GH2[2])], Syzbool):
								return StrPBas(G1,H0,Syzbool,q, true);
							else
								for s from 2 to nops(safe) do
									if maxsafe<safe[s][1][1] or safe[s][1][1]=[] then
										if safe[s][1][1]=[] then
											maxsafe:=safe[s][1][1]:
											print("Performing transformation:"),print(Tran[safe[s][2]][1]),print(Tran[safe[s][2]][1]+Tran[safe[s][2]][2]);
											print("Index of safety:"),print(maxsafe);
											G2:=subs(Tran[safe[s][2]][1]=Tran[safe[s][2]][1]+Tran[safe[s][2]][2],GG):
											H2:=subs(Tran[safe[s][2]][1]=Tran[safe[s][2]][1]+Tran[safe[s][2]][2],HH):
											G1:=[]:
											for i from 1 to nops(GG) do
												G1:=[op(G1),[collect(G2[i][1][4],Vars,distributed),collect(G2[i][2][3],Vars,distributed)]]:
											od:
											H0:=[]:
											for i from 1 to nops(HH) do
												H0:=[op(H0), collect(HH[i][1][4], Vars, distributed)]:
											od:
											GH2:=InvRegAuto(SetGG(G1)):
											G1:=GH2[1]:
											H0:=InvIrregAuto([op(SyzHH(H0)),op(GH2[2])], Syzbool):
											return StrPBas(G1,H0,Syzbool,q, true);
										fi;
										maxsafe:=safe[s][1][1]:
										Superfl:=safe[s][1][2]:
										transk:=safe[s][2]:
									fi;
								od:
								print("Performing transformation:"),print(Tran[transk][1]),print(Tran[transk][1]+Tran[transk][2]);
								print("Index of safety:"),print(maxsafe);
								G2:=subs(Tran[transk][1]=Tran[transk][1]+Tran[transk][2],GG):
								H2:=subs(Tran[transk][1]=Tran[transk][1]+Tran[transk][2],HH):
							#	GH:=InvRegAuto(subs(Tran[transk][1]=Tran[transk][1]+Tran[transk][2],GG)):
							#	G2:=GH[1]:
							#	H2:=GH2[2]:
							#	H2:=InvIrregAuto([op(H2),op(subs(Tran[k][1]=Tran[k][1]+Tran[k][2],HH))],Syzbool):
								G1:=[]:
								for i from 1 to nops(GG) do
									G1:=[op(G1),[collect(G2[i][1][4],Vars,distributed),collect(G2[i][2][3],Vars,distributed)]]:
								od:
								H0:=[]:
								for i from 1 to nops(HH) do
									H0:=[op(H0), collect(HH[i][1][4], Vars, distributed)]:
								od:
								GH2:=InvRegAuto(SetGG(G1)):
								G1:=GH2[1]:
								H0:=InvIrregAuto([op(SyzHH(H0)),op(GH2[2])], Syzbool):
								return StrPBas(G1,H0,Syzbool,q, true);	
							fi;
						fi;								
					else
						#print(DegFV([rr,op(JPP),op(BGG)]));
						if DegFV([rr,op(JPP),op(BGG)])[2]>q and InvDivG(rr[2][2], GG)=[] then
							print("Ideal not quasi stable. Searching for index of safety");
							CountTran:=CountTran+1:
						#	j:=1:
						##	print(1);
							#while not j > nops(GG) do
							#	if InvDivG(GG[j][2][2],GG minus [GG[j]])<>[] then
							#		GG:=subsop(j=NULL,GG):
							#		j:=j-1:
							#	fi;
							#	j:=j+1:
							#od:
							printf("%-1s %1s %1s: %3a\n", Num,of,syzygies,nops(HH));
							printf("%-1s %1s %1s %1s: %3a\n", Num,of,ideal, elements,nops(GG));
							printf("%-1s %1s %1s %1s %1s: %3g\n", Num,of,regular, normal, forms,CountRed);
							printf("%-1s %1s %1s %1s: %3a\n", Num,of,discarded, elements,CountCrit);
#print(2);
							Tran:=findtrans(rr,JPP,BGG,GG, HH,"Ideal");
							safe:=[]:
			###				for k from 1 to 1 do #### saving time -- just take the first obstruction
							for k from 1 to nops(Tran) do
								G2:=subs(Tran[k][1]=Tran[k][1]+Tran[k][2],GG):
								#GH2:=InvRegAuto(subs(Tran[k][1]=Tran[k][1]+Tran[k][2],GG)):
								#G2:=GH2[1]:
								#H2:=GH2[2]:
								G1:=[]:
								for i from 1 to nops(G2) do
									G1:=[op(G1),[collect(G2[i][1][4],Vars,distributed),collect(G2[i][2][3],Vars,distributed)]]:
								od:
								GH2:=InvRegAuto(SetGG(G1)):
								G1:=GH2[1]:
								#print(GH[2]);
								H0:=[]:
								if HH<>[] and not SigOnly then
									H2:=subs(Tran[k][1]=Tran[k][1]+Tran[k][2],HH):
									#H2:=InvIrregAuto([op(H2),op(subs(Tran[k][1]=Tran[k][1]+Tran[k][2],HH))],Syzbool):
									for i from 1 to nops(H2) do
										H0:=[op(H0), collect(H2[i][1][4], Vars, distributed)]:
									od:
									H0:=InvIrregAuto([op(SyzHH(H0)),op(GH2[2])], Syzbool):
									#H0:=[op(SyzHH(H0)),op(GH2[2])]:
								fi;
								#print(G1);
								safe:=[op(safe), [IndexOfSafety(G1,H0,q,Syzbool),k]]:
								##print(safe);
							od:
							maxsafe:=safe[1][1][1]:
							Superfl:=safe[1][1][2]:
							transk:=safe[1][2]:
							if maxsafe=[] then
								print("Performing transformation:"),print(Tran[safe[1][2]][1]),print(Tran[safe[1][2]][1]+Tran[safe[1][2]][2]);
								print("Index of safety:"),print(maxsafe);
								G2:=subs(Tran[safe[1][2]][1]=Tran[safe[1][2]][1]+Tran[safe[1][2]][2],GG):
								G1:=[]:
								for i from 1 to nops(GG) do
									G1:=[op(G1),[collect(G2[i][1][4],Vars,distributed),collect(G2[i][2][3],Vars,distributed)]]:
								od:
								GH2:=InvRegAuto(SetGG(G1)):
								G1:=sort(GH2[1],FHelp):
								H0:=[]:
								if not SigOnly then
									H2:=subs(Tran[safe[1][2]][1]=Tran[safe[1][2]][1]+Tran[safe[1][2]][2],HH):
									for i from 1 to nops(HH) do
										H0:=[op(H0), collect(H2[i][1][4], Vars, distributed)]:
									od:
									H0:=InvIrregAuto([op(SyzHH(H0)),op(GH2[2])], Syzbool):
								#	H0:=[op(SyzHH(H0)),op(GH2[2])]:
								fi;
								return StrPBas(G1,H0,Syzbool,q, true);
							else
								for s from 2 to nops(safe) do
									if maxsafe<safe[s][1][1] or safe[s][1][1]=[] then
										if safe[s][1][1]=[] then
											maxsafe:=safe[s][1][1]:
											print("Performing transformation:"),print(Tran[safe[s][2]][1]),print(Tran[safe[s][2]][1]+Tran[safe[s][2]][2]);
											print("Index of safety:"),print(maxsafe);
											G2:=subs(Tran[safe[1][2]][1]=Tran[safe[1][2]][1]+Tran[safe[1][2]][2],GG):
											G1:=[]:
											for i from 1 to nops(GG) do
												G1:=[op(G1),[collect(G2[i][1][4],Vars,distributed),collect(G2[i][2][3],Vars,distributed)]]:
											od:
											GH2:=InvRegAuto(SetGG(G1)):
											G1:=GH2[1]:
											H0:=[]:
											if not SigOnly then
												H2:=subs(Tran[safe[1][2]][1]=Tran[safe[1][2]][1]+Tran[safe[1][2]][2],HH):
												for i from 1 to nops(HH) do
													H0:=[op(H0), collect(H2[i][1][4], Vars, distributed)]:
												od:
												H0:=InvIrregAuto([op(SyzHH(H0)),op(GH2[2])], Syzbool):
												#H0:=[op(SyzHH(H0)),op(GH2[2])]:
											fi;
											return StrPBas(G1,H0,Syzbool,q, true);
										fi;
										maxsafe:=safe[s][1][1]:
										Superfl:=safe[s][1][2]:
										transk:=safe[s][2]:
									fi;
								od:
								print("Performing transformation:"),print(Tran[transk][1]),print(Tran[transk][1]+Tran[transk][2]);
								print("Index of safety:"),print(maxsafe);
								G2:=subs(Tran[transk][1]=Tran[transk][1]+Tran[transk][2],GG):
								G1:=[]:
								for i from 1 to nops(GG) do
									G1:=[op(G1),[collect(G2[i][1][4],Vars,distributed),collect(G2[i][2][3],Vars,distributed)]]:
								od:
								GH2:=InvRegAuto(SetGG(G1)):
								G1:=GH2[1]:
								H0:=[]:
								if not SigOnly then
									H2:=subs(Tran[transk][1]=Tran[transk][1]+Tran[transk][2],HH):
									for i from 1 to nops(HH) do
										H0:=[op(H0), collect(H2[i][1][4], Vars, distributed)]:
									od:
									H0:=InvIrregAuto([op(SyzHH(H0)),op(GH2[2])], Syzbool):
									#H0:=[op(SyzHH(H0)),op(GH2[2])]:
								fi;
								return StrPBas(G1,H0,Syzbool,q, true);	
							fi;
						fi;	
					
						#	print(1.5);
						#		G2:=subs(Tran[k][1]=Tran[k][1]+Tran[k][2],GG):
						#	#	print(1.7);
						#		H2:=[]:
						#		H0:=[]:
						#		if not SigOnly then
						#			if HH<>[] then
						#				H2:=subs(Tran[k][1]=Tran[k][1]+Tran[k][2],HH):
						#				for i from 1 to nops(H2) do
						#					H0:=[op(H0), collect(H2[i][1][4], Vars, distributed)]:
						#				od:
						#			fi;
						#		fi;
						#	#	print(2.1);
						#		G1:=[]:
						#		for i from 1 to nops(G2) do
						#	#	print(2.2);
						#			G1:=[op(G1),[collect(G2[i][1][4],Vars,distributed),collect(G2[i][2][3],Vars,distributed)]]:
						#	#		print(2.3);
						#		od:
						#	#	print(2.4);
						#		safe:=[op(safe), [IndexOfSafety(SetGG(G1),SyzHH(H0),q,Syzbool),k]]:
						#	#	print(2.6);
						#		print("safe"),##print(safe);
						#	od:
						#	print(3);
						#	maxsafe:=safe[1][1]:
						##	print(4);
							#transk:=safe[1][2]:
						#	print(5);
					#		if maxsafe=[] then
					#			print("Performing transformation:"),print(Tran[safe[1][2]][1]),print(Tran[safe[1][2]][1]+Tran[safe[1][2]][2]);
					#			print("Index of Safety:"),print(maxsafe);
					#			GG:=subs(Tran[safe[1][2]][1]=Tran[safe[1][2]][1]+Tran[safe[1][2]][2],GG):
					#			if not SigOnly then
					#				HH:=subs(Tran[safe[1][2]][1]=Tran[safe[1][2]][1]+Tran[safe[1][2]][2],HH):
					#			fi;
					#			G1:=[]:
					#			for i from 1 to nops(GG) do
					#				G1:=[op(G1),[collect(GG[i][1][4],Vars,distributed),collect(GG[i][2][3],Vars,distributed)]]:
					#			od:
					#			H0:=[]:
					#			if not SigOnly then
					#				for i from 1 to nops(HH) do
					#					H0:=[op(H0), collect(HH[i][1][4], Vars, distributed)]:
					#				od:
					#			fi;
					#			return StrPBas(SetGG(collect(G1,Vars,distributed)),SyzHH(collect(H0,Vars,distributed)),Syzbool,q, true);
					#		else
					#			for s from 2 to nops(safe) do
					#				if maxsafe<safe[s][1] or safe[s][1]=[] then
					#					if safe[s][1]=[] then
					#						maxsafe:=safe[s][1]:
					#						print("Performing transformation:"),print(Tran[safe[s][2]][1]),print(Tran[safe[s][2]][1]+Tran[safe[s][2]][2]);
					#						print("Index of Safety:"),print(maxsafe);
					#						GG:=subs(Tran[safe[s][2]][1]=Tran[safe[s][2]][1]+Tran[safe[s][2]][2],GG):
					#						if not SigOnly then
					#							HH:=subs(Tran[safe[s][2]][1]=Tran[safe[s][2]][1]+Tran[safe[s][2]][2],HH):
					#						fi;
					#						G1:=[]:
					#						for i from 1 to nops(GG) do
					#							G1:=[op(G1),[GG[i][1][4],GG[i][2][3]]]:
					#						od:
					#						H0:=[]:
					#						if not SigOnly then
					#							for i from 1 to nops(HH) do
					#								H0:=[op(H0), HH[i][1][4]]:
					#							od:
					#						fi;
					#						return	StrPBas(SetGG(collect(G1,Vars,distributed)),SyzHH(collect(H0,Vars,distributed)),Syzbool,q, true);
					#					fi;
					#					maxsafe:=safe[s][1]:
					#					transk:=safe[s][2]:
					#				fi;
					#			od:
					#			print("Performing transformation:"),print(Tran[transk][1]),print(Tran[transk][1]+Tran[transk][2]);
					#			print("Index of Safety:"),print(maxsafe);
					#			GG:=subs(Tran[transk][1]=Tran[transk][1]+Tran[transk][2],GG):
					#			HH:=subs(Tran[transk][1]=Tran[transk][1]+Tran[transk][2],HH):
					#			G1:=[]:
					#			for i from 1 to nops(GG) do
					#				G1:=[op(G1),[GG[i][1][4],GG[i][2][3]]]:
					#			od:
					#			H0:=[]:
					#			if not SigOnly then
					#				for i from 1 to nops(HH) do
					#					H0:=[op(H0), HH[i][1][4]]:
					#				od:
					#			fi;
					#			return StrPBas(SetGG(collect(G1,Vars,distributed)),SyzHH(collect(H0,Vars,distributed)),Syzbool,q, true);	
					#		fi;
					#	fi;
						
						
						GH:=InvRegAuto(sort([op(GG),rr],FHelp)):
						#if GH[2]<>[] then
						#	HH:=InvIrregAuto(sort([op(HH),op(GH[2])],FHelp),Syzbool):
						#fi;
   						HH:=[op(HH),op(GH[2])]:
                        GG:=sort(GH[1],FHelp):
						#BGG:=[op(BGG),op(InvBG(rr,GG,HH,Syzbool,q))]:
						BGG:=sort([op(BGG),op(InvBG(rr,GG,HH,Syzbool,q))],FHelp):
                        l2:= nops(BGG): 
					#	if ord=2 then
							for i from 1 to nops(GG) do
								if nops(GG[i][3])>0 then
									for m from 1 to nops(GG[i][3]) do 
										JPP:=[op(JPP),op(InvJP(GG[i],q,Syzbool))]:
										GG:=subsop(i=subsop(3=subsop(1=NULL,GG[i][3]),GG[i]),GG);##GG[i][3]:=subsop(1=NULL, GG[i][3]):
									od:
								fi;
							#	if ord=1 then
							#		JG:=FilljpPot(GG,q,Syzbool):
							#		JPP:=[op(JPP),op(JG[1])]:
							#		GG:=JG[2]:
							#	fi;
								JPP:=sort(JPP,FHelp):
							od:
					#	fi;
						l1:=nops(JPP):
					fi;
		       	else
					CountCrit:=CountCrit +1:
				fi;
	      	else
				CountCrit:=CountCrit +1:
			fi;
		else
			CountCrit:=CountCrit +1:
		fi;
	fi;
	ende:=0:
	while l1=0 and ende<nops(GG) do
		ende:=ende+1:
		
		#if ord=2 then
			for i from 1 to nops(GG) do
				if nops(GG[i][3])>0 then
					for m from 1 to nops(GG[i][3]) do 
						JPP:=[op(JPP),op(InvJP(GG[i],q,Syzbool))]:
						GG:=subsop(i=subsop(3=subsop(1=NULL,GG[i][3]),GG[i]),GG);##GG[i][3]:=subsop(1=NULL, GG[i][3]):
					od:
				fi;
			#	if ord=1 then
			#		JG:=FilljpPot(GG,q,Syzbool):
			#		JPP:=sort([op(JPP),op(JG[1])],FHelp):
			#		GG:=JG[2]:
			#	fi;
				#JPP:=[op({op(JPP)})]:
				JPP:=sort(JPP,FHelp):
			od:
		#fi;
		l1:=nops(JPP):
	od:
	#print(l1+l2);
od:

help:=nops(GG):
i:=1:
while not i > help do
	if InvDivG(GG[i][2][2],GG minus [GG[i]])<>[] then
		GG:=subsop(i=NULL,GG):
		i:=i-1:
		help:=nops(GG):
	fi;
	i:=i+1:
od:
if RestartBool then
	i:=1:
	while not i > nops(HH) do
		if InvIrregGG(HH[i],HH,Syzbool)[1] then
			HH:=subsop(i=NULL,HH):
			i:=i-1:
		fi;
		i:=i+1:
	od:
fi;
maxdg:=maxdeg(GG):
###########################
	secondtime, secondbytes := kernelopts(cputime, bytesused); 
printf("%-1s %1s %1s %1s:         %3a %3a\n", The, cpu, time, is, secondtime-firsttime, sec);
printf("%-1s %1s %1s %1s %1s: %3g\n", Num,of,regular, normal, forms,CountRed);
printf("%-1s %1s %1s %1s: %3a\n", Num,of,discarded, elements,CountCrit);
printf("%-1s %1s %1s %1s %1s: %3a\n",Total, Num,of, coord, trans,CountTran);
printf("%-1s %1s %1s: %3a\n", Num,of,syzygies,nops(HH));
printf("%-1s %1s %1s %1s: %3a\n", Num,of,ideal, elements,nops(GG));
printf("%-1s %1s %1s %1s %1s: %3a\n",Max, degree, of, ideal,element,maxdg);
return [HH,GG];
end:

#### Example ###########

CountTran:=0:
tord:= tdeg(z,y,x);
Vars:=[x,y,z];
SigOnly:=false; ### only keep signatures
Syzbool:=true;
ord:=2; 		# TOP-Lift			
q:=6;   	
H:=[];


#F:=[y^3-x*y*z,y*z^3+y^4,z^2*x+x^2*y,x^2+x*y]: ### PROBLEMS q=5

#F:=[y^2,y*z+y^2,x*z+x*y]: ### q=3
#F:=[y^3-x*y*z+y*z^2,y*z^3+y^2*x^2,z^2*x+x^2*y,x^2+x*y+z*x]: ##q=5

#F:=[y^3-x*y*z+y*z^2,y*z^3+y^2*x^2,z^2*x+x^2*y,x^2+x*y+z*x, x*y-t^2]: ##q=6 ### erste trans enthält weniger elemente die man anschauen muss...
#F:=[t*x*y^2+x^4-y*z^3, x^2*z*y-z*y^3+x*y*z*t, z^2*y^4+x^5*t, z^4*y+x^2*y^2*t]:
#F:=[z^3-y+x,y*z^2+y^2,y^10-x*y*z,x^2+x,y-z+x^2,x^3-y^3,x*y]:


#debug(StrPBas);
#debug(IndexOfSafety);
#HH3,GG3:= op(StrPBas(F,[],Syzbool,q,false)):

#Vars:=[u5,u4,u3,u2,u1,u0]:## katsura 5 - 23 elem q=5
#tord:=tdeg(u0,u1,u2,u3,u4,u5):
#f0:=u5*u5 + u4*u4 + u3*u3 + u2*u2 + u1*u1 + u0*u0 + u1*u1 + u2*u2 + u3*u3 + u4*u4 + u5*u5 - u0:
#f1:=u5*0 + u4*u5 + u3*u4 + u2*u3 + u1*u2 + u0*u1 + u1*u0 + u2*u1 + u3*u2 + u4*u3 + u5*u4 - u1:
#f2:=u5*0 + u4*0 + u3*u5 + u2*u4 + u1*u3 + u0*u2 + u1*u1 + u2*u0 + u3*u1 + u4*u2 + u5*u3 - u2:
#f3:=u5*0 + u4*0 + u3*0 + u2*u5 + u1*u4 + u0*u3 + u1*u2 + u2*u1 + u3*u0 + u4*u1 + u5*u2 - u3:
#f4:= u5*0 + u4*0 + u3*0 + u2*0 + u1*u5 + u0*u4 + u1*u3 + u2*u2 + u3*u1 + u4*u0 + u5*u1 - u4:
#f5:=u5 + u4 + u3 + u2 + u1 + u0 + u1 + u2 + u3 + u4 + u5 - 1:
#F:=[f0,f1,f2,f3,f4,f5];

#Vars:=[u6,u5,u4,u3,u2,u1,u0]:   # KATSURA 6 - 43 elements;; q=7 
#tord:=tdeg(u0,u1,u2,u3,u4,u5,u6):
#f0:=u0^2 - u0 + 2*u1^2 + 2*u2^2 + 2*u3^2 + 2*u4^2 + 2*u5^2 + 2*u6^2:
#f1:=2*u0*u1 + 2*u1*u2 - u1 + 2*u2*u3 + 2*u3*u4 + 2*u4*u5 + 2*u5*u6:
#f2:=2*u0*u2 + u1^2 + 2*u1*u3 + 2*u2*u4 - u2 + 2*u3*u5 + 2*u4*u6:
#f3:=2*u0*u3 + 2*u1*u2 + 2*u1*u4 + 2*u2*u5 + 2*u3*u6 - u3:
#f4:=2*u0*u4 + 2*u1*u3 + 2*u1*u5 + u2^2 + 2*u2*u6 - u4:
#f5:=2*u0*u5 + 2*u1*u4 + 2*u1*u6 + 2*u2*u3 - u5:
#f6:=u0 + 2*u1 + 2*u2 + 2*u3 + 2*u4 + 2*u5 + 2*u6 - 1:
#F:=[f0,f1,f2,f3,f4,f5,f6]:

#Vars:=[u7,u6,u5,u4,u3,u2,u1,u0]:   # KATSURA 7 - 79 elements;; q= 8
#tord:=tdeg(u0,u1,u2,u3,u4,u5,u6,u7):
#f0:=u0^2 - u0 + 2*u1^2 + 2*u2^2 + 2*u3^2 + 2*u4^2 + 2*u5^2 + 2*u6^2 + 2*u7^2:
#f1:=2*u0*u1 + 2*u1*u2 - u1 + 2*u2*u3 + 2*u3*u4 + 2*u4*u5 + 2*u5*u6 + 2*u6*u7:
#f2:=2*u0*u2 + u1^2 + 2*u1*u3 + 2*u2*u4 - u2 + 2*u3*u5 + 2*u4*u6 + 2*u5*u7:
#f3:=2*u0*u3 + 2*u1*u2 + 2*u1*u4 + 2*u2*u5 + 2*u3*u6 - u3 + 2*u4*u7:
#f4:=2*u0*u4 + 2*u1*u3 + 2*u1*u5 + u2^2 + 2*u2*u6 + 2*u3*u7 - u4:
#f5:=2*u0*u5 + 2*u1*u4 + 2*u1*u6 + 2*u2*u3 + 2*u2*u7 - u5:
#f6:=2*u0*u6 + 2*u1*u5 + 2*u1*u7 + 2*u2*u4 + u3^2 - u6:
#f7:=u0 + 2*u1 + 2*u2 + 2*u3 + 2*u4 + 2*u5 + 2*u6 + 2*u7 - 1:
#F:=[f0,f1,f2,f3,f4,f5,f6,f7]:

#Vars:=[w,z,y,x]: ###isaac97 - 12 elements ;;; q=5  
#tord:=tdeg(x,y,z,w): 
#f0:=8*w^2 + 5*x*w + 2*z*w + 3*w + 5*x^2 + 7*y^2 + 7*z^2 - 7*x + 2*x*y - 7*y - 4*w*y - 8*z - 7*x*z - 8*y*z + 8:
#f1:=3*w^2 + 9*w + 4*x^2 + 9*y^2 + 7*z^2 + 7*x - 5*w*x + 2*x*y + 5*y - 3*w*y + 6*y*z + 7*z - 6*w*z - 2*x*z + 5:
#f2:=- 2*w^2 + 9*x*w + 9*y*w + 8*x^2 + 6*y^2 - 4*w + 8*x + 9*x*y + 4*y + 8*z - 7*w*z - 3*x*z - 7*y*z - 6*z^2 + 2:
#f3:=7*w^2 + 5*x*w + 3*y*w + 2*x^2 - 5*w + 4*x + 9*x*y + 6*y - 4*y^2 - 9*z - 5*w*z - 7*x*z - 5*y*z - 4*z^2 + 2:
#F:=[f0,f1,f2,f3];

#Vars:=[h4,h3,h2,h1]: ### chandra 4 - 8 elem q=4
#tord:=tdeg(h1,h2,h3,h4):
#f0:=- 51234*h1^2 - 34156*h1*h2 - 25617*h1*h3 + 1497532*h1 - 1600000:
#f1:=- 170780*h1*h2 - 128085*h2^2 - 102468*h2*h3 + 7743830*h2 - 8000000:
#f2:=- 128085*h1*h3 - 102468*h2*h3 - 85390*h3^2 + 7829220*h3 - 8000000:
#f3:=- 717276*h1*h4 - 597730*h2*h4 - 512340*h3*h4 + 55103405*h4 - 56000000:
#F:=[f0,f1,f2,f3]:

#Vars:=[h5,h4,h3,h2,h1]: ### chandra 5 - 16 elem q=5
#tord:=tdeg(h1,h2,h3,h4,h5):
#f0:=- 256170*h1^2 - 170780*h1*h2 - 128085*h1*h3 - 102468*h1*h4 + 9487660*h1 - 10000000:
#f1:=- 170780*h1*h2 - 128085*h2^2 - 102468*h2*h3 - 85390*h2*h4 + 9743830*h2 - 10000000:
#f2:=- 896595*h1*h3 - 717276*h2*h3 - 597730*h3^2 - 512340*h3*h4 + 68804540*h3 - 70000000:
#f3:=- 1434552*h1*h4 - 1195460*h2*h4 - 1024680*h3*h4 - 896595*h4^2 + 138206810*h4 - 140000000:
#f4:=- 3586380*h1*h5 - 3074040*h2*h5 - 2689785*h3*h5 - 2390920*h4*h5 + 415696344*h5 - 420000000:
#F:=[f0,f1,f2,f3,f4]:

#Vars:=[h6,h5,h4,h3,h2,h1]: ### chandra 6 - 32 elem q=6
#tord:=tdeg(h1,h2,h3,h4,h5,h6):
#f0:=- 256170*h1^2 - 170780*h1*h2 - 128085*h1*h3 - 102468*h1*h4 - 85390*h1*h5 + 11487660*h1 - 12000000:
#f1:=- 1195460*h1*h2 - 896595*h2^2 - 717276*h2*h3 - 597730*h2*h4 - 512340*h2*h5 + 82206810*h2 - 84000000:
#f2:=- 1793190*h1*h3 - 1434552*h2*h3 - 1195460*h3^2 - 1024680*h3*h4 - 896595*h3*h5 + 165609080*h3 - 168000000:
#f3:=- 4303656*h1*h4 - 3586380*h2*h4 - 3074040*h3*h4 - 2689785*h4^2 - 2390920*h4*h5 + 498620430*h4 - 504000000:
#f4:=- 3586380*h1*h5 - 3074040*h2*h5 - 2689785*h3*h5 - 2390920*h4*h5 - 2151828*h5^2 + 499696344*h5 - 504000000:
#f5:=- 33814440*h1*h6 - 29587635*h2*h6 - 26300120*h3*h6 - 23670108*h4*h6 - 21518280*h5*h6 + 5504549820*h6 - 5544000000:
#F:=[f0,f1,f2,f3,f4,f5]:




Vars:=[y5,y4,y3,y2,y1]: #chemequs q=4, 10sec
tord:=tdeg(y1,y2,y3,y4,y5):
f0:=4006491052808310*y1*y2 + 526610813156747*y1 - 473964757724582*y5:
f1:=5313248634539330*y1*y2 + 349185128215561*y1 + 4313634906030440*y2^2 + 7431331204908770*y2*y3^2 + 806989849429520*y2*y3 + 149987275886495*y2*y4 + 132605295558012*y2 - 1047588595859020*y5:
f2:=18874386880707400*y2*y3^2 + 1024812258191070*y2*y3 + 478801780852764*y3^2 + 101454213888243*y3 - 1064282665103090*y5:
f3:=169336224248364*y2*y4 + 19504011585425050*y4^2 - 4730926576002525*y5:
f4:=13941268403110100*y1*y2 + 1832432069217330*y1 + 11318412924828000*y2^2 + 38997679215146700*y2*y3^2 + 4234871305042080*y2*y3 + 787093928408356*y2*y4 + 695877849581671*y2 + 989285552709683*y3^2 + 419243169406760*y3 + 45328426231998600*y4^2 - 203370256289491000000:
F:=[f0,f1,f2,f3,f4]:



#	Vars:=[r,z,y,x]: ## Bronstein-86 ##q=5 1.56sec
#	tord:=tdeg(x,y,z,r):
#	f0:=x*y+z^2-1:
#	f1:=x^2+y^2+z^2-r^2:
#	f2:= x*y*z-x^2-y^2-z+1:
#	F:=[f0,f1,f2]:

#	Vars:=[z, y,x, w]: ## Cyclic_4 ## q=6 278sec
#	tord:=tdeg(w, x, y, z):
#	f0:=w+x+y+z:
#	f1:=w*x+x*y+w*z+y*z:
#	f2:=w*x*y+w*x*z+w*y*z+x*y*z:
#	f3:=w*x*y*z-1:
#	F:=[f0,f1,f2,f3]:

					#	Vars:=[z,y,x]: ###Czapor-86a  q=4 0.078sec
					#	tord:=tdeg(x,y,z):
					#	f0:=5*x^2+8*x*y+9*y^2+4*x*z-6*y*z-z^2+8*x+2*y-7*z+5:
					#	f1:=8*x^2-2*x*y+3*y^2-6*x*z-7*y*z+10*z^2+3*x+10*y-8*z-4:
					#	f2:=10*x^2-2*x*y+9*y^2+6*x*z-y*z-2*z^2-6*x-4*y+5*z-9:
					#	F:=[f0,f1,f2]:

					#	Vars:=[t,z,y,x]: ###Caprasse  q=7 3400sec
					#	tord:=tdeg(x,y,z,t):
					#	f0:=2*y*z*t+x*t^2-x-2*z:
					#	f1:=y^2*z+2*x*y*t-2*x-z:
					#	f2:=-x*z^3+4*y*z^2*t+4*x*z*t^2+2*y*t^3+4*x*z+4*z^2-10*y*t-10*t^2+2:
					#	f3:=-x^3*z+4*x*y^2*z+4*x^2*y*t+2*y^3*t+4*x^2-10*y^2+4*x*z-10*y*t+2:
					#	F:=[f0,f1,f2,f3]:




Vars:=[h6,h5,h4,h3,h2,h1]: ### chandra 6 - 32 elem q=6
tord:=tdeg(h1,h2,h3,h4,h5,h6):
f0:=- 256170*h1^2 - 170780*h1*h2 - 128085*h1*h3 - 102468*h1*h4 - 85390*h1*h5 + 11487660*h1 - 12000000:
f1:=- 1195460*h1*h2 - 896595*h2^2 - 717276*h2*h3 - 597730*h2*h4 - 512340*h2*h5 + 82206810*h2 - 84000000:
f2:=- 1793190*h1*h3 - 1434552*h2*h3 - 1195460*h3^2 - 1024680*h3*h4 - 896595*h3*h5 + 165609080*h3 - 168000000:
f3:=- 4303656*h1*h4 - 3586380*h2*h4 - 3074040*h3*h4 - 2689785*h4^2 - 2390920*h4*h5 + 498620430*h4 - 504000000:
f4:=- 3586380*h1*h5 - 3074040*h2*h5 - 2689785*h3*h5 - 2390920*h4*h5 - 2151828*h5^2 + 499696344*h5 - 504000000:
f5:=- 33814440*h1*h6 - 29587635*h2*h6 - 26300120*h3*h6 - 23670108*h4*h6 - 21518280*h5*h6 + 5504549820*h6 - 5544000000:
F:=[f0,f1,f2,f3,f4,f5]:

Vars:=[c,z,y,x]; ##Noonburg-89
tord:=tdeg(x, y, z, c);
f0:=x^2*z+y^2*z-z*c+1:
f1:= x*y^2+x*z^2-x*c+1:
f2:= x^2*y+y*z^2-y*c+1:
F:=[f0,f1,f2]:
							
ord:=2:						
q:=6:
SigOnly:=true:
Syzbool:=false:
#debug(StrPBas);
HH3,GG3:=op(StrPBas(F,[],Syzbool,q,false)):


##### If you want to check your output copy this to your Maple Worksheet.
for i from 1 to nops(GG3) do
GG3[i][3]:=JVars(GG3[i][2][2]):
GG3[i][1][2],GG3[i][1][3],GG3[i][2][2]:
od:
q2:= q+3:
TestBasis([GG3,HH3],Syzbool,q2); #### q2<q is forbidden.