bisystemsolving[LGP] := proc() `bisystemsolving/LGP`(args) end:
bisystemsolving[LGP_from_symbolic_solution] := proc() `bisystemsolving/LGP_from_symbolic_solution`(args) end:

with(RootFinding):
trans_root_form:=proc(Osolution,Isolution,k)
local i,j,n,tempR,output;
        if k<1 then return Osolution;fi;  
	n:=nops(Isolution);
        output:=[];
        if k>1 then
        	for i from 1 to n do
			tempR:=[op(Osolution),Isolution[i][1]];
                	output:=[op(output),op(trans_root_form(tempR,Isolution[i][2],k-1))]; 
        	od;
        elif k=1 then 
        	for i from 1 to n do
			tempR:=[op(Osolution),Isolution[i]];
                	output:=[op(output),tempR]; 
        	od;
        else
             return Osolution; 
        fi;  
return output;  
end:

###################################################################################
compute_upper_bound:=proc()
local i,j,n,temp,T,tempR,typ,Rbound,ninterval,ii,Isolution,poly,variablelist,symbol;
Isolution:=args[1];
poly:=args[2];
variablelist:=args[3];
if nargs>3 then
   symbol:=args[4];
else
   symbol:=1;
fi;
       ii:=nops(variablelist)-1; 
       tempR:=trans_root_form([],Isolution,ii);
       n:=nops(tempR);
       if n>0 then
          if type(tempR[1][1],'numeric')=true then
		typ:='numeric';
	  else
		typ:='interval'; 
          fi;
       fi;
       Rbound:=1;  
       for i from 1 to n do
	   temp:=tempR[i];
           ninterval:=[];
           for j from 1 to ii do
                if typ='numeric' then
			T:=temp[j]*1000;
                        ninterval:=[op(ninterval),[floor(T)/1000.0,ceil(T)/1000.0]];
                else
                        ninterval:=[op(ninterval),evalf(temp[j])];
                fi;
           od;
           Rbound:=max(Rbound,single_upper_bound(ninterval,poly,variablelist,symbol));
       od;
       return Rbound;
end:

#######################################################################################
single_upper_bound:=proc()
local nintervals,f,variablelist,symbol,intersign,chnintervals,temp,currentterm,ltsign,tsign,degvar,L,TS,tb,tpre,varnum,i,j,termnum,pfu,pfd,nfu,nfd;
nintervals:=args[1];
f:=args[2];
variablelist:=args[3];
if nargs>3 then
   symbol:=args[4];
else
   symbol:=1;
fi;
     varnum:=nops(variablelist);
     termnum:=nops(f);    
     nfu:=0;
     nfd:=0;
     pfu:=0;
     pfd:=0;
     intersign:=[];
     chnintervals:=[];
     tpre:=1;
     for j from 1 to varnum-1 do
            temp:=sign(op(1,op(j,nintervals)));
            intersign:=[op(intersign),temp];    
            if temp=-1 then
              chnintervals:=[op(chnintervals),[-op(2,op(j,nintervals)),-op(1,op(j,nintervals))]];
            else
              chnintervals:=[op(chnintervals), op(j,nintervals)]; 
            fi;
     od;
     for i from 1 to termnum do 
        if termnum>1 then
           currentterm:=op(i,f);
        else
           currentterm:=f;
        fi;
        ltsign:=sign(lcoeff(currentterm));
        tsign:=ltsign;
        for j from 1 to varnum-1 do
            degvar:= degree(currentterm,variablelist[j]) mod 2;
            L:=intersign[j]^degvar;
            tsign:=tsign*L;
        od;
        if tsign=1 then
                pfu:=pfu+ltsign*subs(seq(op(j,variablelist)=chnintervals[j][2],j=1..varnum-1),currentterm); 
                pfd:=pfd+ltsign*subs(seq(op(j,variablelist)=chnintervals[j][1],j=1..varnum-1),currentterm); 
                if degree(currentterm,variablelist[varnum]) mod 2=1 then
                    nfu:=nfu+ltsign*subs(seq(op(j,variablelist)=chnintervals[j][1],j=1..varnum-1),currentterm); 
                    nfd:=nfd+ltsign*subs(seq(op(j,variablelist)=chnintervals[j][2],j=1..varnum-1),currentterm);
                else
                    nfu:=nfu+ltsign*subs(seq(op(j,variablelist)=chnintervals[j][2],j=1..varnum-1),currentterm); 
                    nfd:=nfd+ltsign*subs(seq(op(j,variablelist)=chnintervals[j][1],j=1..varnum-1),currentterm);
                fi;
        else
                pfu:=pfu-ltsign*subs(seq(op(j,variablelist)=chnintervals[j][1],j=1..varnum-1),currentterm); 
                pfd:=pfd-ltsign*subs(seq(op(j,variablelist)=chnintervals[j][2],j=1..varnum-1),currentterm); 
                if degree(currentterm,variablelist[varnum]) mod 2=1 then
                    nfu:=nfu-ltsign*subs(seq(op(j,variablelist)=chnintervals[j][2],j=1..varnum-1),currentterm); 
                    nfd:=nfd-ltsign*subs(seq(op(j,variablelist)=chnintervals[j][1],j=1..varnum-1),currentterm);
                else
                    nfu:=nfu-ltsign*subs(seq(op(j,variablelist)=chnintervals[j][1],j=1..varnum-1),currentterm); 
                    nfd:=nfd-ltsign*subs(seq(op(j,variablelist)=chnintervals[j][2],j=1..varnum-1),currentterm);
                fi;
        fi;    
     od;
