(*
 *		Global controls.
 *)

showproof		=	False;
showreason		=	True;
showHeuristics	=	True;

(*
 *		We represent the simple symmetric sum, [a,b,c],
 *		by the function SS[a,b,c].
 *		This function can have any number of arguments
 *		(but at least 2).
 *		We always print this function as "[a,b,c]".
 *)

Format[SS[x__]]:=SequenceForm["[",Infix[{x},","],"]"];

(*
 *		Fix the builtin "Variables" function to work properly
 *		with inequalities.
 *)

Unprotect[Variables];
Variables[x_> y_]:=Union[Variables[x],Variables[y]];
Variables[x_>=y_]:=Union[Variables[x],Variables[y]];
Variables[x_< y_]:=Union[Variables[x],Variables[y]];
Variables[x_<=y_]:=Union[Variables[x],Variables[y]];
Variables[x_==y_]:=Union[Variables[x],Variables[y]];
Protect[Variables];

(*
 *	Create a function that finds the apparent sign of a term.
 *)

apparentSign[x_?NumberQ] := If[x<0,-1,1,1];
apparentSign[x_?AtomQ] := 1;
apparentSign[x_+y_] :=
	If[apparentSign[x]===apparentSign[y],
		apparentSign[x],
		unknownSign
	  ];
apparentSign[x_ * y_] :=
	If[IntegerQ[apparentSign[x]]
	 &&IntegerQ[apparentSign[y]],
		apparentSign[x]*apparentSign[y],
		unknownSign
	  ];
apparentSign[x_^n_Integer] :=
	If[IntegerQ[apparentSign[x]],
		apparentSign[x]^n,
		unknownSign
	  ];
apparentSign[Sqrt[x_]] := 1;

(*
 *	Move a term to the lhs or rhs depending on its sign.
 *)

Sep[term_]:=
Block[{numf},
	numf = numfactor[term];
	If[!NumberQ[numf],
		numf = Factor[numf];
		numf = First[numf];
	  ];
	
	If[!NumberQ[numf],
		If[apparentSign[term]==1, numf=1];
		If[apparentSign[term]==-1, numf=-1];
	  ];

	If[numf<0,
		right -= term,
		left += term,
		left += term
		];
	];

checksymm[expr_]:=
Block[ {vars},
 vars = Variables[expr];
 Return[symm[expr,vars,checkit,expr]];
	];

checkit[expr1_,expr2_] :=
Block[{temp},
	If[Together[expr1-expr2]==0,Return[True]];
	Return[False];
	];

(*
 *	Move all positive terms to the left
 *	and negative terms to the right.
 *)

signseparate[]:=
(* Implicit arguments: left, right, op. *)
Block[{expr},
 expr = Expand[Together[left-right]];

 left = expr;
 right = 0;

 If[NumberQ[expr] && expr<0,
	right =	-left;
	left =	0;
	Return[];
   ];

 If[AtomQ[expr],Return[]];

 If[Head[expr] =!= Plus,
  		If[apparentSign[expr]==(-1),
			right =	-left;
			left =	0;
  		  ];
  		Return[];
    ];

 (*
 	The expression is a sum of terms.
    Scan these terms and move the
    positive ones to the left
    and the negative ones to the right.
 *)

 left =  0;
 right = 0;
 Map[Sep,expr];
 Return[];
 	];

(*
 *	Collect together terms of the same sign.
 *)

collect[]:=
 (* Implicit arguments: left, right, op. *)
 Block[ {oldleft,oldright},

 oldleft =	left;
 oldright =	right;

 signseparate[];

 If[showproof &&
 	(Together[left-oldleft]=!=0 ||
 	 Together[right-oldright]=!=0),
    	Print["Expanding and collecting terms of the same sign gives:"];
		Print[op[left,right]];
    ];

];

(*
 *	Break an inequality into 3 parts:
 *	the lhs, the operator, and the rhs.
 *	Store in "left", "op", and "right" respectively.
 *)

