SDPTools:=module()
	description "Tools for solving SDP in high presicion";
	export 
		SqrtMatrix,
		SolveSDP,
		BigM;
	local
        RootBisection,
		PotentialFunction,
		Solution_of_Directions,
		Solve_L,
		Solve_U,
		LDL_Decomposition,
		Expresstion_of_PlaneSearch,
		Feasible_Rectangle,
        Feasible_Rectangle2,
		StepsLength_PlaneSearch;
	
	option
		package;



	#============================================================ 
	# Discription: 
	#     For a given semidefinite positive matrix M, the function 
	#     return the square root of  M. 
	# Input:  
	#     M  - a given semidefinite positive matrix 
	# Output: 
	#        - the square root of M 
	#============================================================
	SqrtMatrix:=proc(M::Matrix)
		local i,j,T_,Z_;
   		(T_, Z_) := SchurForm(M, output=['T', 'Z']);
   		for i from 1 to RowDimension(T_) do
			if T_[i,i]<0 then 
				error("Invalid input, the input matrix is not semipositive");
			else
				T_[i,i]:=sqrt(T_[i,i]);
				for j from i+1 to ColumnDimension(T_) do
					T_[i,j]:=0;
				od;
			fi;
		od; 
   		return Z_.T_.Transpose(Z_);
	end proc:



	#=========================================================================================================== 
	# Discreption: 
	#    Return the value of potential function at an iteration point and the gap in primal-dual problem 
	# Inputs: 
	#    x, Z  - an iteration point  
	#    F0    - the matice F0 in the primal problem 
	#    L	   - the matrces F1..Fm in the primal problem
  	#    v     - the parameter that sets the ralative weight of the term involving duality gap and the term which 
 	#            is the deviation from centrality
	# Outputs: 
	#    f     - the value of potential function at an iteration point 
	#    gap   - the gap in primal-dual problem 
	#==========================================================================================================
	PotentialFunction:=proc(x::list,Z::Matrix,F0::Matrix,L::list,v)  
  		local i,j,F,f,n,gap; 
  		F:=F0;   n:=ColumnDimension(F0); 
  		F:=add(x[i]*L[i],i=1..nops(x))+F;
  		gap:=0;
  		for i from 1 to n do
  	 		for j from 1 to n do
      				gap:=gap+F[i,j]*Z[j,i];
    			end do;  
  		end do;
 		f:=(n+v*sqrt(n))*log(gap)-log(Determinant(F))-log(Determinant(Z))-n*log(n);    
  		return gap; 
	end proc: 


	#======================================================= 
	# Discription: 
	#    To get the solution of the suitable direction  
	#     dx and dZ of the objective function(57), P31 
	# Inputs: 
	#    S  - A semidefinite n*n matrix 
	#    D  - A symmetric n*n matrix 
	#    L  - A list of n*n matrice 
	# Outputs: 
	#    dX  - the suitable direction dX  #    dZ  - the suitable dual feasible direction dZ 
	#====================================================== ===============================
  	######################################################################################  
	Solution_of_Directions:=proc(S,D,L::list) 
    		local i,j,p,q,A,AA,b,bb,m,n,dX,dZ,Inv_S,M,N,LL,d;  
    		#print(S);
    		#print(D);
    		m:=nops(L); 
    		#print(m);
    		n := ColumnDimension(D); 
   		# A := Matrix(m);  
    		A:=Matrix(m);
    		b := Vector(m); 
   		# b:=[]; 
   		# M:=[seq(Matrix(n),i=1..m)];
		printf("begin to compute inverse \n");
    		Inv_S := MatrixInverse(S);
    		printf("begin to compute M \n");
    		for i from 1 to m do
      			M[i]:=Inv_S.L[i];
    		end do;
    		printf("End to compute M \n"); 
    		N:=Inv_S.D; 
    		printf("begin to compute A \n");
    		for i from 1 to m do  
      		 	printf("begin %d row",i);
       			for j from 1 to i do 
          			for p from 1 to n do 
            				for q from 1 to n do
              					A[i,j]:=A[i,j]+M[i][q,p]*M[j][p,q];
            				end do;
          			end do;
       			end do; 
       			for p from 1 to n do
         			for q from 1 to n do
           				b[i]:=b[i]-M[i][q,p]*N[p,q];
         			end do;
       			end do; 
    		end do;  
    		for i from 1 to m do
      			for j from i+1 to m do
        			A[i,j]:=A[j,i];
      			end do;
    		end do;  #return A,b;
    		printf("end to compute A");
    		LL,d:=LDL_Decomposition(A);
    		dX:=Solve_L(LL,b);
    		for i from 1 to Dimension(dX) do
      			dX[i]:=dX[i]/d[i];
    		end do;
    		dX:=Solve_U(Transpose(LL),dX);
    		dZ :=-1*D; 
    		for i from 1 to m do  
       			dZ:=dZ-dX[i]*L[i];   
    		end do; 
    		dZ:=evalf(Inv_S.dZ.Inv_S);    
    		return dX,dZ;   
	end proc:





	###############################################################################
	# Discription:
	#   solve eqs A.x=b, where A=LU
	###############################################################################
	Solve_L:=proc(L,b::Vector)
  		local i,j,k,m,n,y;
  		m,n:=Dimension(L);
  		y:=b;
  		for k from 1 to n do
    			if L[k,k]=0 then
      				if y[k]<>0 then
        				error("this system has no solution");
      				end if;
    			else
      				y[k]:=y[k]/L[k,k];
      				y[k+1..n]:=y[k+1..n]-y[k]*L[k+1..n,k];
    			end if;
  		end do;
  		return y;
	end proc:



	#############################################################################
	#
	#############################################################################
	Solve_U:=proc(U,b::Vector)
  		local k,m,n,x;
  		m,n:=Dimension(U);
  		if m<>n then 
    			error("input matrix must be a squre matrix");
  		end if;
  		x:=b;
  		for k from n by -1 to 1 do
    			if U[k,k]=0 then
      				if x[k]<>0 then
        				error("this system has no sultion");
      				end if;
    			else
      				x[k]:=x[k]/U[k,k];
      				x[1..k-1]:=x[1..k-1]-x[k]*U[1..k-1,k];
    			end if;
  		end do;
  		return x;
	end proc:
  
  
	##########################################################################
	# Discription:
	#   Given a semidefinite positive matrix A, return it's LDL' decomposition
	# Author: Feng Guo   Date: 1/9/2009
	##########################################################################
	LDL_Decomposition:=proc(A)
  		local i,j,k,m,n,L,d;
  		m,n:=Dimension(A);
  		if m<>n then
    			error("the input matrix must be a square matrix");
  		end if;
  		L:=Matrix(n); d:=Vector(n);
  		for k from 1 to n do   # assume we have known the k-1 column of L,then 
    			d[k]:=A[k,k];        # compute the k column.
    			for i from 1 to k-1 do
      				d[k]:=d[k]-d[i]*L[k,i]^2;
    			end do;
    			if d[k]<>0 then
      				L[k,k]:=1;
      				for i from k+1 to n do 
        				L[i,k]:=A[i,k];
        				for j from 1 to k-1 do
          					L[i,k]:=L[i,k]-d[j]*L[i,j]*L[k,j];
        				end do;
       	 				L[i,k]:=L[i,k]/d[k];
      				end do;
    			end if;
  		end do;
  		return L,d;
	end proc:  




	#====================================================================================== 
	# Discription: 
	#    Return the objective function for plane search and the feasible rectangle of p,q 
	# Inputs:  
	#    X       - a given strictly feasible point 
	#    Z       - a given strictly dual feasible point 
	#    c,F0,L  - the primal SDP problem data 
	#    dX      - a given suitable direction 
	#    dZ      - a given suitable dual feasible direction 
	# Outputs: 
	#    f       - the objective function for plane search 
	#    rect_pq - the feasible rectangle of p,q 
	#====================================================================================== 
	Expresstion_of_PlaneSearch:=proc(X,Z,F0,c,L::list,dX,dZ,v) 
   		local i,j,m,n,c1,c2,F,Tr_FZ,dF,L1,L2,f,h,rect_pq,M_L,Inv_F,Inv_Z,p,q,vars,s;   
   		F:=F0; 
   		dF:=Matrix(ColumnDimension(Z)); 
   		m := nops(c); 
   		for i from 1 to m do 
      			F := F+X[i]*L[i]; 
      			dF := dF+dX[i]*L[i]; 
   		end do; 
   		F := evalf(F); 
   		Tr_FZ:=0;
   		for i from 1 to ColumnDimension(F) do
     			for j from 1 to ColumnDimension(F) do
       				Tr_FZ:=Tr_FZ+F[i,j]*Z[j,i];
     			end do;
   		end do;
   		c1 := Vector(c).Vector(dX)/Tr_FZ; 
   		s:={};
   		for i from 1 to ColumnDimension(F0) do
     			for j from 1 to ColumnDimension(F0) do
       				if F0[i,j]<>0 then
         				s:=s union {[i,j]};
       				end if;
     			end do;
   		end do;
   		c2:=0;
   		for i from 1 to nops(s) do
     			c2:=c2+F0[s[i][1],s[i][2]]*dZ[s[i][2],s[i][1]];
   		end do;
   		c2 := c2/Tr_FZ; 
   		printf("begin to compute the Eigenvalues \n"); 
   		printf("begin to compute M_L \n"); 
   		Inv_F:=MatrixInverse(F);
   		M_L :=evalf(Inv_F.dF);
   		printf("end to compute M_L and begin to compute the Eigenvalues\n"); 
   		L1 := Eigenvalues(M_L,output='list'); # elements are complex because of machine erro 
   		L1:=[seq(Re(L1[i]),i=1..nops(L1))]; # get real part 
   		Inv_Z:=MatrixInverse(Z);
   		L2 := Eigenvalues(evalf(Inv_Z.dZ),output='list'); # elements are complex because of machine erro 
   		L2:=[seq(Re(L2[i]),i=1..nops(L2))]; # get real part 
   		printf("end to compute the Eigenvalues \n"); 
   		n := nops(L1); 
   		f := (n+v*sqrt(n))*log(1+c1*p+c2*q); 
   		for i from 1 to n do  
      			f := f-log(1+p*L1[i])-log(1+q*L2[i]); 
   		end do; 
   		vars:=[p,q];
   		return f,c1,c2,L1,L2,vars; 
	end proc:




	#====================================================================================== 
	# Discription: 
	#    Given a list of real values, this function return the minimal positive value and  
	#    the maximal negative value 
	# Input: 
	#    L   - a list of real values 
	# Output: 
	#    min - the minimal positive value 
	#    max - the maximal negative value 
	#====================================================================================== 
	Feasible_Rectangle2:=proc(L::list) 
   		local i,min,max,L_negative,L_positive; 
   		L_negative:=[];  L_positive:=[]; 
   		for i from 1 to nops(L) do  
      			if L[i] < 0 then 
         			L_negative := [op(L_negative),L[i]]; 
      			else 
         			L_positive := [op(L_positive),L[i]]; 
      			end if; 
   		end do; 
   		min := L_negative[1]; max := L_positive[1]; 
   		for i from 2 to nops(L_negative) do  
      			if min < L_negative[i] then 
         			min := L_negative[i]; 
      			end if; 
   		end do; 
   		for i from 2 to nops(L_positive) do  
      			if max > L_positive[i] then 
         			max := L_positive[i]; 
      			end if; 
   		end do; 
   		return min, max; 
	end proc:




	#============================================================================== 
	# Discription: 
	#    Return the feasible rectangle of p and q 
	# Inputs: 
	#    list_u          - the list of eigenvalues u_1..u_n 
	#    list_v          - the list of eigenvalues v_1..v_n 
	# Outputs: 
	#    [p_min,p_max]   - 
	#    [q_min,q_max]   _ 
	#============================================================================== 
	Feasible_Rectangle:=proc(list_u,list_v) 
   		local i,n,u_min,u_max,v_min,v_max,p_min,p_max,q_min,q_max; 
   		if nops(list_u)<>nops(list_v) then 
      			error("Invalid input, the dimension of list_u and list_v must be equal"); 
   		end if; 
   		n:=nops(list_u); 
   		u_min:=list_u[1];u_max:=list_u[1]; 
   		v_min:=list_v[1];v_max:=list_v[1]; 
   		for i from 2 to n do 
      			if list_u[i]<u_min then  
         			u_min:=list_u[i];   
      			end if; 
      			if list_u[i]>u_max then  
         			u_max:=list_u[i];   
      			end if; 
      			if list_v[i]<v_min then  
         			v_min:=list_v[i];   
      			end if; 
      			if list_v[i]>v_max then  
         			v_max:=list_v[i];   
      			end if; 
   		end do; 
   		if u_max>0 then p_min:=-1/u_max; else p_min:=-infty end if; 
   		if u_min<0 then p_max:=-1/u_min; else p_max:=infty end if; 
   		if v_max>0 then q_min:=-1/v_max; else q_min:=-infty end if; 
   		if v_min<0 then q_max:=-1/v_min; else q_max:=infty end if; 
   		return [p_min,p_max,q_min,q_max]; 
	end proc: 




	#================================================================================================== 
	# Discription 
	#    Return the root of f(x) lying between the interval L using the bisection method 
	# Inputs: 
	#    f    - the function 
	#    L    - L[1],L[2] are the left and right end point of the interval respectively 
	#    tol  - the tolerence 
	# Output: 
	#    mid  - the root of f(x) within the error tolerence 
	#==================================================================================================  
	RootBisection :=proc(f,L::list,tol) 
   		local min,max,mid,r; 
   		if nops(L)<>2 then  
      			error("invalid input, dimension of L must be two, min,max, the interval"); 
   		end if; 
   		if L[2]<L[1] then  
      			error("invalid input, L[1] must be smaller than or equal to L[2] "); 
   		end if; 
   		if nops(indets(f))<>1 then
     			error("Invalid input, f must be a univariable polynomial");
   		end if;
   		r:=op(indets(f));
   		min := L[1]; max := L[2]; 
   		mid := (min+max)/2; 
   		while abs(subs(r=mid,f))>tol do 
     			if subs(r=mid, f)*subs(r=min, f)>0 then 
        			min := mid; 
     			else  
        			max := mid; 
     			end if; 
     			mid := (min+max)/2; 
   		end do; 
   		return mid; 
	end proc:




  	#=========================================================================================== 
	# Discription:  
	#    Return the length of each step in plane search. This function involves computing minimal 
	#    value of f. When the hessen matrix at an iteration point is positive definite, we use  
	#    the Newton direction, otherwise we use the steepest decent direction. Then we compute the  
	#    one dimnension search of r, in x_k+1=x_k+r*dx, to get the update point. 
	# Inputs: 
	#    f         - the optimization function for computing the length of step 
	#    v0        - the initial point, usually we choose [0,0] 
	#    tol       - the tolerence in scaling the value of f between two steps 
	#    max_iter  - the maximal number of iterations 
	#    c1,c2     - parameters in f which are computed in Expresstion_of_PlaneSearch 
	#    L1,L2     - eigenvaluies u1..un, v1..vn 
	# Outputs: 
	#    v[1],v[2] - the length of step respectivly for dX, dZ 
	#=========================================================================================== 
	###################################################################################################################### 	 
	StepsLength_PlaneSearch:=proc(f,v0,tol,max_iter,c1,c2,L1,L2,vars) 
      		local i,j,Jac_f,Hess_f,A,b,list_var,v,dv,dValue_f,value_f,L,min,s1,s2,s,L_r,rect_r,rr,r,p,q; 
      		p:=vars[1]; q:=vars[2];
      		list_var := vars; # the variables in f 
      		s1:=[seq(1+p*L1[i],i=1..nops(L1))];  
      		s2:=[seq(1+q*L2[i],i=1..nops(L2))];  
      		s:=[1+c1*p+c2*q,op(s1),op(s2)];  # get the expressions of p and q in the log function  
      		Jac_f := VectorCalculus[Jacobian]([f],list_var); 
      		Hess_f := VectorCalculus[Hessian](f,list_var); 
      		v := convert(v0,Vector); dValue_f := 1; 
      		A :=eval(Hess_f,[p=v[1],q=v[2]]); 
      		b :=convert(eval(Jac_f,[p=v[1],q=v[2]]),Vector); 
      		value_f := eval(f,[p=v[1],q=v[2]]); 
      		for i from 1 to max_iter do 
        		# printf("the %d iteration begins",i);
        		# print(`dValue,tol`,dValue_f,tol);
         		if abs(dValue_f) > tol then 
             			A := evalf(eval(Hess_f,[p=v[1],q=v[2]]));   
             			L := evalf(Eigenvalues(A,output='list'));  
             			L := [seq(Re(L[i]),i=1..nops(L))]; # the eigenvalus of Hessen matrix A 
             			min:=L[1];  
             			for j from 2 to nops(L) do  # get the minimal value of L 
                			if L[j]<min then 
                   				min:=L[j]; 
                			end if; 
             			end do;  
             			if min<=0 then   # Hessen matrix is not positive definite, Newton direction is invalide, we use steepest decent direction 
             				#   printf("the Hessen matrix of this step is not positive, and steepest decent direction is adopted.\n");     
                			dv := eval(Jac_f,[p=v[1],q=v[2]]); # the steepest decent direction 
                			dv:=convert(dv,Vector);   
       					#  print(`dv`,evalf(dv));
             			else 
             				# printf("the Hessen matrix of this step is positive, and Newton direction is adopted.\n");     
                 			b :=convert(eval(Jac_f,[p=v[1],q=v[2]]),Vector);  
                 			dv :=LinearAlgebra[LinearSolve](A,b);    # the Newton direction 
             			end if;   # 1 dimension search, to get the optimal r 
             			s:=subs(p=(v-r*dv)[1],q=(v-r*dv)[2],s); 
             			s:=[seq(expand(evalf(s[i])),i=1..nops(s))];# get the expressions of r in the log function  
             			L_r:=[seq(-(s[i]-(coeff(s[i],r))*r)/coeff(s[i],r),i=1..nops(s))];  
             			rect_r:=Feasible_Rectangle2(L_r);  # get the feasible interval of r   
             			rr := RootBisection(diff(evalf(eval(f,[p=(v-r*dv)[1],q=(v-r*dv)[2]])),r),[0,rect_r[2]],10^(-10));     
             			v := evalf(v-rr*dv);      # update                     
             			s:=[1+c1*p+c2*q,op(s1),op(s2)];     
             			dValue_f := evalf(value_f-eval(f,[p=v[1],q=v[2]]));  
             			print(`dValue_f= `,dValue_f); 
             			value_f := evalf(eval(f,[p=v[1],q=v[2]])); 
         		else  
             			break; 
         		end if; 
      			# printf("the %d iteration ends", i);
       		end do; 
       		return v[1],v[2]; 
	end proc: 




	#=========================================================================================================                    
	# Discription: 
	#    Return the solution of SDP problem using method1 for computing dX and dZ 
	# Inpuuts: 
	#    c        - paramenters in primal problem 
	#    X0,Z0    - the initial point which are strickly positive 
	#    L        - list, F_0,F_1,..F_m 
	#    v        - parameter in potential function 
	#    max_iter - the maximal number of iterations 
	# Outputs: 
	#    X, Z     - the solution of SDP problem 
	#======================================================================================================== 
	SolveSDP:=proc(c,x0,Z0::Matrix,L::list,max_iter::integer,{v::integer:=5}) 
    		local i,n,m,rho,F0,F,S,D,X,Z,dX,dZ,g,p_,q_,value_potential,gap,c1,c2,L1,L2,num_iter,l,vars; 
    		F0:=L[1]; X:=x0; Z:=Z0;
    		m:=nops(L)-1; num_iter := 1; 
    		for i from 1 to max_iter  do 
     			 printf("The %d th iteration begins", num_iter); 
      			F:=evalf(add(X[i-1]*L[i],i=2..m+1)+L[1]); # get F(x)  
      			n:=ColumnDimension(F);
      			rho:=(n+v*sqrt(n))/(add(X[i]*c[i],i=1..m)+Trace(L[1].Z));# get the rho, S, D for Method 1  
      			D:=evalf(rho*evalf(F.Z.F)-F); 
      			S:=F; 
      			printf("begin to compute the dX and dZ \n");     # return S,D,[seq(L[i],i=2..m+1)];
      			dX,dZ:=Solution_of_Directions(S,D,[seq(L[i],i=2..m+1)]); #return dX,dZ;#if i=max_iter then return dX,dZ; fi; 
      			printf("end to compute the dX and dZ \n"); 
      			printf("begin to compute the Expresstion_of_PlaneSearch \n"); 
      			g,c1,c2,L1,L2,vars:=Expresstion_of_PlaneSearch(X,Z,L[1],c,[seq(L[i],i=2..m+1)],dX,dZ,v);
      			printf("end to compute the Expresstion_of_PlaneSearch \n"); 
      			printf("begin to compute the StepsLength_PlaneSearch \n");
      			l:=Feasible_Rectangle(L1,L2);
      			print(l);
      			if l[2]=infty then l[2]:=0; end if;
      			p_, q_:= StepsLength_PlaneSearch(g,[l[2]/2,0],0.000001,100,c1,c2,L1,L2,vars);       
     			printf("end to compute the StepsLength_PlaneSearch \n"); 
      			X:=evalf(X+p_*[seq(dX[i],i=1..m)]); 
      			Z:=evalf(Z+q_*dZ); 
      			gap:= evalf(PotentialFunction(X,Z,L[1],[seq(L[i],i=2..m+1)],v)); 
      			print(`gap= `, gap); 
     			printf("The %d th iteration ends \n", num_iter); 
      			num_iter := num_iter + 1; 
    		end do; 
    		return X,Z; 
	 end proc: 




	#=========================================================================================== 
	# Discription 
	#    Transform the original primal-dual SDP problem to a new one in case that neither  
	#    a strickly feasible x nor a strickly feasible Z is known. More details can be found 
	#    in the paper(Semidefinite Programming, p40). Return the new data of the new problem 
	#    and the strikly feasible x and Z  
	# Inputs: 
	#    L,c          - the F_0..F_m, c of the original primal problem
	#    m1,m2        - the initial M1,M2 
	# Outputs: 
	#    L_new,c_new  - the new F_0..F_m, c of the new primal problem 
	#    X_new,Z_new  - initial strickly feasible primal and dual points for the new problem 
	#========================================================================================== 
	###################################################################################### 
	BigM:=proc(L::list,c::list,m1,m2) 
   		local i,j,k,m,n,M,M1,M2,s,var,A,b,X,List,min,z,Z_new,t,L_new,X_new,c_new; 
   		var := []; L_new := []; 
   		m := nops(L)-1; 
   		n := ColumnDimension(L[1]); 
   		Z_new := Matrix(n+2);  
   		for i from 2 to m+1 do 
     			s := [];  
     			for j from 1 to n do    
       				for k from 1 to n do    
         				if (L[i])[j,k]=1 then       
           					s:=[op(s),[j,k]];       
         				end if;      
       				end do;  
     			end do;    
     			for j from 1 to nops(s) do    
       				Z_new[(s[j])[1],(s[j])[2]]:=c[i-1]/nops(s);  
     			end do;  
   		end do;       
   		List := evalf(Eigenvalues(Z_new[1..n,1..n],output='list'));  
   		List := [seq(Re(List[i]),i=1..nops(List))]; 
   		min:=List[1];  
   		for j from 2 to nops(List) do  # get the minimal value of L 
      			if List[j]<min then 
         			min:=List[j]; 
      			end if; 
   		end do;  
   		if min < 0 then 
      			z := -min+1; 
   		else 
      			z := 1; 
   		end if; 
   		Z_new[1..n,1..n] := Z_new[1..n,1..n]+z*IdentityMatrix(n); 
   		M1 := Trace(Z_new[1..n,1..n])+1; 
   		M1 := 10*M1;  
   		if M1<m1 then   
     			M1:=m1; 
   		end if;
   		Z_new[n+1,n+1] := z; 
   		Z_new[n+2,n+2] := M1-Trace(Z_new[1..n,1..n]); 
   		List := evalf(Eigenvalues(L[1],output='list')); #let x_0 =0; 
   		List := [seq(Re(List[i]),i=1..nops(List))]; 
   		min:=List[1];  
   		for j from 2 to nops(List) do  # get the minimal value of L 
      			if List[j]<min then 
         			min:=List[j]; 
      			end if; 
   		end do;  
   		if min < 0 then 
      			t := -min+1; 
   		else 
      			t := 1; 
   		end if; 
   		M2 := Trace(L[1])+100;   
   		if m2=0 then   
     			M2 := 10*M2;   
   		else   
     			M2:=m2;  
   		end if; 
   		print(`M1 is`, M1);   
   		print(`M2 is`,M2); 
   		L_new:=[];
   		for i from 1 to m+1 do   
     			L_new:=[op(L_new),Matrix(n+2)];   
     			for j from 1 to n do  
       				for k from 1 to n do  
         				(L_new[i])[j,k]:=(L[i])[j,k];
      				end do;   
     			end do;  
   		end do; 
   		(L_new[1])[n+1,n+1] := M2-Trace(L[1]); 
   		for i from 2 to m+1 do 
      			(L_new[i])[n+1,n+1] := -Trace(L[i]); 
   		end do; 
   		M := Matrix(n+2); 
   		for i from 1 to n do 
      			M[i,i] := 1; 
   		end do; 
   		M[n+2,n+2] := 1;  
   		L_new := [op(L_new), M];  
   		X_new := [seq(0,i=1..m)]; 
   		X_new := [op(X_new), t]; 
   		c_new := [seq(0,i=1..nops(c)+1)];
   		for i from 1 to nops(c) do 
      			c_new:=subsop(i=c[i],c_new);
   		end do;
   		c_new:=subsop(nops(c)+1=M1,c_new);
   		return L_new,c_new,X_new,Z_new;  
	end proc: 