if symbol=1 then
     TS:=[pfu,pfd,nfu,nfd]; 
     tb:=1;
     for i from 1 to 4 do
         temp:=map(t->op(2,t),Isolate(TS[i]));
         if temp<>[] then
            tb:=max(abs(temp[1]),abs(temp[-1]),tb);
         fi;
     od;
     temp:=2*ceil(tb);
     if temp=tb then
        temp:=temp+1;
     fi;
else
    temp:=2*max(rootbound(pfu,variablelist[-1]),rootbound(pfd,variablelist[-1]),rootbound(nfu,variablelist[-1]),rootbound(nfd,variablelist[-1]));
fi;
return temp;
end:

####################################################################################
compute_sep_bound:=proc()
local i,n,sep,temp,rr,typ,rootlist;
if nargs>=1 then
   rootlist:=args[1];
else
   ERROR("The function has at least one parameter!");
fi;
if nargs>1 then
   typ:=args[2];
else
   typ:='numeric';
fi;
n:=nops(rootlist);
sep:=10;
if n=1 or n=0 then 
   return sep;
fi;
for i from 1 to n-1 do
   if typ='interval' then
      temp:=rootlist[i+1][1]-rootlist[i][2];
   elif typ='numeric' then
      temp:=rootlist[i+1]-rootlist[i];
   else
      ERROR("The type should be 'interval' or 'numeric'.");
   fi; 
   sep:=min(sep,temp);
od;
rr:=sep;
temp:=1;
while floor(sep)=0 do
   sep:=10*sep;
   temp:=temp*10;
od;
	sep:=floor(sep)/temp;
rr:=sep;
return rr;
end: #compute_sep_bound


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

myabs:=proc(A)
if type(A,'numeric')=true then
   return abs(A);
elif type(A,list)=true then
   return max(abs(A[1]),abs(A[2]));
fi;
end:

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

myminus:=proc(A,B)
if type(A,'numeric')=true and type(B,'numeric')=true then
   return A-B;
elif type(A,list)=true and type(B,list)=true then
   return [A[1]-B[2],A[2]-B[1]];
else
   ERROR("The type of the input parameters isn't fit the function: myminus, please check!");
fi;
end:

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

mynumber_multiply:=proc(num,A)
if type(A,'numeric')=true then
   return num*A;
elif type(A,list)=true then
   return [num*A[1],num*A[2]];
fi;
end:

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

rootgroup_division:=proc()
local i,j,k,nX,nY,temp,tempR,Rg,X,Y,r,s,typ;
if nargs<4 then
   ERROR("There are at least 4 parameters for this function");
else
   X:=args[1];
   Y:=args[2];
   r:=args[3];
   s:=args[4];
fi;
if nargs>=5 then
   typ:=args[5];
else
   typ:='numeric';
fi;
nX:=nops(X);
if nX=0 then 
	ERROR("The first parameter can be empty list!");
fi;
nY:=nops(Y);   
Rg:=[];
tempR:=[];
j:=1;
k:=0;
i:=1;
while i<=nY do
  temp:=myminus(Y[i],X[j]);  
  if myabs(temp)<r/2 then
     tempR:=[op(tempR),mynumber_multiply(1/s,temp)];
     i:=i+1;
  else
     k:=i;
     Rg:=[op(Rg),tempR];
     j:=j+1;
     tempR:=[];
  fi;
  if j>nX then break;fi;
od;
Rg:=[op(Rg),tempR];
j:=j+1;
if j<=nX then
   for i from j to nX do
       Rg:=[op(Rg),[]];
   od;
fi; 
if Rg=[[]] then Rg:=[];fi;
return Rg;
end: #rootgroup_division

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

mapxy:=proc()
global r_sep_bound;
# r_sep_bound is the root bound, we need to compute r many times and choose the minimal one.
local X,Y,bl,bc,n,i,T,r1,typ,temp;
if nargs<2 then
   ERROR("There are at least two parameters!");