getop[inequal_]:=
Block[{op},
 op = Head[inequal];

 left =  Part[inequal,1];
 right = Part[inequal,2];

 If[op =!= Less  &&
    op =!= LessEqual &&
    op =!= Greater  &&
    op =!= GreaterEqual &&
    op =!= Equal,
		Print["Expression, ",inequal," is not an inequality."];
        Print["Main operator is: ", op];
        Error[inequal];
   ];

 Return[op];
	];

homogineq[inequal_,varlist_]:=
Block[{left,right,op,temp,vleft,vright,vl,vr,result},

 If[showproof,showreason=True];

 op = getop[inequal];

 normalorder[];

 clearfractions[];

 result = checkobvious[];
 If[result==True || result==False,Return[result]];

 gcdize[];

 collect[];

 vleft =	vec[left,varlist];
 vright =	vec[right,varlist];

 (* clear fractions by multiplying both
    sides by the lcm of the denominators *)

 vl =	Apply[Plus,vleft];
 vl =	Denominator[Together[vl]];
 vr =	Apply[Plus,vright];
 vr =	Denominator[Together[vr]];
 temp =	PolynomialLCM[vl,vr];
 vleft =	temp*vleft;
 vright =	temp*vright;

 If[showreason,
    Print["Expressing in terms of symmetric polynomials gives:"];
	Print[op[Apply[Plus,vleft],Apply[Plus,vright]]];
   ];

 result = majorize[vleft,vright];

 If[supressbanner=!=True && result=!=True,
    Print["left  side: ",vleft];
	Print["right side: ",vright];
   ];

 Return[result];

	];

ApplyMuirhead[reason_,Lindex_,Rindex_]:=
(* implicit inputs and outputs:
	leftsum, rightsum, left, right, change
*)
Block[{lentry,rentry,con,
	rawleft,rawright,lfactor,rfactor},
	lentry =	left[[Lindex]];
	rentry =	right[[Rindex]];
	lfactor =	numfactor[lentry];
	rawleft =	lentry/lfactor;
    rfactor =	numfactor[rentry];
    rawright =	rentry/rfactor;
	con =		Min[lfactor,rfactor];
	If[showreason,
	    Print[con*rawleft>=con*rawright];
		If[showHeuristics,
			Print["... by Heuristic ",reason,"."]];
	  ];
	rightsum =	rightsum - con*rawright;
	leftsum =	leftsum - con*rawleft;
	left =		listerize[leftsum];
	right =		listerize[rightsum];
	change =	True;
	];

(*
	majorize takes two vectors of symmetric sums [a,b,c],
	and returns true if the left one majorizes the right one
	or if the inequality can be proven false.
	It returns false if right majorizes left and returns
	notsure if it is not sure.
*)

