(*:Copyright: Copyright 1988-1999 Wolfram Research, Inc. *)

(*:Mathematica Version: 3.0 *)

(*:Package Version: 1.0 *)

(*:Name: Algebra`AlgebraicInequalities` *)

(*:Title: Generic Cylindrical Decomposition *) 

(*:Author: Adam Strzebonski *)

(*:Keywords: algebraic inequalities, cylindrical decomposition *)

(* :Sources: A. Strzebonski, An Algorithm for Systems of Strong Polynomial
             Inequalities, The Mathematica Journal, Vol.4 Iss. 4, Fall 1994,
	     74-77. 
*)

(*:Summary:
This package defines SemialgebraicComponents, which gives at least one point in
each connected component of an open semialgebraic set.
*)

BeginPackage["Algebra`AlgebraicInequalities`"]

SemialgebraicComponents::usage = "SemialgebraicComponents[{I1,...,In},{x1,...,xm}] gives at least one point in each connected component of the open semialgebraic set defined by inequalities {I1,...,In}. Both sides of each inequality are polynomials in variables {x1,...,xm} with rational coefficients."

Begin["Algebra`AlgebraicInequalities`Private`"]

Unprotect[SemialgebraicComponents]

allnonzero[fs_,vars_,aa_] :=
           Module[{a,i,bb,gs,sols={}},
	          Do[a=aa[[i]];
		     If[Head[a]===List,
		        gs=fs/.vars[[1]]->a[[1]];
			bb=allnonzero[gs,Rest[vars],a[[2]]];
			a=a[[1]],
			gs=fs/.vars[[1]]->a;
			bb=ratpoints[gs]];
		     sols=Append[sols,{a,bb}],
		     {i,Length[aa]}];
		  sols]

allpos[fs_,vars_,aa_] :=
       Module[{a,i,bb,gs,sols={}},
	      Do[a=aa[[i]];
		 If[Head[a]===List,
		    gs=fs/.vars[[1]]->a[[1]];
		    bb=allpos[gs,Rest[vars],a[[2]]];
		    sols=Union[sols,Prepend[#,a[[1]]]&/@bb],
		    gs=fs/.vars[[1]]->a;
		    If[Select[gs,#<=0&]=={}, sols=Append[sols,{a}]]],
		 {i,Length[aa]}];
	      sols]
			
ratpoints[fs_] := (*fs relatively prime and squarefree*)
       Module[{i,l,rr,a,aa,bb,fl=0},
              rr=System`Private`RealRoots[fs][[1]];  
              l=Length[rr]; 
              Do[If[rr[[i]][[1]]==rr[[i]][[2]], 
                    a=rr[[i]][[1]]; fl=1;
                    If[i==1, aa=a-1, aa=(rr[[i-1]][[1]]+7a)/8];
                    If[i==l, bb=a+1, bb=(rr[[i+1]][[2]]+7a)/8];
                    While[Length[System`Private`RealRoots[fs,
			                                {{aa,bb}}][[2]]]>1,
                          aa=(aa+7a)/8; bb=(bb+7a)/8];
                    rr=ReplacePart[rr,{aa,bb},i]],
                 {i,l}];
              If[fl==1, 
                 Do[If[rr[[i]][[2]]>rr[[i+1]][[1]],
                       a=(rr[[i+1]][[1]]+rr[[i]][[2]])/2;
                       rr=ReplacePart[rr,{rr[[i]][[1]],a},i];
                       rr=ReplacePart[rr,{a,rr[[i+1]][[2]]},i+1]],
                    {i,l-1}]];
	      If[l>0,
	         aa=Reverse[Append[Table[rr[[i]][[1]],{i,l}],rr[[l]][[2]]]],
		 aa={0}];
              aa]

sfrp[fs_] :=
    Module[{gs,hs,ks,i,j,f,g,d},
	   hs=Select[fs,!NumberQ[#]&];
	   hs=Transpose[FactorSquareFreeList[#]][[1]]&/@hs;
	   hs=Select[Union@@hs,!NumberQ[#]&];
	   gs={}; 
	   While[Length[hs]>0, 
	         ks={hs[[1]]}; hs=Rest[hs];
		 Do[i=1;
		    While[i<=Length[ks];
		          d=PolynomialGCD[ks[[i]],hs[[j]]];
	                  If[!NumberQ[d],
		             f=Cancel[ks[[i]]/d];
			     g=Cancel[hs[[j]]/d];
		             hs=ReplacePart[hs,g,j];
			     ks=Prepend[ReplacePart[ks,f,i],d];
			     i++];
			  i++],
	            {j,Length[hs]}];
		 gs=Union[gs,Select[ks,!NumberQ[#]&]];
                 hs=Select[hs,!NumberQ[#]&]];
	   gs]

samplepts[fs_,vars_] :=
        Module[{gs,hs={},l,i,j,aa,y,g,n=Length[vars]},
	       gs=sfrp[fs]; 
	       If[n==1,
	          ratpoints[gs],
		  y=vars[[n]];
		  l=Length[gs];
		  Do[Do[If[i==j,
		          hs=Append[hs,Last[CoefficientList[gs[[i]],y]]];
		          g=Resultant[gs[[i]],D[gs[[i]],y],y],
			  g=Resultant[gs[[i]],gs[[j]],y]];
		       hs=Append[hs,g],
		       {j,i,l}],{i,l}];
		  aa=samplepts[hs,Drop[vars,-1]]; 
		  allnonzero[gs,vars,aa]]]

semialg[fs_,vars_] :=
	Module[{gs,n=Length[vars]},
	       If[Select[fs,#<=0&]!={}, Return[{}]];
	       gs=Select[fs,!NumberQ[#]&];
	       If[gs=={}, Return[Table[0,{i,n}]]];
	       allpos[gs,vars,samplepts[gs,vars]]]

rationalpolyQ[expr_, vars_] := 
  Switch[Head[expr],
         Plus | Times, 
         And@@(rationalpolyQ[#, vars]&/@(List@@expr)),
	 Power,
	 Length[expr]==2 && IntegerQ[expr[[2]]] && Positive[expr[[2]]] &&
	 rationalpolyQ[expr[[1]], vars],
	 Integer | Rational,
	 True,
	 _,
	 MemberQ[vars, expr, 1]]
	       
SemialgebraicComponents[ineqs_,vars_]/;(Head[ineqs]===List) &&
    (Head[vars]===List) && (Variables[vars]===Sort[vars]) &&
    And@@(((Head[#]===Greater || Head[#]===Less) && Length[#]==2)&/@ineqs) &&
    And@@((rationalpolyQ[#[[1]], vars]===True && 
           rationalpolyQ[#[[2]], vars]===True)&/@ineqs) := 
    Sort[semialg[If[Head[#]===Less,#[[2]]-#[[1]],#[[1]]-#[[2]]]&/@ineqs,vars]]  

SemialgebraicComponents[ineq_,vars_]/;(Head[ineq]===Greater || 
    Head[ineq]===Less) && (Head[vars]===List) :=
    SemialgebraicComponents[{ineq},vars]
    
SemialgebraicComponents[ineqs_,var_]/;(Head[var]=!=List) :=
    Flatten[SemialgebraicComponents[ineqs,{var}]]

End[]  (* Algebra`AlgebraicInequalities`Private` *)

Protect[SemialgebraicComponents]

EndPackage[]  (* Algebra`AlgebraicInequalities` *)

            