else
   X:=args[1];
   Y:=args[2];
fi;
if nargs>2 then
   bl:=args[3];
else
   bl:=false;
fi;
if nargs>3 then
   bc:=args[4];
else
   bc:=true;
fi;
n:=nops(X);
if nops(Y)<>n then
	ERROR("the number of two lists don't match!");
fi;

i:=1;
while i<=n and Y[i]=[] do
i:=i+1;
od;
if i>n then #all Y[i] are empty
	return [];
else
	temp:=Y[i];
fi;
if n>0 and type(temp[1],'numeric')=true then
	typ:='numeric';
else
	typ:='interval';
fi;
T:=[];
for i from 1 to n do
   if Y[i]<>[] then
        if bc=true then   
        	r1:=compute_sep_bound(Y[i],typ);
                r_sep_bound:=min(r_sep_bound,r1);
        fi; 
   	T:=[op(T),[X[i],Y[i]]];
   else
        if bl=false then
	        T:=[op(T),[]];  
        fi;   
   fi;
od;
return T;
end:

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


SingleDimensionalRootSolve:=proc(Isolution,U,rlist,slist,typ)
local i,k,n,T,tempR,tempS,crlist,cslist,output;
n:=nops(rlist);
if nops(slist)<>n then
   ERROR("rlist and slist are not matching!");
fi;
if n=1 then 
	tempR:=rootgroup_division(Isolution,U,rlist[1],slist[1],typ);
        if nops(tempR)<>nops(Isolution) then
		ERROR("The root number is not match, the problem may be the root bound not correct.");
        fi;
        T:=mapxy(Isolution,tempR,true,false);	
else
        T:=[];
        tempR:=rootgroup_division(map(t->op(1,t),Isolution),U,rlist[1],slist[1],typ);
        crlist:=rlist[2..-1];
        cslist:=Rlist[2..-1];
        k:=nops(tempR);
        for i from 1 to k do
               if tempR[i]<>[] then 
		       tempS:=SingleDimensionalRootSolve(Isolution[i][2],tempR[i],crlist,cslist,typ);			       	       
                       T:=[op(T),[Isolution[i][1],tempS]];
               else
	       	       T:=[op(T),[]];
               fi;     
        od;
fi;
output:=[];
for i from 1 to nops(T) do
    if T[i]<>[] then
	output:=[op(output),T[i]];
    fi;
od;
return output;
end:

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

compute_y_bound:=proc(f1,f2,x,y)
local temp,L,R;
L:=map(t->op(2,t),Isolate(resultant(f1,f2,x)));
if L<>[] then
	temp:=max(abs(L[1]),abs(L[-1]));
fi;
if temp<1 then 
	R:=1;
else
        R:=ceil(temp);
        if R=temp then R:=R+1;fi;
fi;
return R;
end: 

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

  
`bisystemsolving/LGP`:=proc()
local polysystem,varlist,pre,cpre,ipoly,typ,r,r2,R,Tsolve,i,j,k,n,st,temp,X,Y,Isolution,resh,tp,s,bsymbolic,tt;
if nargs<2 then
  error("There are at least two parameters: 1. polynomial system; 2. variable list.");
fi;
polysystem:=args[1];
varlist:=args[2];
if nargs>2 then 
   if args[3]='symbolic' then
	bsymbolic:=true;
        pre:=0; 
   elif type(args[3],'integer') then
        bsymbolic:=false;
	pre:=args[3];
   else
	ERROR("The third parameter should be 'symbolic' or a positive integer!");
   fi;
else
   bsymbolic:=false;
   pre:=10;
   typ:='numeric';
fi;
if nargs>3 then 
   typ:=args[4];
else
   typ:='numeric'; 
fi;
n:=nops(polysystem);
if n<>2 then ERROR("The input system should contain two bivariate polynomials only!");fi;
if nops(varlist)<>n then
	ERROR("The system and variable list are not matching!");
fi;
st:=time():
resh:=resultant(polysystem[1],polysystem[2],varlist[2]);
if [op(indets(resh))]=[] then return [];fi;
cpre:=[pre+3];
Isolution:=Isolate(resh,output='interval',digits=cpre[1],format=list);
if Isolution=[] then 
 return [];
