(* :Name: findSign *)

(* :Package Version: 1.0 *)

(* :Mathematica Version: 2.2 *)

(* :Title: Find the sign of a given polynomial *)

(* :Author: Yuji Itaya
			e-mail: itaya@alice.asahi-u.ac.jp *)

(* :Context: utility`findSign` *)

(* :Requirements: No special system required *)

(* :History: Version 0.1 March 11, 1994 
             Version 0.2 Dec 20, 1994 
             Version 0.8 Feb 1, 1995
             Version 0.9 Feb. 25, 1995
             Version 0.95 Feb. 28, 1995
             Version 0.96 July 15, 1996
             Version 1.0 August 6, 1996
*)

(* :Warning: This packages warns, when the expression 
    has the possibilty of zero division.
*)
	
(* :Copyright: (c) 1995,1996  Yuji Itaya *)

(* :Examples:
	findSign[1+p,positive->{p}]
	findSign[1+p-n,positive->{p},negative->{n},
			structureDisplay->True]
	findSign[1+p-n+x^4,positive->{p},negative->{n}]
	findSign[1+p,positive->{p},tracing->On]
	findSign[1+x^4+p-n-2x^2+y^6,positive->{p},
			negative->{n}]
	findSign[a+5+x^2,positive->{a-b,b-5}]
	findSign[a+x^2+y^2+y^4-c-2x y,positive->{a-b,b-c}]
	findSign[Exp[x]]
	findSign[Log[x],positive->{x},nonpositive->{x-1/2}]
*)

BeginPackage["utility`findSign`"]
findSign::usage =
	"findSign[polynomial,option] returns the sign of a
	 polynomial. Parameters and expressions the signs
     of which are known beforehand can be specified
	 in the option."
signQ::usage=
	"signQ[polynomial] decidies the sign of a
     polynomial."
tracing::usage = 
	"tracing is an option for findSign. With
     tracing->On, intermediate results are reported."
structureDisplay::usage=
	"structureDisplay is an option for findSign.
	 With structureDisplay->True, the structure of
     a polynomial is displayed."
positive::usage=
	"option for findSign: a list of positive parameters
     and expressions"
nonpositive::usage=
	"option for findSign: a list of nonpositive
	 parameters and expressions"
zero::usage=
	"option for findSign: a list of zero parameters and
     expressions"
nonzero::usage=
	"option for findSign: a list of nonzero parameters
	 and expressions"
negative::usage=
	"option for findSign: a list of negative parameters
	 and expressions"
nonnegative::usage=
	"option for findSign: a list of nonnegative
     parameters and expressions"

Protect[positive];
Protect[nonpositive];
Protect[zero];
Protect[nonzero];
Protect[negative];
Protect[nonnegative];

Protect[tracing];

Options[findSign]=
	{tracing->Off,structureDisplay->False,
	 positive->{},nonpositive->{},zero->{},nonzero->{},
     negative->{},nonnegative->{}}

Begin["`Private`"]

signs={"positive","nonnegative","zero","nonpositive",
		"negative","nonzero"};

allSigns=Append[signs,"indeterminate"];

SetAttributes[findSign,HoldFirst];
SetAttributes[signQ,HoldFirst];
(* Note that the attribute of "sign" is not set 
   HoldFirst; Otherwise, no pattern matching in "sign"
   is performed.*)