end module;





	
ConvexHull:=module()
	description "package for computing the convex hull of a given set of points";
	export 
		kth_subset,
		allSubsets,	
		inConvexHull;
	local
        convexHull,
		creatHyperplane,
		creatSimplex,
		partition,
		furthestPoint,
		visibleSet,
		creatFacets,
		simplicialConvexHull,
		certification,
		nChoosek,
		getR,
		affineDims,
		affineTrans;
	option 
		package;
	


 	#===================================================================   
	#Description  
	#  Given d points which are independant in d dimensional vector space,  
	#  return the representation of the hyperplane containing the points.   
	#Input   
	#  S     --- the d points  
	#Output  
	#  x     --- the representation of the hyperpalne  
	#            (x[1..d].p=x[d+1])  
	#===================================================================  
	creatHyperplane:=proc(S::set)  
  		local i,j,d,A,x;  
  		d:=nops(S);  
  		A:=Matrix(d,d+1);  
  		for i from 1 to d do  
    			A[i,1..d]:=Vector[row](S[i]);  
    			A[i,d+1]:=-1;  
  		end do;  
  		x:=LinearSolve(A,Vector(d));   
  		for i from 1 to d do  
    			if x[i]<>0 then   
      				x:=x/x[i]; #return v[i];  
      				break;  
    			end if;  
  		end do;   
  		x:=x/(add(x[i]^2,i=1..d))^(1/2);   
  		return x;  
	end proc:





	#====================================================================   
	#Description   
	#  Given d+1 points which are affine independent, return   
	#  the simplex of these d+1 points. We assume that the point  
	#  are represented by a list which contains its coordinates   
	#Input  
	#  S    --- the d+1 points  
	#Output  
	#  F_V  --- the vertices of the d+1 facets of the simplex  
	#  H    --- the hyperplane equations of the d+1 facets of the simplex  
	#           H[i][1..d] is the norm vector whose 2-norm is 1,  
	#           H[i][1..d].x=H[i][d+1] is the hyperplane equation  
	#===================================================================  
	creatSimplex:=proc(S::set)  
  		local i,j,d,x,H,A,b,F_V,F,p;  
  		d:=nops(S)-1;  
  		for i from 1 to d+1 do  
    			if nops(S[i])<>d then   
      				error("Invalid inputs");  
    			end if;  
  		end do;  
  		H:=[]; F_V:=[];  
  		for i from 1 to d+1 do   
    		F:=S; p:=S[i];#return s[i];  
    		F:=F minus {p};#return s[i];  
    		x:=creatHyperplane(F);  
    		if evalf(Vector([op(p),-1]).x)>0 then  
      			x:=-x;  
    		end if;  
    		H:=[op(H),x]; F_V:=[op(F_V),F];  
  		end do;  
  		return F_V,H;  
	end proc:    
    

	#================================================================== 
	#Description 
	#  Given a set of points and a set of facets, this function assigns 
	#  the points to each facet's outside set. 
	#Inputs 
	#  P      --- the unassigned points 
	#  F_V,H  --- the representation of the facets(its vertices and hyperplanes) 
	#Outputs 
	#  S      --- each element of S is the outside set of the relevant facet   
	#================================================================== 
	partition:=proc(P,F_V,H) 
  		local i,j,d,d_p,S,s; 
  		d:=nops(F_V); 
  		if d-nops(H)<>0 then 
    			error("The inputs are invalid, the number of facets and hyperplanes must be equal"); 
  		end if; 
  		d_p:=nops(P); 
  		S:=[]; 
  		for i from 1 to d do 
    			s:={}; 
    			for j from 1 to d_p do 
      				if evalf(Vector[row]([op(P[j]),-1]).H[i])>0 then 
        				s:=s union {P[j]}; 
     				end if; 
    			end do; 
    			S:=[op(S),s];
  		end do; 
  		return S; 
	end proc: 
 




	#======================================================================= 
	#Description  
	#  Given a set of points and a hyperplane, return the furthest point to  
	#  the plane. We assume the hyperplane is represented by a vector H and the  
	#  hyperplane's equation is H[1..Dim-1].x=H[Dim] 
	#inputs 
	#  P         --- a set of points 
	#  H         --- the equation of the hyperplane, represented by a vector 
	#Output 
	#  P[index]  --- the furthest point to the hyperplane 
	#====================================================================== 
	furthestPoint:=proc(P,H) 
  		local i,max,ind,d; 
  		d:=nops(P); 
  		max:=evalf(Vector[row]([op(P[1]),-1]).H); 
  		ind:=1; 
  		for i from 1 to d do 
    			if max < evalf(Vector[row]([op(P[i]),-1]).H) then 
      				max:=evalf(Vector[row]([op(P[i]),-1]).H); 
      				ind:=i; 
    			end if; 
  		end do; 
  		return P[ind]; 
	end proc: 
 
 
	#================================================================================ 
	#Description 
	#  Given a point P and a set of facets. The point is in the outside set of one  
	#  facets. This function finds the visible set(denoted V) of P, and the boundary 
	#  of V.  
	#Inputs  
	#  P       --- a point 
	#  F_V,H   --- the representation of the facets(its vertices and hyperplanes) 
	#  adjM    --- the adjacency matrix of the facets 
	#  index   --- the P is in the outside set of the index-th facet 
	#Outputs 
	#  VS      --- the visible set to P 
	#  ridges  --- the set of horizon ridges 
	#================================================================================ 
	visibleSet:=proc(P,F_V,H,adjM,ind) 
		local i,d,VS,UN,ridges,s,VF,index2,s1,s2; #VF:= visited facet; UN:=unvisited neighbors; VS:= visible set 
  		d:=nops(F_V); 
  		if d<>nops(H) or d<>RowDimension(adjM) or d<>ColumnDimension(adjM) then 
    			error("The inputs are invalid, Dimensions must be equal"); 
  		end if; 
  		VS:={ind}; UN:={}; ridges:={};  
  		VF:=[seq(1,i=1..d)];
  		VF:=subsop(ind=0,VF);
  		for i from 1 to d do   #initialize the unvisited neighbors 
    			if adjM[ind,i]=1 and i<>ind then 
      				UN:=UN union {i}; 
    			end if;     
  		end do; #return UN; 
  		for i from 1 to nops(F_V[ind]) do  # initialize the horizon ridges 
    			s1:=F_V[ind] minus {F_V[ind][i]}; 
    			ridges:=ridges union {s1}; 
  		end do; #return UN,ridges; 
  		while UN<>{} do  
    			index2:=UN[1]; VF:=subsop([index2]=0,VF); #mark the visited facet 
    			if evalf(Vector[row]([op(P),-1]).H[index2])>0 then #determine visible set 
      				VS:=VS union {index2};  
      				s1:=F_V[index2]; s2:={};  # add the new ridges 
      				for i from 1 to nops(F_V[index2]) do  
        				s2:=s2 union {s1 minus {s1[i]}}; 
      				end do; 
      				s1:=ridges minus ((ridges union s2) minus s2);  
      				ridges:= (ridges union s2) minus s1; 
      				for i from 1 to d do    
        				if adjM[index2,i]=1 and VF[i]=1 then 
          					UN:=UN union {i}; # add new unvisited neighbors 
        				end if; 
       				end do; 
     			end if;   
    			UN:=UN minus {index2}; # delete the visited neighbors 
  		end do; 
  		return VS,ridges; 
	end proc: 
   
   
	#============================================================================= 
	#Description 
	#  For each ridge R, this function creats a new facet from R and p and links the  
	#  new facet to its neighbors 
	#Inputs 
	#  P               --- the furthest point to be considered 
	#  VS              --- the visible set to point p 
	#  ridges          --- the horizon ridges 
	#  adjM            --- the adjacency matrix 
	#  F_V,H           --- the representation of the facets  
	#Outputs 
	#  adjM_new        --- the new adjacency matrix after deleting the visible set 
	#  F_V_new,H_new   --- the representation of the facets after deleting the VS   
	#============================================================================= 
	creatFacets:=proc(P,VS,ridges,adjM,F_V,H) 
  		local i,j,d,m,n,x,s,p,adjM_temp,adjM_new,F_V2,H2,F_V_new,H_new,M; 
  		d:=nops(F_V); n:=nops(ridges); m:=nops(VS); 
  		if d-nops(H)<>0 or d-RowDimension(adjM)<>0 or d-ColumnDimension(adjM)<>0 then 
    			error("Invalid inputs, the nember of the facets must be equal"); 
  		end if; 
  		if m=0 then 
    			error("visible set must have at least one element"); 
  		end if; 
  		p:=add(F_V[VS[1]][i],i=1..nops(F_V[VS[1]]))/nops(F_V[VS[1]]);# used to determine th direction of norm vector 
  		F_V2:=F_V; H2:=H;  
  		for i from 1 to n do  # creat new facet 
    			s:=ridges[i] union {P};  
    			x:=creatHyperplane(s); 
    			if evalf(Vector[row]([op(p),-1]).x)>0 then  
      				x:=-x; 
    			end if; 
    			F_V2:=[op(F_V2),s]; H2:=[op(H2),x];  
  		end do;   
  		adjM_temp:=Matrix(d+n);    #creat new adjacency matrix 
  		adjM_temp[1..d,1..d]:=adjM;  
  		for i from 1 to d+n do 
    			for j from d+1 to d+n do 
      				s:=F_V2[j] minus ((F_V2[i] union F_V2[j]) minus F_V2[i]); 
      				if nops(s)=nops(F_V2[1])-1 then 
        				adjM_temp[i,j]:=1; 
      				end if; 
    			end do; 
  		end do; 
  		for i from d+1 to d+n do 
    			for j from 1 to d do 
      				adjM_temp[i,j]:=adjM_temp[j,i]; 
     			end do; 
  		end do; 
  		F_V_new:=[]; H_new:=[];  
  		adjM_new:=Matrix(d+n-m); M:=Matrix(d+n,d+n-m); 
  		#delete the facets in VS from F_V, H and M 
  		for i from 1 to VS[1]-1 do 
    			F_V_new:=[op(F_V_new),F_V2[i]]; 
    			H_new:=[op(H_new),H2[i]]; 
    			M[1..d+n,i]:=adjM_temp[1..d+n,i]; 
  		end do; 
  		for i from 1 to m-1 do 
    			for j from VS[i]+1 to VS[i+1]-1 do 
      				F_V_new:=[op(F_V_new),F_V2[j]]; 
      				H_new:=[op(H_new),H2[j]]; 
      				M[1..d+n,j-i]:=adjM_temp[1..d+n,j]; 
    			end do; 
 	 	end do;   
  		for i from VS[m]+1 to d+n do 
    			F_V_new:=[op(F_V_new),F_V2[i]]; 
    			H_new:=[op(H_new),H2[i]];      
    			M[1..d+n,i-m]:=adjM_temp[1..d+n,i];    
  		end do;    
  		for i from 1 to VS[1]-1 do 
    			adjM_new[i,1..d+n-m]:=M[i,1..d+n-m]; 
  		end do; 
  		for i from 1 to m-1 do 
    			for j from VS[i]+1 to VS[i+1]-1 do 
      				adjM_new[j-i,1..d+n-m]:=M[j,1..d+n-m]; 
    			end do; 
  		end do; 
  		for i from VS[m]+1 to d+n do 
    			adjM_new[i-m,1..d+n-m]:=M[i,1..d+n-m]; 
  		end do; 
  		return adjM_new,F_V_new,H_new; 
	end proc: 
 

	#================================================================================================== 
	#Description 
	#  In n dimensional vector space, given a set of points which are in general position(so it has n+1  
	#  affine independant points), this function computes the convex hull. So the convex hull is a  
	#  simplicial complex. 
	#Inputs 
	#  S       ---the set of lists which represent the points 
	#  indS    ---the n+1 independant points 
	#Outputs 
	#  F_V     ---a list of lists of lists which are the vertices of the facets 
	#  H       ---a list of normal vectors of the facets having the same order within F_V   
	#================================================================================================== 
	simplicialConvexHull:=proc(S,indS) 
  		local i,j,d,F,F_V,H,s,outS,P,VS,ridges,adjM,ind,iter; 
  		F_V,H:=creatSimplex(indS);  
  		d:=nops(F_V); 
  		adjM:=Matrix(d,(i,j)->1); 
  		ind:=0; iter:=0; 
  		outS:=partition(S,F_V,H);  
  		for i from 1 to d do  
    			if outS[i]<>{} then  
      				ind:=i;  
      				break;  
    			end if;  
  		end do;  
  		while ind>0 do    
    			iter:=iter+1;  
    			P:=furthestPoint(outS[ind],H[ind]);  #print(ind,P); 
    			VS,ridges:=visibleSet(P,F_V,H,adjM,ind);
    			adjM,F_V,H:=creatFacets(P,VS,ridges,adjM,F_V,H); #print(H);#return adjM,F_V,H; 
    			outS:=partition(S,F_V,H);  
    			d:=nops(outS);  ind:=0; 
    			for i from 1 to d do   
      				if  outS[i]<>{} then  
        				ind:=i;  
        				break;  
      				end if;  
    			end do;  
  		end do; 
  		return F_V,H; 
	end proc: 
 
	#=============================================================================================== 
	#Description 
	#  certify the convex hull 
	#=============================================================================================== 
	certification:=proc(S,H) 
  		local i,j,m,n,s,s_; 
  		m:=nops(S);  n:=nops(H); 
  		s:={};   
  		for i from 1 to m do 
    			s_:={S[i]}; 
    			for j from 1 to n do 
      				if evalf(Vector[row]([op(S[i]),-1]).H[j])>0 then 
        				s_:=s_ union {j}; 
      				end if; 
    			end do; 
    			if s_<>{S[i]} then 
      				s:=s union {s_}; 
    			end if; 
  		end do; 
  		if s={} then 
    			#printf("The convex hull is computed right \n"); 
  		else 
    			printf("The convex hull is computed wrong, and the points which are not in the convex hull is following \n"); 
    			return s; 
  		end if; 
	end proc: 
 
 
	#=========================================================================================================== 
	#Description 
	#  Given integer n and k(k<=n), return all the subsets of [1...n], which have k elements. 
	#Inputs 
	#  n, k   --- described above 
	#Output 
	#  S      --- a list of lists which are subsets having k elements 
	#=========================================================================================================== 
	########################################################################################

	nChoosek:=proc(n,k)
  		return n!/(k!*(n-k)!);
	end proc:

	kth_subset:=proc(min,max,k)
  		local i,j,M,m,n,num,L;
  		n:=max-min+1;
  		if k>n then
   			error("Invalid inputs, k is larger than max-min+1");
  		end if;
  		m:=nChoosek(n,k);
  		M:=Matrix(m,k);
  		num:=[seq(0,i=1..k)];
  		for i from 1 to k do
    			m:=nChoosek(n-i,k-i);
    			M[1..m,i]:=min+i-1;
    			num[i]:=m;
  		end do;
  		for i from k by -1 to 1 do
    			for j from min+i to max-(k-i) do
      				m:=nChoosek(max-j,k-i);
      				M[num[i]+1..num[i]+m,i]:=j;
      				M[num[i]+1..num[i]+m,i+1..k]:=M[num[i]-m+1..num[i],i+1..k];
      				num[i]:=num[i]+m;
    			end do;
  		end do; 
  		L:=[];
  		for i from 1 to RowDimension(M) do
    			L:=[op(L),convert(M[i,1..k],list)];
  		end do;
  		return L;
	end proc:
 
	#===================================================================================== 
	#Description 
	#  Return all the subsets of [1...n], except the empty subset. 
	#===================================================================================== 
	allSubsets:=proc(n)   
  		local i,L; 
 		L:=[]; 
  		for i from 1 to n do 
    			L:=[op(L),op(kth_subset(1,n,i))]; 
  		end do; 
  		return L; 
	end proc: 
 
	#======================================================================================
	#Description
	#  Return k different random numbers from the range 1..n 
	#======================================================================================
	getR:=proc(n,k)
  		local i,j,L,roll,flag;
  		if type(n,integer)=false or type(k,integer)=false or n<=0 or k<0 then
    			error("Invalid inputs, n must be a positive integer and k be nonnegative integer");
  		end if;
  		if k>n then
    			error("Invalid inputs, k must be no larger than n");
  		end if;
  		L:=[];
  		roll:=rand(1..n);
  		for i from 1 to k do
    			flag:=0;
    			L:=[op(L),roll()];
    			while flag=0 do
      				flag:=1;
      				for j from 1 to i-1 do
        				if L[i]=L[j] then
          					L:=subsop(i=roll(),L);
	      					flag:=0;
          					break;
        				end if;
      				end do;
    			end do;
  		end do;
  		return L;
	end proc:

  
	#====================================================================================== 
	#Description 
	#  Given a set of points in n dimensional vector space, return the affine dimension of  
	#  the set and a affine independant points. 
	#Input 
	#  S    --- a set of lists which represent the points  
	#Outputs 
	#  r    --- the affine dimension of the set 
	#  indS --- a set of r-affinely independant points  
	#====================================================================================== 
 	affineDims:=proc(S::set) #detect dimension from high to low using random chosen points
  		local i,j,k,m,n,N,s,r,M,indS,ind,R,s2; 
  		m:=nops(S); n:=nops(S[1]); 
  		for i from 2 to m do 
    			if n<>nops(S[i]) then 
      				error("Invalid inputs, all the points must have the same dimension"); 
    			end if; 
  		end do; 
  		s:=[op(S)]; 
		M:=Matrix(n);
		r:=0; i:=2; j:=1;
		indS:={s[1]};
		while i<=nops(s) do
			M[j,1..n]:=Vector[row](s[i]-s[1]); 
			if Rank(M)=r+1 then
				r:=r+1;	j:=j+1;
				indS:=indS union {s[i]};
				if r=n then
					break;
				fi;
				next;
			else
				M[j,1..n]:=0; 
			fi;
			i:=i+1;
		od;	 
  		return r, indS;
	end proc:

 
	#=========================================================================================== 
	#Description 
	#  Given a set points with r-affine dimension, this fuction affinely transformate thes points 
	#  into a r dimensional vector space. We make r vectors from the r+1 points in indS as the  
	#  basis of S and denote A as the basis matrix. If p=A.x, then f(A):=IdentityMatrix(r).x. So  
	#  f is a affine transformation. With the matrix A and prebase, we can tranformate a point in  
	#  S to r dimensional space and recover it back to n dimensional spsce. 
	#Inputs 
	#  S           --- a set of points with r-affine dimension 
	#  indS        --- a subset of r-affinely independant points 
	#  r           --- affine dimension 
	#Outputs   
	#  TRS,TrIndS  --- the points after transformation in S and indS, respectively 
	#  A           --- the basis matrix of S with the r vectors in indS as the basis 
	#  preBase     --- the r vectors composed of the r+1 points in indS 
	#=========================================================================================== 
	affineTrans:=proc(S,indS,r) 
  		local i,j,A,b,x,M,n,TrS,TrIndS,preBase; 
  		if nops(indS)<>r+1 then  
    			error("Invalid inputs, the independant set must have r+1 elements"); 
  		end if; 
  		if nops(S)=0 or nops(indS)=0 then 
    			error("Invalid inputs, the inputs sets must be nonempty"); 
  		end if; 
  		n:=nops(S[1]); 
  		TrS:=[op(S)]; TrIndS:=[op(indS)];  
  		M:=IdentityMatrix(r); 
  		for i from 1 to nops(TrS) do 
    			A:=Matrix(n,r); 
    			for j from 1 to r do 
      				A[1..n,j]:=Vector(TrIndS[j+1]-TrIndS[1]); 
    			end do; 
    			x:=LinearSolve(A,Vector(TrS[i]-TrIndS[1])); 
    			TrS:=subsop(i=convert(M.x,list),TrS); 
  		end do; 
  		preBase:=TrIndS; 
  		TrIndS:=subsop(1=[seq(0,i=1..r)],TrIndS); 
  		for i from 2 to r+1 do 
    			TrIndS:=subsop(i=convert(M[1..r,i-1],list),TrIndS); 
  		end do; 
  		return TrS,TrIndS,A,preBase; 
	end proc: 
    
 
	#========================================================================================= 
	#Description 
	#  Given a set points which are n dimensional, return the its convex hull. If these points 
	#  sre r-affine dimensional, with r<n, then we make affine transformation on the points into 
	#  r dimensional vector space, and computed the convex hull of the transfated points. 
	#========================================================================================= 
	convexHull:=proc(S::set) 
  		local i,j,m,n,r,indS,end1,end2,ind,TrS,TrIndS,A,preBase,F_V,H; 
  		m:=nops(S); 
 	 	if m=0 then  
    			error("Invalid input, the set must be nonempty set"); 
  		end if; 
  		if m=1 then  
   			# printf("The S has only one point, so the convex hull is itself"); 
    			return S; 
  		end if; 
  		n:=nops(S[1]); 
  		for i from 2 to m do 
    			if n<>nops(S[i]) then 
      				error("Invalid input, the points must have same dimension"); 
    			end if; 
  		end do; 
  		r,indS:=affineDims(S);  
  		if r=1 then   #if these points are in the same line, return the end points 
    			end1:=indS[1]; end2:=indS[2];  
    			for i from 1 to m do 
                    for j from 1 to n do
                        if (S[i]-end2)[j]<>0 then
                            if ((S[i]-end2)[j]/(end1-end2)[j])>1 then
                                end1:=S[i];
                                break;
                            fi;
                            if  ((S[i]-end2)[j]/(end1-end2)[j])<0 then
                                end2:=S[i];
                                break;
                            fi;
                        fi;
                    end do;
    			end do; 
   			# printf("These points sre 1-affine independant, the convex hull is a line segment and the end points are as follows"); 
    			return end1,end2; 
  		elif r=n then      # if these points are n affine dimensional, return the facets  
    			F_V,H:=simplicialConvexHull(S,indS); 
    			certification(S,H); 
   			# printf("These points are %d -affine independant",n); 
    			return F_V,H; 
  		else   # if their affine dimension is less than n 
    			TrS,TrIndS,A,preBase:=affineTrans(S,indS,r);  
    			F_V,H:=simplicialConvexHull({op(TrS)},{op(TrIndS)}); 
    			certification(TrS,H); 
   			# printf("These points are %d -affine independant and the convex hull is computed after affine transformation to %d dimension vector space",r,r); 
    			return F_V,H; 
  		end if;   
	end proc: 
 
       



	#======================================================================================== 
	#Description 
	#  Given a set of points S and a list of points, this function returns the points in L which  
	#  is in the convex hull  of S. Assume the points are n dimensional, if their affine dimension 
	#  r isless than n, we affinely transformate these points into r dimensional vector space, 
	#  and judge if p belong the convex hull therein. 
	#Inputs:
    	#  S, L  --- sets of points
	#  LL    --- ponits in L which is in the convex hull of S		
	#======================================================================================== 
	inConvexHull:=proc(S::set,L) 
  		local i,j,F_V,H,n,r,m,indS,k,S_min,S_max,ind,TrP,x,TrS,TrIndS,A,preBase,LL,flag,l; 
  		m:=nops(S);LL:=[]; 
  		if m=0 then 
    			error("Invalid inputs, S must be a nonempty set"); 
  		end if; 
  		if m=1 then  
   			# printf("The S has only one point, so the convex hull is itself"); 
    			for i from 1 to nops(L) do 
      				if L[i]=S[1] then 
        				LL:=[op(LL),L[i]]; 
      				end if; 
    			end do;
    			return LL:
  		end if; 
  		n:=nops(S[1]); 
  		for i from 2 to m do 
    			if n<>nops(S[i]) then 
      				error("Invalid input, the points must have same dimension"); 
    			end if; 
  		end do; 
  		printf("computing the affine dimension \n"); 
  		r,indS:=affineDims(S);
  		m:=nops(L); 
  		printf("computing the convex hull \n"); 
  		if r=1 then           # if these points are in the same line 
    			S_min,S_max:=convexHull(S); 
                	for i from 1 to m do
                    		if L[i]=S_min or L[i]=S_max then
                        		LL:=[op(LL),L[i]];
                    		else
                        		for j from 1 to n do
                            			if (L[i]-S_max)[j]<>0 then
                                			if ((L[i]-S_max)[j]/(S_min-S_max)[j])>=0 and ((L[i]-S_max)[j]/(S_min-S_max)[j])<=1 then
                                    				LL:=[op(LL),L[i]];
                                			fi;
                                			break;
                            			fi;
                        		od;
                    		fi;
                	od;
                	return LL;
  		elif r=n then      # if these points are n-affine dimensional 
    			F_V,H:=simplicialConvexHull(S,indS); 
    			for i from 1 to m do  
      				flag:=0;  
      				for j from 1 to nops(H) do 
        				if evalf(Vector[row]([op(L[i]),-1]).H[j])>0 then   
          					flag:=1;
          					break; 
        				end if; 
      				end do;     
      				if flag=0 then     
        				LL:=[op(LL),L[i]];
      				end if;
    			end do;
    			return LL: 
  		else   #if these points are r-affine dimensional, do affine transformation 
    			#printf("affine dimension is less than n"); 
    			l:={};  
    			for i from 1 to m  do 
      				if (affineDims((indS union {L[i]})))[1]=r+1 then  
      					# if p is not in the affine hull of S, p is not in their convex hull, obviously. 
        				l:=l union {L[i]}; 
      				end if;     
    			end do;  
    			l:={op(L)} minus l;# return l; 
    			l:=[op(l)]; 
    			m:=nops(l);
    			TrS,TrIndS,A,preBase:=affineTrans(S,indS,r); 
    			F_V,H:=simplicialConvexHull({op(TrS)},{op(TrIndS)});  
   	 		for i from 1 to m do  
      				flag:=0; 
      				x:=LinearSolve(A,Vector(l[i]-preBase[1])); 
      				TrP:=IdentityMatrix(r).x;  
      				for j from 1 to nops(H) do 
        				if evalf(Vector[row]([op(convert(TrP,list)),-1]).H[j])>0 then 
          					flag:=1; 
          					break; 
        				end if;  
      				end do; 
      				if flag=0 then      
        				LL:=[op(LL),l[i]]; 
      				end if;
    			end do;
    			return LL;
  		end if;      
	end proc: 


	