fi;
    tt:=time();
    r:=compute_sep_bound(Isolution,interval);
    if r<=0 then
        printf("WARNING! The output intervals of the function 'Isolate' are not disjoint!");
	temp:=cpre[1];
	while r<=0 do
		temp:=temp+3;
		Isolution:=Isolate(resh,output='interval',digits=temp,format=list);
		r:=compute_sep_bound(Isolution,interval);		
	od;
    fi; 
     
    if degree(polysystem[1],varlist[2])>=degree(polysystem[2],varlist[2]) then
    	temp:=compute_upper_bound(Isolution,polysystem[2],varlist,2); 
    else
        temp:=compute_upper_bound(Isolution,polysystem[1],varlist,2);
    fi;
    if r<1 then 
        R:=temp+2;
    else
        R:=temp+2*ceil(r);
    fi;
    if r>10^(4-pre) then 
    	tp:=R/r;
    else
	if r>2/10^pre then
		tp:=R/(r-1/10^pre);
	else
		tp:=2*R/r;
		cpre:=[ceil(abs(log[10](r/2)))+3];
	fi;
    fi;		
    if tp>1 then 
       s:=1/ceil(tp);
    else
       s:=1;
    fi;       
    ipoly:=resultant(subs(varlist[1]=X-s*Y,varlist[2]=Y,polysystem[1]),subs(varlist[1]=X-s*Y,varlist[2]=Y,polysystem[2]),Y);
    if [op(indets(ipoly))]=[] then return [];fi;
if bsymbolic then
    return [resh,ipoly,r,s];
else
    tp:=ceil(evalf(log[10](2/s)));
    cpre:=[op(cpre),max(cpre[1]-3,tp+pre)];
    Tsolve:=Isolate(ipoly,output='interval',digits=cpre[2],format=list);
    r2:=compute_sep_bound(Tsolve,interval);
    tp:=max(ceil(evalf(log[10](2/r2))),tp+pre,cpre[1]); 
    if tp>pre+3 then
    	Isolution:=Isolate(resh,output='interval',digits=tp,format=list);
    fi;     
    Isolution:=SingleDimensionalRootSolve(Isolution,Tsolve,[r],[s],'interval');
    Tsolve:=trans_root_form([],Isolution, 2);
fi;
if typ='numeric' then
    Tsolve:=map(t->[evalf(t[1][1]+t[1][2])/2,evalf(t[2][1]+t[2][2])/2],Tsolve);  	  
fi;
printf("The computing time is %a seconds!",time()-st);
return Tsolve;
end: # LGP


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

`bisystemsolving/LGP_from_symbolic_solution`:=proc()
local h,T,r,s,sep,pre,typ,rpre,temp,solution1d,solution2d,Isolution; 
if nargs<1 then
  error("There are at least one parameter: [h,T,r,s].");
fi;
h:=args[1][1];
T:=args[1][2];
r:=args[1][3];
s:=args[1][4];
if nargs>1 then 
	pre:=args[2];
else
	pre:=10;
        typ:='numeric';
fi;
if nargs>2 then 
   typ:=args[3];
else
   typ:='numeric'; 
fi;
temp:=ceil(evalf(log[10](2/s)));
rpre:=pre+temp;
solution1d:=Isolate(h,output=typ,digits=rpre,format=list);
if solution1d=[] then 
 return [];
fi;
solution2d:=Isolate(T,output=typ,digits=rpre,format=list);
if typ='interval' then 
	sep:=compute_sep_bound(solution2d,interval);
	if sep<=1/10^rpre then
		solution1d:=Isolate(h,output=interval,digits=ceil(abs(log[10](sep/2))),format=list);
	fi;
fi;
Isolution:=SingleDimensionalRootSolve(solution1d,solution2d,[r],[s],typ);
Isolution:=trans_root_form([],Isolution, 2);
return Isolution;
end: # LGP_from_symbolic_solution 
####################################################################


print("*******************************************************************");
printf("Coded by CHENG Jin-San: jcheng@amss.ac.cn\n"); 
printf("Copyright@ Key Laboratory of Mathematics Machanization, Chinese Academy of Sciences, 2008.\n");
with(bisystemsolving);
printf("If you need help about how to use the functions, please type 'HELP();'."); 
HELP:=proc()
local bshowhelp;
bshowhelp:=true:
if bshowhelp then
	print("****************************HELP***********************************");
	printf("1.  a. LGP([f(x,y),g(x,y)],[x,y],precision in digits,solution type).\n");
	printf("       precision in digits: an integer p, the output precision for roots 10^(-p), default is 10.\n");
	printf("       solution type: 'interval' or 'numeric' to point out the requirement for output solution, default is 'numeric'.\n");
	printf("    b. LGP([f(x,y),g(x,y)],[x,y],'symbolic').\n");
	printf("       This case just outputs a symbolic solution for the system [h(x),T(X),r,s].\n");
	print("*******************************************************************");
	printf("2.  LGP_from_symbolic_solution(T,precision in digits,solution type).\n");
	printf("    Where T:=LGP([f(x,y),g(x,y)],[x,y],'symbolic'); precision in digits,solution type are similar as the first function!");
	print("*******************************************************************");
fi;
end: