-- Hdepth 1.0
-- Copyright (C) Bogdan Ichim and Andrei Zarojanu
-- This program is free software, it is distributed in the hope that
-- it will be useful, but WITHOUT ANY WARRANTY.
-- Please refer to Hdepth and the related article in any
-- publication for which it has been used.

Define ByN(S,T)    -- define the sorting function
  If sum(S) < sum(T) Then Return True;
  Elif sum(S)=sum(T) Then
	For i:=1 To Len(S) Do
		If S[i]>T[i] Then
			Return True;
		Elif S[i]<T[i] Then
			Return False;
		EndIf;
	EndFor;
  
Elif sum(S) > sum(T) Then Return False;
  EndIf;
    Return True;
EndDefine; -- ByN




 Define rho(L,G) -- rho function
TopLevel CurrentRing;
  N:=NumIndets(CurrentRing);
b:=0;
For k:=1 To N  Do
  If L[k]=G[k] Then b:=b+1; 
 EndIf;
 	EndFor;
 Return b;
	 EndDefine; -- rho

Define divide(A,B) -- check if monom(A) divides monom(B)
TopLevel CurrentRing;
N:=NumIndets(CurrentRing);
For k:=1 To N  Do
	If A[k] > B[k] Then Return false;
	EndIf;
EndFor;
Return true;
EndDefine; --divide

Define monom(L) -- creates the monomial from a list 
TopLevel CurrentRing;
  N:=NumIndets(CurrentRing);
m:=1;
L1:=indets(CurrentRing);
For k:=1 To N Do
m:=m*L1[k]^L[k];
EndFor;
Return m;
EndDefine; --monom
	





    Define FETC(s,T,G)  -- build B_{<s} ( the list of monomials with the rho function < s )
	TopLevel CurrentRing;
  N:=NumIndets(CurrentRing);
L1:=[0];
	For i:=1 To Len(T) Do
		If rho(T[i],G) < s Then
		L1:=ConcatLists([L1,[T[i]]]);
EndIf;
EndFor;
      Return tail(L1); 
    EndDefine; -- FETC



 Define FPC(a,F)  -- build B_{=s} ( the list of monomials which stand over a and have rho function = s)
	TopLevel CurrentRing;

L2:=[0];

For i:=1 To Len(F) Do
	If divide(a,F[i]) Then L2:=ConcatLists([L2,[F[i]]]);
	EndIf;
EndFor;

Return tail(L2); 
      EndDefine; -- FPC

Define TestNatural(P) -- check if P \in N
TopLevel CurrentRing;
  Coeffs := coefficients(P);
  For i:= 1 To Len(Coeffs) Do
    If Coeffs[i] < 0 Then Return False;
    EndIf;
  EndFor;
  
 Return True;
EndDefine;  -- TestNatural

Define BuildInterval(B,A,st, k,Ref s) -- build the polynomial corresponding to the interval [A,B] 
TopLevel CurrentRing;
N:=NumIndets(CurrentRing);
  
  While k >= 0 Do
    If k=N Then s:=s+monom(st); k:=k-1;
    EndIf;
    If st[k+1] < B[k+1] Then
      st[k+1]:=st[k+1]+1;
      k:=k+1;
 Else
 	st[k+1]:=A[k+1];
      k:=k-1;
  EndIf;
  EndWhile;

  EndDefine; -- BuildInterval


  

Define Deduction(P1,P2,E) -- We test if the list E changes after the deduction of P2 from P
TopLevel CurrentRing;
  L1:=support(P1);
L2:=support(P2);
For i:=1 To len(L2) Do
	If L2[i] IsIn L1 Then nimic:=0;
	Else b:=log(L2[i]);  E:=diff(E,[b]);
	EndIf;
EndFor;
Return E;
EndDefine; --Deduction

Define Deduction2(P1,b,F) -- We test if the list F changes after the deduction of P2 from P
TopLevel CurrentRing;
  L1:=support(P1);
	If monom(b) IsIn L1 Then nimic:=0;
	Else F:=diff(F,[b]);
	EndIf;
Return F;
EndDefine; --Deduction2



Define CHD(P,E,F) -- Check Hilbert Depth (backtracking)
TopLevel CurrentRing;
  If TestNatural(P) = false Then Return false;
EndIf;
 N:=NumIndets(CurrentRing);
  st:=NewList(N,1);

  If len(E) = 0 Then   Return true;
Else
	For i:=1 To len(E) Do
		C:=FPC(E[i],F);
		If len(C) = 0 Then Return false;
		EndIf;
		For j:=1 to Len(C) Do
			P2:=0;
		  D:=E[i]-st;
			BuildInterval(C[j],D,D,0,Ref P2);
		
			P1:= P - P2;
		 	t:= CHD(P1,Deduction(P1,P2,E),Deduction2(P1,C[j],F));
		  If t= true Then   PrintLn P2;
		    Return true;
			EndIf;
		EndFor;
	EndFor;
Return false;
EndIf;
EndDefine; -- CHD




Define Start(g,s,P) -- start function
TopLevel CurrentRing;
  TopLevel ByN;
F:=[0];
G:=log(g);
L:=support(P);
T:=NewList(Len(L));
For k:=1 To len(L)  Do
  T[k]:=log(L[k]);
	If rho(T[k],G)=s Then F:=ConcatLists([F,[T[k]]]); -- create F = (m | rho(m) = s)
	EndIf;
EndFor;
 E:=FETC(s,T,G); -- create  E = (m | rho(m) < s)
 reverse(ref E);
  F:=tail(F);
Return CHD(P,E,F);
EndDefine; --Start