end module;





CertifiTools:=module()
	description "package for computing and certifying the lower bound r of rational functions by SOS";
	local
        	nChoosek,
        	kth_subset,
		projSOS,
		roundApproxi,
		cutoff;
	export
		nVars_kth_monos,
		nVars2kth_monos,
		deg_mono,
		allDegree,
		getMonos,
		coeffs_monos,
		getSDP,
		certifiSOS,
		modifiSOS,
		isSemiPositive;
	option
		package;





	#=========================================================================== 
	#Description 
	#  Given n variables, return all monomials of total degree k 
	#Inputs 
	#  n     --- the number of variables 
	#  k     --- the total degree 
	#Output 
	#  monos --- a list of lists whose elements are the degree of each variable  
	#=========================================================================== 
	nVars_kth_monos:=proc(n::integer,k::integer) 
  		local i,j,N,l,m,monos; 
  		if k<0 then 
    			error("Invalid input, degree k must be nonnegative"); 
  		end if; 
  		if k=0 then
    			return [[seq(0,i=1..n)]];
  		end if;
  		if n=1 then
    			return [[k]];
  		end if;
  		N:=ConvexHull[kth_subset](1,n+k-1,n-1);#return N; 
  		m:=nops(N); 
  		monos:=[]; 
  		for i from 1 to m do 
    			l:=[N[i][1]]; 
    			for j from 2 to n-1 do 
      				l:=[op(l),N[i][j]-N[i][j-1]]; 
    			end do; 
    			if n>1 then
      				l:=[op(l),n+k-N[i][n-1]]; 
    			end if; #return l;
    			for j from 1 to nops(l) do 
      				l:=subsop(j=l[j]-1,l); 
    			end do; 
    			monos:=[op(monos),l]; 
  		end do; 
  		return monos; 
	end proc: 
 
 
	#=========================================================================== 
	#Description 
	#  Given n variables, return all monomials of total degree upto k 
	#Inputs 
	#  n     --- the number of variables 
	#  k     --- the total degree 
	#Output 
	#  L     --- a list of lists(degree of each variable)  
	#=========================================================================== 
	nVars2kth_monos:=proc(n::integer,k::integer) 
  		local i,L; 
  		if k<0 then 
    			error("Invalid input, degree k must be nonnegative"); 
  		end if; 
  		L:=[]; 
  		for i from 0 to k do 
    			L:=[op(L),op(nVars_kth_monos(n,i))]; 
  		end do; 
  		return L; 
	end proc:     
 
	#=========================================================================== 
	#Description 
	#  Given a monomial and a list of variables, return a list whose elements are  
	#  the degree of each variable 
	#Inputs 
	#  mono    --- a monomial 
	#  vars    --- a list of variables 
	#Output 
	#  L       --- the degree list 
	#=========================================================================== 
	deg_mono:=proc(mono,vars::list) 
  		local i,n,L,flag; 
  		if type(mono,polynom)=false or type(mono,`+`) then 
   			error("Invalid inputs, mono must be a monomial"); 
  		end if; 
  		n:=nops(vars); 
  		if n<=0 then 
   			error("Invalid inputs, the number of variables must be more than 0 "); 
  		end if; 
  		L:=[seq(0,i=1..n)]; 
  		for i from 1 to n do 
    			flag:=0; 
    			while flag=0 do 
      				if divide(convert(mono,rational),(vars[i])^(L[i])) then 
        				L:=subsop(i=L[i]+1,L); 
      				else 
        				flag:=1; 
        				L:=subsop(i=L[i]-1,L); 
      				end if; 
    			end do; 
  		end do; 
  		return L; 
	end proc:  
 
 
	#============================================================================ 
	#Description 
	#  Given a polynomial and a list of variables, return all the degree list of the  
	#  monomials apearing in f and the highest total degree 
	#Inputs 
	#  f     --- a polynomial 
	#  vars  --- a list of variables 
	#Output 
	#  L     --- the degree list 
	#  d     --- the total degree of f 
	#============================================================================ 
	allDegree:=proc(f,vars::list) 
  		local i,m,n,L,max,d,l; 
  		if type(f, polynom)=false then 
    			error("Invalid inputs, f must be a polynomial"); 
  		end if; 
  		if type(f,`+`)=false then  
    			L:=deg_mono(f,vars); 
    			d:=add(L[i],i=1..nops(L));
    			L:=[L];
    			return L,d; 
  		end if; 
  		n:=nops(vars); 
  		if n<=0 then  
    			error("Invalid inputs, variable number must be more than 0 "); 
  		end if; 
  		L:=[op(expand(f))]; 
  		m:=nops(L);  l:=[];
  		for i from 1 to m do 
    			l:=[op(l),deg_mono(L[i],vars)]; 
  		end do;   
  		max:=1; 
  		for i from 2 to m do 
    			if add(l[i][j],j=1..n)>add(l[max][j],j=1..n) then 
      				max:=i; 
    			end if; 
  		end do; 
  		d:=add(l[max][i],i=1..n); 
  		return l,d; 
	end proc: 
 
 
	#======================================================================== 
	#Description 
	#  Given polynomials f and g, return a list of monomials appearing  
	#  in the SOS relaxation: f-r*g=SOS
	#Inputs
	#  f     --- a polynomial
	#  vars  --- a list of variables
	#Output
	#  s     --- a list of monomials apearing in the SOS relaxation
	#======================================================================== 
	getMonos:=proc(f,g,vars::list) 
  		local i,j,m,n,L_f,L_g,d,L,s,mono; 
  		L_f,d:=allDegree(f,vars);  #return L_f,d; 
  		L_g,n:=allDegree(g,vars); #return L_f,L_g;
  		if d<n then 
			d:=n; 
		end if;
  		if type(d,even)=false then  
    			error("Invalid inputs, max(deg(f),deg(g)) must be even for r to be finite"); 
  		end if;
  		n:=nops(vars);
  		if n<=0 then
    			error("Invalid inputs, the number of variables must be more than 0"); 
  		end if; 
  		L:=nVars2kth_monos(n,d/2); 
  		m:=nops(L); 
  		s:=ConvexHull[inConvexHull]({op(L_f),op(L_g)},2*L); 
 		s:=1/2*s;
  		for i from 1 to nops(s) do 
    			mono:=1; 
    			for j from 1 to n do 
      				mono:=mono*(vars[j])^(s[i][j]); 
    			end do; 
      			s:=subsop(i=mono,s); 
  		end do; 
 	 	return s; 
	end proc: 
 
 
	#=================================================================
	# Discription:
	#    Given a polynomial f and a list of monomials, coeffs_monos returns
	#    a list of coefficients of the monomials in f
	# Inputs:
	#    f      --- a given polynomial
	#    L      --- a given list of monomials
	# Outputs:
	#    l      --- an list of coefficients 
	#================================================================= 
	coeffs_monos:=proc(f,L)
  		local m,n,i,j,l,L_f,vars;
  		if type(f, polynom)=false then
    			error("Invalid inputs, f must be a polynomial");
  		end if;
  		m:=nops(L);
  		for i from 1 to m do
    			if type(L[i],polynom)=false then
      				error("Invalid inputs, L must be a list of monomials");
    			end if;
  		end do;
  		l:=[seq(0,i=1..m)];
  		if type(f,`*`) then
    			for i from 1 to m do
      				if f/coeffs(f,[op(indets(L[i]))])=L[i] then
        				l:=subsop(i=f/L[i],l);
      				end if;
    			end do;
    			return l;
  		end if; 
  		vars:={};
  		for i from 1 to m do
    			vars:= vars union indets(L[i]);
  		end do;
  		vars:=[op(vars)];
  		L_f:=[op(expand(f))]; #return L_f,L,vars;
  		n:=nops(L_f);
  		for i from 1 to m do
    			for j from 1 to n do
      				if L_f[j]/coeffs(L_f[j],vars)=L[i] then
        				l:=subsop(i=l[i]+coeffs(L_f[j],vars),l);
      				end if;
    			end do;
  		end do;
  		return l;
	end proc:



	#=======================================================================
	# Discription:
	#    Return the data of the SDP problem for solving a SOS relaxation: f-r*g=SOS 
	# Inputs:
	#    f, g --- two given polynomials
	#    vars --- the variables which appear in f and g
	#    SDPmonomials -- the monomials list in f-r*g=m_d^T.W.m_d
	#    SDPM1 --- upper bound of the trace of the Gram matrix
	#    SDPM2 --- uppper bound of the moment matrix
	# Outputs:
	#    L,c,x0,Z0 --- the data for solving the SDP
	#    monos     --- monomials vector 
	#========================================================================
	getSDP:=proc(f::polynom, g::polynom, vars::list, {SDPmonomials:=[]}, {SDPM1::posint:=0}, {SDPM2::posint:=0},$)  
   		local i,j,k,l,m,n,L,c,s,x0,Z0,monos,M1,M2;       
   		s:={};  
   		if SDPmonomials=[] then    
     			monos:=getMonos(f,g,vars); print(monos);
   		else 
     			monos:=SDPmonomials;
   		end if; 
        monos:=convert(monos,Vector);
   		n:=Dimension(monos);
   		for i from 1 to n do
     			for j from 1 to n do 
       				s:=s union {monos[i]*monos[j]};
     			end do;  
   		end do;
   		s:=[op(s)]; 
   		l:=coeffs_monos(g,s);
   		L:=[seq(Matrix(n+2),i=1..nops(s))];
   		for i from 1 to n do 
     			for j from 1 to n do
       				for k from 1 to nops(s) do
         				if s[k]=monos[i]*monos[j] then  
           					L[k][i,j]:=1;
         				end if;
         				L[k][n+1,n+1]:=l[k];
	     				L[k][n+2,n+2]:=-l[k];
       				end do;   
     			end do;
   		end do; 
   		L:=[Matrix(n+2),op(L)];
   		L[1][n+1,n+1]:=-1; 
   		L[1][n+2,n+2]:=1;
   		c:=coeffs_monos(f,s); 
		if SDPM1=0 then
			M1:=100*max(op(map(i->abs(i),[coeffs(f),coeffs(g)])));
		else 
			M1:=SDPM1;
		end if;
		if SDPM2=0 then
			M2:=100*max(op(map(i->abs(i),[coeffs(f),coeffs(g)])));
		else 
			M2:=SDPM2;
		end if; 
   		L,c,x0,Z0:=SDPTools[BigM](L,c,M1,M2);   
     	return L,c,x0,Z0,monos;  
	end proc:




	#====================================================================================================================
	#Description
	#  Given a rational function f/g, return the lower bound r, if f-r*g=SOS exsits. When you run CertifiSOS at first time, 
    #  it is no need to input the SDP data and it can get r if the iterations are enough. If it does return
	#  a certified r, or you want to get a better one, run the second calling squence with the result of the first running.
	#  Then it will continue to compute r.
	#Inputs:
	#  f, g          ---  polynomials
	#  iter_num      ---  number of iterations
	#  F0_Fm, c_T    ---  data of the SDP from the SOS relaxation
	#  x_SF, Z_SF    ---  the strictly feasible points 
	#  SOSmonomials  ---   monomials vector apearing in the SOS relaxation, get from the first running
	#  SOSM1 --- upper bound of the trace of the Gram matrix
	#  SOSM2 --- upper bound of the trace of the moment matrix
	#Outputs
	#  L, c, X, Z    ---  described above, to be uesed the next running
	#  minval        ---  the cerfitified lower bound
	#  Q             ---  Gram matrix
	#  monos         ---  described above
	#====================================================================================================================
	certifiSOS:=proc(f::polynom, g::polynom, iter_num::posint, {SOSmonomials::Vector:=[]}, {SOSM1::posint:=0}, {SOSM2::posint:=0}, {c_T::list:=[]}, {F0_Fm::list:=[]}, {x_SF::list:=[]}, {Z_SF::Matrix:=0},$)
   		local i,j,k,monos,m,M1,M2,X,Z,c,L,minval,resid,Q,coeffs_,polys,tol,t,LL,d; 
		if SOSM1=0 then
                         M1:=100*max(op(map(i->abs(i),[coeffs(f),coeffs(g)])));
		else 
			M1:=SOSM1;
                end if;
                if SOSM2=0 then
                        M2:=100*max(op(map(i->abs(i),[coeffs(f),coeffs(g)])));
		else 
			M2:=SOSM2;
                end if; 
     		if SOSmonomials=[] then	
     			L,c,X,Z,monos:=getSDP(f,g,[op({op(indets(f)),op(indets(g))})],SDPM1=M1,SDPM2=M2);
     		else 
			monos:=SOSmonomials;
			if c_T=[] or F0_Fm=[] or x_SF=[] or Z_SF=0 then
				L,c,X,Z,monos:=getSDP(f,g,[op({op(indets(f)),op(indets(g))})],SDPmonomials=monos,SDPM1=M1,SDPM2=M2);
			else
				L:=F0_Fm; c:=c_T; X:=x_SF; Z:=Z_SF;
			end if;
     		end if;	 # return c,X,Z,L;
   		Q:=Matrix(Dimension(monos),Dimension(monos));
   		for k from 1 to iter_num do
     			t:=time();
			printf("************************************************************************************ \n"); 
			printf("************************************************************************************ \n");
     			printf("the %d iteration begins \n",k); 
     			X,Z:=SolveSDP(c,X,Z,L,1,v=5);
     			print(`the SDP objective is  `,-Trace(L[1].Z));
     			for i from 1 to Dimension(monos) do  
       				for j from 1 to i do 
           				Q[i,j]:=Z[i,j];
       				end do;
     			end do;
     			for i from 1 to Dimension(monos) do
       				for j from i+1 to Dimension(monos) do
           				Q[i,j]:=Q[j,i];
       				end do;
     			end do;
     			minval:=Z[Dimension(monos)+1,Dimension(monos)+1]-Z[Dimension(monos)+2,Dimension(monos)+2];
     			resid:=norm(expand(f-minval*g-Transpose(monos).Q.monos),2); 
     			print(`At this iteration the lower bound  is`); 
     			print(minval);
     			printf("At this iteration the residual error before projection is\n ");
     			print(resid);  
     			minval,Q:=projSOS(f,g,minval,monos,Q);
     			LL,d:=isSemiPositive(Q);
     			if LL=[] then
       				print(`This lower bound is not certified and need more iterations \n`);
       				print(`The time for this iteration is `,time()-t);
       				next;
     			else
       				print(`The Gram matrix Q is semidefinite positive and this lower bound is certified. \n`);
       				LL:=Transpose(monos).LL; 
       				tol:=evalf(expand(f-minval*g-add(d[i]*LL[i]^2,i=1..Dimension(d))));
       				print(`tolerance is`,tol);
       				print(`The time for this iteration is`,time()-t);
     			end if;
     			printf("the %d iteration ends \n",k);
	 		printf("************************************************************************************ \n");
			printf("************************************************************************************ \n");
   		end do; 
   		return minval,L,c,X,Z,Q,monos; 
	end proc:







	#=========================================================================
	#Description
	#  When we use certiSOS to get a high accuracy Gram matrix, if it is not 
	#  positive semidefinite, modifiSOS can round and cutoff the digits in Gram 
	#  matrix and minval to try to get a positve semidefinite matrix
	#Inputs:
	#  f,g,Z,monos    --- results from the certifiSOS
	#  digit1         --- if the difference between Z[i,j] and the nearest integer 
	#                     number is less than 10^(digit1), then return the nearest integer
	#  digit2         --- retain Z[i,j] digit2 figures after the decimal point
	#Outputs
	#  if it gets a positive semidefinite Gram matrix, then return:
	#  minval         --- certified lower boud
	#  Q              --- positive semidefinite matrix
	#  LL             --- list of polynomials, satisifies  f-minval*g-add(d[i]*LL[i]^2,i=1..Dimension(d))=0
	#  d              --- coefficients list 
	#=========================================================================
	modifiSOS:=proc(f,g,Z,monos::Vector,digit1,digit2)
  		local i,j,tol,resid,minval,Z1,m,n,Q,LL,d; 
  		n:=Dimension(monos); 
  		m:=RowDimension(Z); 
  		Z1:=Matrix(m); 
  		for i from 1 to m do  
    			for j from 1 to m do 
      				Z1[i,j]:=cutoff(roundApproxi(Z[i,j],digit1),digit2); 
    			end do; 
  		end do; 
  		minval:=Z1[n+1,n+1]-Z1[n+2,n+2];
  		Q:=Z1[1..n,1..n];
  		resid:=norm(expand(f-minval*g-Transpose(monos).Q.monos),2); 
  		print(`The residual error after modification and before projection is`); 
  		print(resid);
  		minval,Q:=projSOS(f,g,minval,monos,Q); 
  		LL,d:=isSemiPositive(Q);
  		if LL=[] then
    			print(`This lower bound is not certified and need more iterations \n`);
  		else
    			print(`The Gram matrix Q is semidefinite positive and this lower bound is certified. \n`);
    			LL:=Transpose(monos).LL; 
    			tol:=evalf(expand(f-minval*g-add(d[i]*LL[i]^2,i=1..Dimension(d))));
    			print(`tolerance is`,tol);
  		end if;
  		return minval,Q,LL,d; 
	end proc:






	#==================================================================
	#Description
	#   
	#==================================================================
	projSOS:=proc(f,g,lower,monos::Vector,Q)
  		local i,j,k,d,m,n,rf,rg,rQ,rlower,c1,c2,mv,eqs,W;
  		rf:=convert(f,rational,exact);
  		rg:=convert(g,rational,exact);
  		rlower:=convert(lower,rational,exact);
  		rQ:=map(i->convert(i,rational,exact),Q);
  		d:=expand(rf-rlower*rg-Vector[row](monos).rQ.Vector(monos));
  		if d=0 then 
    			return rlower,(Transpose(rQ)+rQ)/2; 
  		else
    			if type(d,`+`) then
      				d:=[op(d)];
    			else
      				d:=[d];
    			end if;
    			c1:=[seq(coeffs(d[i]),i=1..nops(d))];
    			c2:=[seq(d[i]/coeffs(d[i]),i=1..nops(d))];
    			eqs:=[seq([],i=1..nops(c2))]; 
    			mv:=[];
    			n:=Dimension(monos);
    			for i from 1 to n do
      				for j from 1 to n do 
        				for k from 1 to nops(c2) do
          					if monos[i]*monos[j]=c2[k] then 
               						eqs:=subsop(k=[op(eqs[k]),[i,j]],eqs);
          					end if;
        				end do;
      				end do;
    			end do; 
    			W:=Matrix(n);
    			for i from 1 to nops(eqs) do
      				for j from 1 to nops(eqs[i]) do
        				W[eqs[i][j][1],eqs[i][j][2]]:=c1[i]/nops(eqs[i]);
      				end do;
    			end do; 
    			W:=W+rQ; 
    			W:=(Transpose(W)+W)/2;
    			return rlower,W;
  		end if;
	end proc:


	#==================================================================
	#Description
	#  Given a matrix A, this function judges that if A is semidefinite
	#  positive matrix. If it is, return the LDL decomposition.
	#==================================================================
	isSemiPositive:=proc(A::Matrix)
  		local i,j,k,m,n,L,d,M;
  		m,n:=Dimension(A);
  		if m<>n then
    			error("the input matrix must be a square matrix \n");
  		end if;
  		M:=Transpose(A); 
  		M:=M-A;
  		for i from 1 to n do
    			for j from 1 to n do
      				if M[i,j]<>0 then
        				printf("The input matrix is not a symmetrical \n");
        				return false;
      				end if;
    			end do;
  		end do; 
  		L:=Matrix(n); d:=Vector(n);
  		for k from 1 to n do   # assume we have known the k-1 column of L,then 
    			for i from k to n do 
      				L[i,k]:=A[i,k];
      				for j from 1 to k-1 do
        				L[i,k]:=L[i,k]-d[j]*L[i,j]*L[k,j];
      				end do;
    			end do;
    			if L[k,k]>0 then 
      				d[k]:=L[k,k];
      				L[k..n,k]:=L[k..n,k]/d[k];
    			elif Norm(L[k..n,k])>0 then
      				printf("It is not a semidefinite matrix \n");
      				return [],[];   
    			end if;  
  		end do;
  		return L,d;
	end proc:




	#====================================================================================
	#Description
	#  Given a float number f, if the difference between f and the nearest integer number
	#  id less than 10^(digit), then return the nearest integer, else return f.
	#====================================================================================
	roundApproxi:=proc(f,digit) 
   		local num,den,temp; 
   		if type(f,float)=false and type(f,integer)=false then
     			error("Invalid, the 1th argument requires to be a float or integer"); 
   		end if;
   		if whattype(digit)<>integer or digit<0 then
     			error("Invalid input, the 2th argument required to be a nonnegative integer");
   		end if;
   		temp:=abs(f);
   		num:=numer(convert(temp,rational));
   		den:=denom(convert(temp,rational)); #return num,den;
   		if evalf(irem(num,den)/den)<10^(-digit) then
     			if f>=0 then
       				return iquo(num,den);    
     			else
       				return -iquo(num,den); 
     			end if; 
   		elif evalf(1-irem(num,den)/den)<10^(-digit) then
     			if f>=0 then
       				return iquo(num,den)+1;
     			else 
       				return -(iquo(num,den)+1);  
     			end if;
   		else 
     			return f;
   		end if;
	end proc:






	#=======================================================================
	#Description
	#  Given a float point, retain its digit figures after the decimal point
	#=======================================================================
	cutoff:=proc(f,digit)
   		local num,den,temp; 
   		if type(f,float)=false and type(f,integer)=false then
     			error("Invalid, the 1th argument requires to be a float or integer");
   		end if;
   		if whattype(digit)<>integer or digit<0 then   
     			error("Invalid input, the 2th argument required to be a nonnegative integer"); 
   		end if; 
   		if type(f,integer) then
     			return f;
   		end if;
   		temp:=abs(f);
   		temp:=convert(temp*10^(digit),rational); 
   		num:=numer(temp); den:=denom(temp); 
		temp:=iquo(num,den)*10^(-digit); 
   		if f>=0 then
     			return temp; 
   		else 
     			return -temp; 
   		end if; 
	end proc:




end module;