majorize[lefty_,righty_] :=
Block[{lfactor,rfactor,rawright,rawleft,len,result,termsize,p,d,m1,m2,con,
	 leftsum,rightsum,whileresult,change,temp,lcount,lentry,
	 rcount,rentry,res2,every,some,Li,Ri,applied},

 left = lefty;
 right = righty;

 	(*
 	 Heuristic 0:  If there is one term on each side,
		  then the result is determined by Muirhead's Theorem.
	 *)

 If[Length[left]==1 && Length[right]==1,
 	left =		First[left];
	right =		First[right];

	result =	major[left,right];

	If[showreason && result==True,
	    Print["This result follows from the majorization theorem."];
	  ];

	Return[result];
	];

 If[showreason,
    Print["This follows from the following majorizations:"];
	];

While[Length[right]>0,
 leftsum =	Apply[Plus,left];
 rightsum =	Apply[Plus,right];

 change = False;
 
 (* Heuristic 1:
    See if there is some term on the left that does not majorize
    anything on the right.  If that is the case, then see if we
    can apply Schur's theorem, either

			[p+2d,0,0]  + [p,d,d]  > 2 [p+d,d,0]		case a

    or		[p+d,p+d,0] + [2p,d,d] > 2 [p+d,d,p]		case b, d<=2p

    or		[p+d,p+d,0] + [2p,d,d] > 2 [p+d,d,p]		case c, d>2p
 *)

 temp =		First[left];
 temp =		temp/numfactor[temp];
 termsize =	Length[temp];

 If[termsize==3,
  For[i=1,i<=Length[left],i++,
	lfactor =	numfactor[left[[i]]];
    rawleft =	left[[i]]/lfactor;

    For[j=1, j<=Length[right],j++,
		rfactor =	numfactor[right[[j]]];
		rawright =	right[[j]]/rfactor;
		If[major[rawleft,rawright],Break[]];
       ];
    If[j<=Length[right],Continue[]];

    If[Part[rawleft,2]==Part[rawleft,3], (* i.e. [x,d,d] *)
		d =	Part[rawleft,2];
	    p =	Part[rawleft,1];
	    m1 =	Coefficient[leftsum,SS[p+2*d,0,0]];
	    m2 =	Coefficient[rightsum,SS[p+d,d,0]];
	    If[m1=!=0 && m2>1,
	    	con = Min[lfactor,m1,Floor[m2/2]];
			leftsum =  leftsum-con*SS[p+2*d,0,0]-con*SS[p,d,d];
			rightsum = rightsum-2*con*SS[p+d,d,0];
			If[showreason,
		    	Print[con*SS[p+2*d,0,0]+con*SS[p,d,d]>=2*con*SS[p+d,d,0]];
				If[showHeuristics,Print["... by Heuristic 1a."]];
			  ];
			left =	listerize[leftsum];
			right =	listerize[rightsum];
			change = True;
			Break[];
	    (* else *),
	    If[EvenQ[p],
	   		p =	p/2;
			m1 =	Coefficient[leftsum,SS[p+d,p+d,0]];
			m2 =	Coefficient[rightsum,SS[p+d,d,p]];
			If[m1!=0 && m2>1,
				con = Min[lfactor,m1,Floor[m2/2]];
				leftsum = leftsum-con*SS[p+d,p+d,0]-con*SS[2p,d,d];
				rightsum = rightsum-2*con*SS[p+d,d,p];
				If[showreason,
			 	 	Print[con*SS[p+d,p+d,0]+con*SS[2p,d,d]>=2*con*SS[p+d,d,p]];
					If[showHeuristics,Print["... by Heuristic 1b."]];
			 	  ];
				left =	listerize[leftsum];
				right =	listerize[rightsum];
				change = True;
				Break[];
			  ];
	      ]; (* end if EvenQ *)
	   ]; (* end if m1!=0 *)
     ];	(* end if [x,d,d] *)
    If[Part[rawleft,1]==Part[rawleft,2], (* i.e. [d,d,x] *)
		d =	Part[rawleft,1];
	    p =	Part[rawleft,3];
	    If[EvenQ[p],
	    	p = p/2;
		    m1 =	Coefficient[leftsum,SS[p+d,p+d,0]];
		    m2 =	Coefficient[rightsum,SS[p+d,d,p]];
		    If[m1=!=0 && m2>1,
		    	con = Min[lfactor,m1,Floor[m2/2]];
				leftsum =  leftsum-con*SS[p+d,p+d,0]-con*SS[d,d,2p];
				rightsum = rightsum-2*con*SS[p+d,d,p];
				If[showreason,
			    	Print[con*SS[p+d,p+d,0]+con*SS[d,d,2p]>=2*con*SS[p+d,d,p]];
					If[showHeuristics,Print["... by Heuristic 1c."]];
				  ];
				left =	listerize[leftsum];
				right =	listerize[rightsum];
				change = True;
				Break[];
		     ]; (* end if m1!=0 *)
		   ]; (* end if EvenQ *)
     ];	(* end if [d,d,x] *)
    ];	(* end for i *)
  If[i<=Length[left],Continue[]];

  ]; 	(* end if termsize==3 *)

(*
	 Heuristic 2: If there is only one term on the right,
		 then the terms on the left must majorize it away.
 *)

 If[Length[right]==1,
	len =		Length[left];

	For[i=len, i>0, i--,
		If[Length[right]==0,Return[True]];
		rfactor =	numfactor[right[[1]]];
		If[rfactor<=0,Return[True]];
		rawright =	right[[1]]/rfactor;
	    lfactor =	numfactor[left[[i]]];
	    rawleft =	left[[i]]/lfactor;
	    If[major[rawleft,rawright],
	    	ApplyMuirhead["2",i,1];
	      ];
	   ];
	];

 If[leftsum==0 && rightsum=!=0,Return[False]];
 If[rightsum==0,Return[True]];

 (*
 	Heuristic 3:
    If there is a term on the right that is majorized by
    precisely one term on the left, then we may as well
    remove that term now.
 *)

  applied = False;
  For[j=1, j<=Length[right], j++,
    rfactor =	numfactor[right[[j]]];
    rawright =	right[[j]]/rfactor;
    lcount =	0;

    For[i=1, i<=Length[left], i++,
		lfactor =	numfactor[left[[i]]];
		rawleft =	left[[i]]/lfactor;
		If[major[rawleft,rawright],
	    	lcount++;
			Li = i;
	      ];
       ];

    If[lcount==1,
    	ApplyMuirhead["3",Li,j];
    	applied = True;
		Break[];
      ];
   ];
 If[applied,Continue[]];

 If[leftsum==0 && rightsum=!=0,Return[False]];
 If[rightsum==0,Return[True]];

 (* Heuristic 4:
    If there is a term on the left that majorizes
    precisely one term on the right, then we may as well
    remove that term now.
 *)

  applied = False;
  For[j=1, j<=Length[left], j++,
    lfactor =	numfactor[left[[j]]];
    rawleft =	left[[j]]/lfactor;
    rcount =	0;

    For[i=1, i<=Length[right], i++,
		rfactor =	numfactor[right[[i]]];
		rawright =	right[[i]]/rfactor;
		If[major[rawleft,rawright],
	    	rcount++;
	    	Ri = i;
	      ];
       ];

    If[rcount==1,
    	ApplyMuirhead["4",j,Ri];
    	applied = True;
		Break[];
      ];
   ];
 If[applied,Continue[]];

 If[leftsum==0 && rightsum=!=0,Return[False]];
 If[rightsum==0,Return[True]];

 (* Heuristic 5:
    If every term on the left majorizes every term on the right,
    then we can pair them off any way we want.
 *)

  every = True;
  For[j=1, j<=Length[right], j++,
    rfactor =	numfactor[right[[j]]];
    rawright =	right[[j]]/rfactor;

    For[i=1, i<=Length[left], i++,
		lfactor =	numfactor[left[[i]]];
		rawleft =	left[[i]]/lfactor;
		If[!major[rawleft,rawright],
			every = False;
			Break[];
		  ];
       ];
     If[!every,Break[]];
     ];

  If[every,
  	ApplyMuirhead["5",1,Length[right]];
    ];

 If[leftsum==0 && rightsum=!=0,Return[False]];
 If[rightsum==0,Return[True]];

 (* Heuristic 6:
    If some term on the left majorizes some term on the right,
    then we can remove this pair.
 *)

  For[j=Length[right], j>0, j--,
    rfactor =	numfactor[right[[j]]];
    rawright =	right[[j]]/rfactor;

	some = False;
    For[i=1, i<=Length[left], i++,
		lfactor =	numfactor[left[[i]]];
		rawleft =	left[[i]]/lfactor;
		If[major[rawleft,rawright],
			ApplyMuirhead["6",i,j];
			some = True;
			Break[];
		  ];
   	   ];
   If[some, Break[]];
   ];

If[leftsum==0 && rightsum!=0,Return[False]];
If[rightsum==0,Return[True]];

If[!change, Break[]];

];

 If[Length[right]==0,Return[True]];
 If[Length[left]>1 && Length[right]>1,Return[notsure]];
 Return[notyetimplemented];

];
 