(* findSign preprocesses the options and calls signQ. *)
findSign[polynomial_,opt___]:=
	Module[{structureD},
		(* Note that trace, positiveSymbols,...,
           nonzeroSymbols,and level are not local
           variables in this Module. *)

		(* preprocess the options *)
		trace=tracing/.{opt}/.Options[findSign];

		structureD=structureDisplay/.
			{opt}/.Options[findSign];

		(* parameters the sign of which are known *)
		positiveSymbols=Select[positive/.
			{opt}/.Options[findSign],AtomQ];
		nonnegativeSymbols=Select[nonnegative/.
			{opt}/.Options[findSign],AtomQ];
		zeroSymbols=Select[zero/.
			{opt}/.Options[findSign],AtomQ];
		nonpositiveSymbols=Select[nonpositive/.
			{opt}/.Options[findSign],AtomQ];
		negativeSymbols=Select[negative/.
			{opt}/.Options[findSign],AtomQ];
		nonzeroSymbols=Select[nonzero/.
			{opt}/.Options[findSign],AtomQ];

		(* expressions the sign of which are known *)
		positiveExpressions=Select[positive/.
			{opt}/.Options[findSign],!(AtomQ[#])&];
		nonnegativeExpressions=Select[nonnegative/.
			{opt}/.Options[findSign],!(AtomQ[#])&];
		zeroExpressions=Select[zero/.
			{opt}/.Options[findSign],!(AtomQ[#])&];
		nonpositiveExpressions=Select[nonpositive/.
			{opt}/.Options[findSign],!(AtomQ[#])&];
		negativeExpressions=Select[negative/.
			{opt}/.Options[findSign],!(AtomQ[#])&];
		nonzeroExpressions=Select[nonzero/.
			{opt}/.Options[findSign],!(AtomQ[#])&];

		level=0;  (* level is a variable to count 
					 the depth of the sign
					 determination *)

		(* call signQ and return its returned sign *)
		If[TrueQ[structureD],
			(* if structureD is True, the sign and 
			   the structure of the polynomail is
			   returned.*)
			Return[signQ[polynomial]],
			(* if structureD is False, only the sign
			   is returned.*)
			Return[First[signQ[polynomial]]]
		];
	];

(* signQ finds the sign of a polynomial by five methods.
   signQ calls sign or signQ recursively. *)

signQ[x_]:=
	Module[{s,y,f,list,s1,s2,ss,x1,x2,xs,i,j,
			usedExpression,varInX,sol,found},
		level=level+1;
		debugPrint[level,"Evaluating the sign of ",x];
                                
		(* STEP 1: factorization method *)
		debugAllPrint[level,"STEP 1"];
		f=Factor[x];
		If[x=!=f,
			s=sign[f];
			If[MemberQ[signs,s[[1]]],
				debugPrintSign[level,s];
				level=level-1;Return[s]]
		];
                                
		(* STEP 2: Sturm's method *)
		debugAllPrint[level,"STEP 2"];
		f=sturmMethod[x];
		If[First[f]==0,
			s={sign[Last[f]][[1]],{x}};
			debugPrintSign[level,s];
			level=level-1;
			Return[s]];
                                
		(* STEP 3: the method of the sum of squares *)
		debugAllPrint[level,"STEP 3"];
		If[Head[x]===Plus,
			list=listDegreeVariables[x];
			If[Length[list]!=0,
				y=Max[Transpose[list][[2]]];
				If[EvenQ[y] && (y!=0),
					For[i=1,i<=Length[list],i++,
						{x1,x2}=
					   		polyDecomposed[x,list[[i]]];
						debugPrint[level,x,
							" is decomposed into ",
							x1," and ",x2];
						s1=signQ[x1];
						s2=signQ[x2];
						s=sign[s1[[1]]+s2[[1]]];
						ss=s[[1]];
						If[MemberQ[signs,ss],Break[]]];
					If[MemberQ[signs,ss],
						debugPrintSign[level,
							{ss,{Plus,s1,s2}}];
						level=level-1;
						Return[{ss,{Plus,s1,s2}}]]
					]
				];  
			];
                                
 
		(* STEP 4: the substitution method *)
		debugAllPrint[level,"STEP 4"];
		found=0; (* if the sign is found, found is set
				    to 1 *)
		varInX=Variables[x];

		(* eliminate the variables using the
		   positive expressions *)
		Do[usedExpression=positiveExpressions[[i]];
			sol=Solve[usedExpression==
				"positive " <> ToString[usedExpression],
									varInX[[j]]];
			If[Head[sol]=!=Solve,
				xs=(x/.sol)[[1]];
				If[xs=!=x,
					positiveExpressions=
					   Drop[positiveExpressions,{i,i}];
					s=signQ[xs];
					positiveExpressions=
						Insert[positiveExpressions,
							usedExpression,i];
					If[MemberQ[signs,s[[1]]],
						found=1;Break[]];
				];
			],
			{j,1,Length[varInX]},
			{i,1,Length[positiveExpressions]}];

		If[found==1,
			If[MemberQ[signs,s[[1]]],
				debugPrintSign[level,s];
				level=level-1;
				Return[s]
			]
		];

		(* eliminate the variables using the
		   nonnegative expressions *)
		Do[usedExpression=nonnegativeExpressions[[i]];
			sol=Solve[usedExpression==
				"nonnegative " <> 
					ToString[usedExpression],
						varInX[[j]]];
			If[Head[sol]=!=Solve,
				xs=(x/.sol)[[1]]; 
				If[xs=!=x,
					nonnegativeExpressions=
						Drop[nonnegativeExpressions,
							{i,i}];
					s=signQ[xs];
					nonnegativeExpressions=
						Insert[nonnegativeExpressions,
							usedExpression,i];
				
					If[MemberQ[signs,s[[1]]],
						found=1;Break[]];
				];
			],
			{j,1,Length[varInX]},
			{i,1,Length[nonnegativeExpressions]}];

		If[found==1,
			If[MemberQ[signs,s[[1]]],
				debugPrintSign[level,s];
				level=level-1;
				Return[s]
			]
		];

		(* eliminate the variables using the
		   zero expressions *)
		Do[usedExpression=zeroExpressions[[i]];
			sol=Solve[usedExpression==
				"zero " <> ToString[usedExpression],
								varInX[[j]]];
			If[Head[sol]=!=Solve,
				xs=(x/.sol)[[1]]; 
				If[xs=!=x,
					zeroExpressions=
						Drop[zeroExpressions,{i,i}];
					s=signQ[xs];
					zeroExpressions=
						Insert[zeroExpressions,
							usedExpression,i];
					If[MemberQ[signs,s[[1]]],
						found=1;Break[]]
				];
			],
			{j,1,Length[varInX]},
			{i,1,Length[zeroExpressions]}];

		If[found==1,
			If[MemberQ[signs,s[[1]]],
				debugPrintSign[level,s];
				level=level-1;
				Return[s]
			]
		];

		(* eliminate the variables using the
		   nonpositive expressions *)
		Do[usedExpression=nonpositiveExpressions[[i]];
			sol=Solve[usedExpression==
				"nonpositive " <> 
					ToString[usedExpression],
						varInX[[j]]];
			If[Head[sol]=!=Solve,
				xs=(x/.sol)[[1]]; 
				If[xs=!=x,
					nonpositiveExpressions=
						Drop[nonpositiveExpressions,
							{i,i}];
					s=signQ[xs];
					nonpositiveExpressions=
						Insert[nonpositiveExpressions,
							usedExpression,i];
					If[MemberQ[signs,s[[1]]],
						found=1;Break[]]
				];
			],
			{j,1,Length[varInX]},
			{i,1,Length[nonpositiveExpressions]}];

		If[found==1,
			If[MemberQ[signs,s[[1]]],
				debugPrintSign[level,s];
				level=level-1;
				Return[s]
			]
		];

		(* eliminate the variables using the
		   negative expressions *)
		Do[usedExpression=negativeExpressions[[i]];
			sol=Solve[usedExpression==
				"negative " <> ToString[usedExpression],
									varInX[[j]]];
			If[Head[sol]=!=Solve,
				xs=(x/.sol)[[1]]; 
				If[xs=!=x,
					negativeExpressions=
						Drop[negativeExpressions,
							{i,i}];
					s=signQ[xs];
					negativeExpressions=
						Insert[negativeExpressions,
							usedExpression,i];
					If[MemberQ[signs,s[[1]]],
						found=1;Break[]]
				];
			],
			{j,1,Length[varInX]},
			{i,1,Length[negativeExpressions]}];

		If[found==1,
			If[MemberQ[signs,s[[1]]],
				debugPrintSign[level,s];
				level=level-1;
				Return[s]
			]
		];

		(* eliminate the variables using the
		   nonzero expressions *)
		Do[usedExpression=nonzeroExpressions[[i]];
			sol=Solve[usedExpression==
				"nonzero " <> 
					ToString[usedExpression],
						varInX[[j]]];
			If[Head[sol]=!=Solve,
				xs=(x/.sol)[[1]]; 
				If[xs=!=x,
					nonzeroExpressions=
						Drop[nonzeroExpressions,{i,i}];
					s=signQ[xs];
					nonzeroExpressions=
						Insert[nonzeroExpressions,
							usedExpression,i];
					If[MemberQ[signs,s[[1]]],
						found=1;Break[]]
				];
			],
			{j,1,Length[varInX]},
			{i,1,Length[nonzeroExpressions]}];

		If[found==1,
			If[MemberQ[signs,s[[1]]],
				debugPrintSign[level,s];
				level=level-1;
				Return[s]
			]
		];

		(*STEP 5*)        
		debugAllPrint[level,"STEP 5"];
		s=sign[x];
		If[MemberQ[signs,s[[1]]],
			debugPrintSign[level,s];
			level=level-1;Return[s]];

		(* RESERVE FOR OTHER STEPS *)


		(* When every method fails, the sign is
		   indeterminate *)
		s={"indeterminate",{x}};
		debugPrintSign[level,s];
		level=level-1;
		Return[s]
	];


(* routine for printing messages in debugging *)
debugPrint[level_,out__]:=
	If[(trace===On||trace===All),
		Print[StringJoin[Table["  ",{level-1}]],out]];

debugAllPrint[level_,out__]:=
	If[trace===All,Print[StringJoin[Table["  ",
		{level-1}]],out]];
	
(* routine for printing the sign of the polynomial *)
debugPrintSign[level_,s_]:=debugPrint[level,
	"its sign->",s[[1]]];

(* finding the sign of the polynomial by simple
   rules *)
sign[x_?AtomQ]:=
	Module[{pos},
		If[Positive[x],Return[{"positive",{x}}]];
		If[TrueQ[x==0],Return[{"zero",{x}}]];
		If[Negative[x],Return[{"negative",{x}}]];
		If[MemberQ[positiveSymbols,x],
			Return[{"positive",{x}}]];
		If[MemberQ[nonnegativeSymbols,x],
			Return[{"nonnegative",{x}}]];
		If[MemberQ[zeroSymbols,x],Return[{"zero",{x}}]];
		If[MemberQ[nonpositiveSymbols,x],
			Return[{"nonpositive",{x}}]];
		If[MemberQ[negativeSymbols,x],
			Return[{"negative",{x}}]];
		If[MemberQ[nonzeroSymbols,x],
			Return[{"nonzero",{x}}]];
		If[MemberQ[signs,x],Return[{x,{x}}]];
		(* if x is a sign followed by an expression,
		   take its sign, and determine it *)
		pos=StringPosition[x," ",1];
		If[pos!={},
			If[MemberQ[signs,
						StringTake[x,pos[[1,1]]-1]],
					Return[{StringTake[x,pos[[1,1]]-1],
						{ToString[StringDrop[x,
							pos[[1,1]]]]}}]]];
		Return[{"indeterminate",{x}}] 
	];
                                
sign[-x_]:=
	Module[{s=signQ[x],p},
		p=Position[allSigns,s[[1]]];
		Return[Part[
			{{"negative",{Times,-1,s}},
			 {"nonpositive",{Times,-1,s}},
             {"zero",{Times,-1,s}},
	         {"nonnegative",{Times,-1,s}},
             {"positive",{Times,-1,s}},
             {"nonzero",{Times,-1,s}},
			 {"indeterminate",{Times,-1,s}}},
			  p[[1,1]] 
		  ]
		]
	];


sign[x_+y_]:=
	Module[{s1=signQ[x],s2=signQ[y],p1,p2},
		p1=Position[allSigns,s1[[1]]];
		p2=Position[allSigns,s2[[1]]];
		Return[Part[
		  {
			{{"positive",{Plus,s1,s2}},
			 {"positive",{Plus,s1,s2}},
			 {"positive",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
		 	 {"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}}},
			{{"positive",{Plus,s1,s2}},
			 {"nonnegative",{Plus,s1,s2}},
			 {"nonnegative",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}}},
			{{"positive",{Plus,s1,s2}},
			 {"nonnegative",{Plus,s1,s2}},
			 {"zero",{Plus,s1,s2}},
			 {"nonpositive",{Plus,s1,s2}},
			 {"negative",{Plus,s1,s2}},
			 {"nonzero",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}}},
			{{"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"nonpositive",{Plus,s1,s2}},
			 {"nonpositive",{Plus,s1,s2}},
			 {"negative",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}}},
			{{"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"negative",{Plus,s1,s2}},
			 {"negative",{Plus,s1,s2}},
			 {"negative",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}}},
			{{"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"nonzero",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}}},
			{{"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}},
			 {"indeterminate",{Plus,s1,s2}}}
		  },
		  p1[[1,1]],p2[[1,1]] 
		  ]
		]
	];

sign[x_ y_]:=
	Module[{s1=signQ[x],s2=signQ[y],p1,p2},
		p1=Position[allSigns,s1[[1]]];
		p2=Position[allSigns,s2[[1]]];
		Return[Part[
		  {
			{{"positive",{Times,s1,s2}},
			 {"nonnegative",{Times,s1,s2}},
			 {"zero",{Times,s1,s2}},
			 {"nonpositive",{Times,s1,s2}},
			 {"negative",{Times,s1,s2}},
			 {"nonzero",{Times,s1,s2}},
			 {"indeterminate",{Times,s1,s2}}},
			{{"nonnegative",{Times,s1,s2}},
			 {"nonnegative",{Times,s1,s2}},
			 {"zero",{Times,s1,s2}},
			 {"nonpositive",{Times,s1,s2}},
			 {"nonpositive",{Times,s1,s2}},
			 {"indeterminate",{Times,s1,s2}},
			 {"indeterminate",{Times,s1,s2}}},
			{{"zero",{Times,s1,s2}},
			 {"zero",{Times,s1,s2}},
			 {"zero",{Times,s1,s2}},
			 {"zero",{Times,s1,s2}},
			 {"zero",{Times,s1,s2}},
			 {"zero",{Times,s1,s2}},
			 {"zero",{Times,s1,s2}}},
			{{"nonpositive",{Times,s1,s2}},
			 {"nonpositive",{Times,s1,s2}},
			 {"zero",{Times,s1,s2}},
			 {"nonnegative",{Times,s1,s2}},
			 {"nonnegative",{Times,s1,s2}},
			 {"indeterminate",{Times,s1,s2}},
			 {"indeterminate",{Times,s1,s2}}},
			{{"negative",{Times,s1,s2}},
			 {"nonpositive",{Times,s1,s2}},
			 {"zero",{Times,s1,s2}},
			 {"nonnegative",{Times,s1,s2}},
			 {"positive",{Times,s1,s2}},
			 {"nonzero",{Times,s1,s2}},
			 {"indeterminate",{Times,s1,s2}}},
			{{"nonzero",{Times,s1,s2}},
			 {"indeterminate",{Times,s1,s2}},
			 {"zero",{Times,s1,s2}},
			 {"indeterminate",{Times,s1,s2}},
			 {"nonzero",{Times,s1,s2}},
			 {"nonzero",{Times,s1,s2}},
			 {"indeterminate",{Times,s1,s2}}},
			{{"indeterminate",{Times,s1,s2}},
			 {"indeterminate",{Times,s1,s2}},
			 {"indeterminate",{Times,s1,s2}},
			 {"indeterminate",{Times,s1,s2}},
			 {"indeterminate",{Times,s1,s2}},
			 {"indeterminate",{Times,s1,s2}},
			 {"indeterminate",{Times,s1,s2}}}
		  },
		  p1[[1,1]],p2[[1,1]] 
		  ]
		]
	];


sign[x_^n_]:=
	Module[{s1=signQ[x],s2=signQ[n],p1,p2},
		p1=Position[allSigns,s1[[1]]];
		p2=Position[allSigns,s2[[1]]];
		warning[Part[
		  {
			{"","","","","","",""},
			{"","","","might divide by zero",
				"might divide by zero",
				"might divide by zero",
				"might divide by zero"},
			{"","","","divide by zero, or zero by zero",
				"divide by zero",
				"might divide by zero",""},
			{"","","","might divide by zero",
				"might divide by zero",
				"might divide by zero",
				"might divide by zero"},
			{"","","","","","",""},
			{"","","","","","",""},
			{"","","","divide by zero, or zero by zero",
				"divide by zero",
				"might divide by zero",""} 		 
		  },
		  p1[[1,1]],p2[[1,1]] 
		  ]
		];		
		Return[Part[
		{
			{{"positive",{Power,s1,s2}},
			 {"positive",{Power,s1,s2}},
			 {"positive",{Power,s1,s2}},
			 {"positive",{Power,s1,s2}},
			 {"positive",{Power,s1,s2}},
			 {"positive",{Power,s1,s2}},
			 {"positive",{Power,s1,s2}}},
			{{"nonnegative",{Power,s1,s2}},
			 {"nonnegative",{Power,s1,s2}},
			 {"positive",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}}},
			{{"zero",{Power,s1,s2}},
			 {"nonnegative",{Power,s1,s2}},
			 {"positive",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}}},
			{If[IntegerQ[n],If[EvenQ[n],
						{"nonnegative",{Power,s1,s2}},
						{"nonpositive",{Power,s1,s2}}],
			 		{"indeterminate",{Power,s1,s2}}],
			  If[IntegerQ[n],If[EvenQ[n],
						{"nonnegative",{Power,s1,s2}},
						{"nonpositive",{Power,s1,s2}}],
					{"indeterminate",{Power,s1,s2}}],
			  {"positive",{Power,s1,s2}},
			  {"indeterminate",{Power,s1,s2}},
			  {"indeterminate",{Power,s1,s2}},
			  {"indeterminate",{Power,s1,s2}},
			  {"indeterminate",{Power,s1,s2}}},
			{If[IntegerQ[n],If[EvenQ[n],
						{"positive",{Power,s1,s2}},
						{"negative",{Power,s1,s2}}],
			 		{"indeterminate",{Power,s1,s2}}],
			   If[IntegerQ[n],If[EvenQ[n],
						{"nonnegative",{Power,s1,s2}},
						{"nonpositive",{Power,s1,s2}}],
					{"indeterminate",{Power,s1,s2}}],
			 {"positive",{Power,s1,s2}},
			   If[IntegerQ[n],If[EvenQ[n],
						{"positive",{Power,s1,s2}},
						{"negative",{Power,s1,s2}}],
					{"indeterminate",{Power,s1,s2}}],
			   If[IntegerQ[n],If[EvenQ[n],
						{"positive",{Power,s1,s2}},
						{"negative",{Power,s1,s2}}],
					{"indeterminate",{Power,s1,s2}}],
			 {"indeterminate",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}}},
			{{"indeterminate",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}},
			 {"positive",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}}},
			{If[IntegerQ[n],If[EvenQ[n],
					{"nonnegative",{Power,s1,s2}},
					{"indeterminate",{Power,s1,s2}}],
				{"indeterminate",{Power,s1,s2}}],
			 If[IntegerQ[n],If[EvenQ[n],
					{"nonnegative",{Power,s1,s2}}],
					{"indeterminate",{Power,s1,s2}}],
			 {"positive",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}},
			 {"indeterminate",{Power,s1,s2}}}
		  },
		  p1[[1,1]],p2[[1,1]] 
		  ]
		]
	];


(* routine for printing a warning message *)
warning[s_String]:=If[s!="",Print["Warning: ",s]];	


(* Sturm's method *)
omitZeroes[x_] := 
	Module[{m=x},
		m=Select[m,#=!="zero"&];
		m=m//.
			{a___,"nonnegative","nonnegative",b___}->
			{a,"nonnegative",b};
		m=m//.
			{a___,"nonpositive","nonpositive",b___}->
			{a,"nonpositive",b};
		m=m//.
			{a___,"nonnegative","positive",b___}->
			{a,"positive",b};
		m=m//.
			{a___,"positive","nonnegative",b___}->
			{a,"positive",b};
		m=m//.
			{a___,"negative","nonpositive",b___}->
			{a,"negative",b};
		m=m//.
			{a___,"nonpositive","negative",b___}->
			{a,"negative",b};
		Return[m];
	];

countChanges[m_] := 
	Module[{i,len=Length[m],num=0} ,
		Do[ If[m[[i]]=!= m[[i+1]], num++], {i,len-1} ];
		num
	];
	
topCoeff[f_,x_]:=Coefficient[f,x,Exponent[f,x]];

signedTopCoeff[f_,x_]:=
	Module[{d=Exponent[f,x]}, Coefficient[f,x,d] (-1)^d ];

sturmSequence[f_,g_,x_]:=
	Module[ { seq = {f}, p=f, q=g},
		While[If[q!=0,True,False,True],
			seq = Append[seq,q];
			{p,q} = {q,
				Expand[-PolynomialRemainder[p,q,x]]};
		];
		seq
	];

countRoots[f_,x_]:=
	Module[{seq = sturmSequence[f,D[f,x],x],s1,s2} ,
		s1=Map[signQ[#][[1]]&,
				Map[signedTopCoeff[#,x]&,seq]];
		s2=Map[signQ[#][[1]]&,
				Map[topCoeff[#,x]&,seq]];
		If[MemberQ[s1,"indeterminate"] || 
			MemberQ[s2,"indeterminate"], Return[-1]];
		s1=omitZeroes[s1];
		If[MemberQ[s1,"nonnegative"] ||
			MemberQ[s1,"nonpositive"],Return[-1]];
		s2=omitZeroes[s2];
		If[MemberQ[s2,"nonnegative"] ||
			MemberQ[s2,"nonpositive"],Return[-1]];
		countChanges[s1] -countChanges[s2]
	];

sturmMethod[f_]:=
	Module[{poly,symbol},
		poly=Expand[f];
		symbol=Variables[poly];
 		symbol=Complement[symbol,positiveSymbols,
				nonnegativeSymbols,	zeroSymbols,
				nonpositiveSymbols,negativeSymbols,
				nonzeroSymbols];	
		If[Length[symbol]!=1,Return[{-1,Null}]];
		If[Not[PolynomialQ[poly,symbol[[1]]]],
			Return[{-1,Null}]];
		Return[{countRoots[poly,symbol[[1]]],
			topCoeff[poly,symbol[[1]]]}];
	];

(* decompose a polynomial to a square and the rest *)
polyDecomposed[poly_,{symbol_,degree_}]:=
	Module[{coeff,polySmall},
		coeff=Coefficient[poly,symbol,degree];
		If[(degree<2)||
				(Coefficient[poly,symbol,degree-1]===0),
			Return[{Expand[poly-coeff*symbol^degree],
				coeff*symbol^degree}]];
		polySmall=coeff (symbol+
					Coefficient[poly,symbol,degree-1]/
						(degree coeff))^degree;
		Return[{Expand[poly-polySmall],polySmall}];
	];


listDegreeVariables[polynomial_]:=
	Module[{symbols,poly},
		poly=Expand[polynomial];
		symbols=Variables[poly];
		symbols=Complement[symbols,positiveSymbols,
				nonnegativeSymbols,zeroSymbols,
				nonpositiveSymbols,negativeSymbols,
				nonzeroSymbols];
		If[poly==0,Return[{0,Null}]];
		
		Return[Transpose[{symbols,Map[Exponent[poly,#]&,
				symbols]}]];
	];

(* knowledge-based *)
sign[Exp[x_]]:={"positive",{Exp,x}};

sign[Log[x_]]:=
	Module[{s=signQ[x],ss},
		If[s[[1]]=="positive",
			ss=signQ[x-1];
			Switch[ss[[1]],				
				"positive",Return[{"positive",{Log,s}}],
        								   (* x>1 *)
				"nonnegative",
					Return[{"nonnegative",{Log,s}}],
										   (* x>=1 *)
				"zero",Return[{"zero",{Log,s}}],
         				      			   (* x=1 *)
				"nonpositive",
					Return[{"nonpositive",{Log,s}}],
       									   (* 0<x<=1 *)
				"negative",Return[{"negative",{Log,s}}],
     									   (* 0<x<1 *)
				_,Return[{"indeterminate",{Log,s}}] 
			],
	   		(* Else *)
			Return[{"indeterminate",{Log,s}}]
		]
	];

sign[Log[b_,x_]]:=sign[Log[x]/Log[b]];

End[]; (* end the private context *)

EndPackage[] (* end the package context *)