(*
 *	Turn n1[a1,b1,c1]+n2[a2,b2,c3]+n3[a3,b3,c3] into the list
 *	{n1[a1,b1,c1] , n2[a2,b2,c3] , n3[a3,b3,c3]}.
 *)

listerize[expr_]:=
Block[{temp},
 	If[AtomQ[expr],
		If[expr===0, Return[{}], Return[{expr}]];
	  ];
 If[Head[expr]===Plus,
 	Return[ReplacePart[expr,List,0]];
   ];
   	
 Return[{expr}];
 ];

(* major takes two symmetric sums, such as [a,b,c] and [d,e,f]
   and returns true if the left one majorizes the right one.
*)

major[left_,right_] :=
Block[{lfactor,rfactor,lterm,rterm,lsum,rsum,len,result,ineqList},
	(* Extract a and {x,y,z} from a{x,y,z} *)
	{lfactor,lterm} =	FactorTermsList[left];
	{rfactor,rterm} =	FactorTermsList[right];
	If[lfactor < rfactor,Return[funnyratio]];

	(* turn SS[a,b,c] into normal list {a,b,c} *)

 	lterm = lterm/.{SS[x__]->{x}};
 	rterm = rterm/.{SS[x__]->{x}};

	If[Length[lterm] != Length[rterm],
		Print["Length error."];
	 	Print[lterm," ",rterm," ",lfactor," ",rfactor];
		Return[lengthsdiffer];
	  ];
	len = Length[lterm];

	(* compute the partial sums of each one
	and check for majorization *)

	lsum = Drop[FoldList[Plus,0,lterm],1];
	rsum = Drop[FoldList[Plus,0,rterm],1];

	(* form sequence of inequalities *)
	ineqList = MapThread[GreaterEqual,{lsum,rsum}];
	
	(* "and" them all together *)
	
	result = Apply[And,ineqList];
	Return[result];
	];

(*
	Call a function "fn" once for each
	permutation of the variables in "expr".
	The list of variables is given in "varlist".
	The function is called as
		fn(permutedExpr,arg)
	and returns True to continue, False to abort.
*)
symm[expr_,varlist_,fn_,arg_]:=
Block[{val,permList,len,i,perm,permExpr,permRule},
	permList = Permutations[varlist];
	len = Length[permList];
	For[i=1,i<=len,i++,
		perm = permList[[i]];
		permRule = MapThread[Rule,{varlist,perm}];
		permExpr = expr/.permRule;
		val = fn[permExpr,arg];
		If[val==False,Return[False]];
	   ];
	Return[True];
	];

numfactor[term_]:=
Block[{flist},
	flist = FactorTermsList[term];
	Return[First[flist]];
	];


vec[exp_,varlist_] :=
Block[ {resultvec,vector,term,size,sym,nn,
 	co,coef,ratio,newsize,
	powlist,powterm,powco,expr},

	expr = exp;
	resultvec = {};
	expr =		Expand[Together[expr]];
	nn = 		Length[varlist];

	While[expr=!=0,
    	If[isterm[expr],
	    	size = 1,
	    	size = Length[expr]
	      ];
		If[isterm[expr],
	    	term = expr,
	    	term = First[expr]
	      ];
		coef =	numfactor[term];
		term =	term/coef;
		sym =	symmetrize[term,varlist];
		co =	Coefficient[sym,term];
		If[isterm[sym],
		  	co = nn!,
		  	co = numfactor[First[sym]]
		  ];
		If[co==0,
	  		Print["co=0 term=",term," sym=",sym,
	  			" expr=",expr," coef=",coef];
	  	  ];
		ratio =	coef/co;

		vector =	Exponent[term,varlist,Identity];

		powlist =	MapThread[Power,{varlist,vector}];
		powterm =	Apply[Times,powlist];
		powco =		Together[term/powterm];

		vector =	Apply[SS,Reverse[Sort[vector]]];
		vector =	ratio*powco*vector;
		PrependTo[resultvec,vector];
		expr =		Expand[expr - ratio*sym];
		If[isterm[expr],
	    	newsize = 1,
	   		newsize = Length[expr]
	   	  ];
		If[newsize >= size && expr =!= 0,
			Print["internal error, expression got larger"];
			Print["term=",term," ratio=",ratio," coef=",coef," co=",co];
			Print["size=",size," newsize=",newsize," sym=",sym];
			Print[expr];
			error[expr];
		  ];
    	];

    Return[resultvec];
	];

isterm[expr_] :=
Block[{},
	If[AtomQ[expr],Return[True]];
	If[Head[expr]=!=Plus,Return[True]];
	Return[False];
	];

symmetrize[term_,varlist_]:=
Block[ {result},

	cumulate[tem_,arg_]:= {result+=tem;Return[True];};

	result = 0;
	symm[term,varlist,cumulate,0];
	Return[result];
	];

test[expr_] :=
Block[ {varlist,result},
	varlist = Variables[expr];
	Print[" "];
	Print["----------"];
	Print[expr];
	result = homogineq[expr,varlist];
	If[result=!=True,
		Print["********* COULD NOT PROVE THIS INEQUALITY *********"];
	  ];
	];

clearfractions[] :=
(* Implicit arguments: left, right, op. *)
Block[ {dl,dr,prod,pr},
	left =	Together[left];
	right =	Together[right];
	dl = Denominator[left];
	dr = Denominator[right];
	prod = PolynomialLCM[dl,dr];

 	If[prod=!=1,
   		left =	Numerator[left]*(prod/dl);
		right =	Numerator[right]*(prod/dr);
		left =	Expand[Together[left]];
		right =	Expand[Together[right]];
		If[showproof,
	 		pr = Factor[prod];
	   		Print["Multiplying both sides by ",pr," gives:"];
	   		Print[op[left,right]];
	   	  ];
	  ];
	];

gcdize[] :=
(* Implicit arguments: left, right, op, showproof. *)
Block[ {temp},

	If[left==0 || right==0, Return[]];

	temp = PolynomialGCD[left,right];

	If[NumberQ[temp] && temp<0,
 		Print["Oops: gcd was negative."];
	  ];

 	(* Should really confirm that temp is positive. *)

 	If[temp==1,Return[]];

	left =	Expand[Together[left/temp]];
	right =	Expand[Together[right/temp]];
	If[showproof,
	    temp = Factor[temp];
		Print["Dividing both sides by ", temp, " gives:"];
		Print[op[left,right]];
	    ];
	If[!NumberQ[temp],collect[]];
	];

normalorder[] :=
(* Returns left, right and op making
   left side larger than right. *)
(* Implicit inputs: left, right, op *)
Block[{},

	If[op===Less || op===LessEqual,
    	{left,right} = {right,left};
		If[op===Less,op=Greater];
		If[op===LessEqual,op=GreaterEqual];
	  ];

 	];

checkobvious[] :=
 (* Implicit inputs: left, right, op, showreason. *)
 Block[ {},

 collect[];

 If[right==0 && op=!=Equal,
    If[showreason,
       Print["This is true because a sum of positive terms must be positive."];
	  ];
    Return[True];
   ];

 If[left==0 && right=!=0 && op=!=Equal,
    If[showreason,
      Print["This is false because a sum of positive terms must be positive."];
      ];
     Return[False];
   ];

 Return[undecided];

	];