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

(*:Mathematica Version: 3.0 *)

(*:Package Version: 1.0 *)

(*:Name: Algebra`InequalitySolve` *)

(*:Title: Solving of univariate inequalities *) 

(*:Author: Adam Strzebonski *)

(*:Keywords: inequalities *)

BeginPackage["Algebra`InequalitySolve`"]

InequalitySolve::usage = "InequalitySolve[expr, x] gives the solution set of an expression containing logical connectives and univariate polynomial equations and inequalities in the variable x. InequalitySolve[expr, {x1, ..., xn}] gives the solution set of an expression containing logical connectives and linear equations and inequalities in the variables {x1, ..., xn}."

InequalitySolve::nineq = "`1` is not a formula constructed with univariate polynomial equations and inequalities in `2`."

InequalitySolve::nlin = "`1` is not a formula constructed with equations and inequalities in `2`, which are linear or are polynomial and have rational number coefficients."

InequalitySolve::npi = "A nonpolynomial equation or inequality encountered. The solution set may be incorrect."

InequalitySolve::mepr = "In attempting to compare solutions `1` $MaxExtraPrecision `2` was encountered. The solutions were assumed equal. Increasing the value of $MaxExtraPrecision may help resolve the uncertainty."

InequalitySolve::real = "In attempting decide whether a solution `1` is real $MaxExtraPrecision `2` was encountered. The solution was assumed real. Increasing the value of $MaxExtraPrecision may help resolve the uncertainty."

InequalitySolve::nind = "Variables `1` and `2` are not independent."

Begin["Algebra`InequalitySolve`Private`"]

Unprotect[InequalitySolve]

mymax[a___] := 
   Block[{m=Max[a], t},
         If[Head[m]=!=Max, 
            m, 
            t=m[[1]]; 
            m=Expand[#-t]&/@m;
            If[Head[m]=!=Max, Expand[m+t], Expand[#+t]&/@m]]]

mymin[a___] := 
   Block[{m=Min[a], t},
         If[Head[m]=!=Min, 
            m, 
            t=m[[1]]; 
            m=Expand[#-t]&/@m;
            If[Head[m]=!=Min, Expand[m+t], Expand[#+t]&/@m]]]
            
(* Takes formula constructed with x>a, x>=a, x<a, x<=a, x==a,  *)
(* and logical operations, and returns a simplified version.   *)
(* If the simplified version is an interval or a point returns *)
(* {a, b, ineq1, ineq2} for interval a ineq1 x ineq2 b.        *)
simpand[form_, x_] := 
   Module[{lform, gr, le, grh, leh, eq, tmp, sol},
      If[Head[form]=!=And, 
         tmp=Head[form];
         If[tmp===Equal, 
            sol={form[[2]], form[[2]], Equal, Equal},
            If[tmp===Less || tmp===LessEqual,
               sol={-Infinity, form[[2]], Less, tmp},
               If[tmp===Greater,
                  sol={form[[2]], Infinity, Less, Less},
                  sol={form[[2]], Infinity, LessEqual, Less}]]];
         Return[sol]];
      lform=List@@form;
      eq=Select[lform, Head[#]===Equal&];
      If[Length[eq]>0,
         lform=Complement[lform, eq];
         eq=(#[[2]])&/@eq;
         If[Length[eq]>1, 
            sol=Equal@@eq;
            If[sol===False, Return[False]];
            If[sol===True, eq={eq[[1]]}, eq=List@@sol]];
         sol=And@@(Head[#][eq[[1]]-#[[2]], 0]&/@lform);
         If[sol===False, Return[False]];
         If[sol===True,
            If[Length[eq]==1,
               Return[{eq[[1]], eq[[1]], Equal, Equal}],
               Return[And@@(x==#&/@eq)]]];
         sol=And@@Select[lform, Head[#][eq[[1]]-#[[2]], 0]=!=True&]; 
         sol=sol && And@@(x==#&/@eq);
         Return[sol]]; 
      gr=Select[lform, Head[#]===Greater || Head[#]===GreaterEqual&];
      le=Select[lform, Head[#]===Less || Head[#]===LessEqual&];
      gr=(#[[2]])&/@gr;
      le=(#[[2]])&/@le;
      gr=mymax[gr];
      le=mymin[le];
      If[Head[gr]=!=Max && Head[le]=!=Min,
         grh=Head/@Select[lform, #[[2]]==gr&];
         leh=Head/@Select[lform, #[[2]]==le&];
         If[FreeQ[grh, Greater], grh=LessEqual, grh=Less];
         If[FreeQ[leh, Less], leh=LessEqual, leh=Less];
         If[gr-le>0, Return[False]];
         If[Expand[gr-le]==0, 
            If[leh===Less || grh===Less, 
               Return[False],
               Return[{gr, gr, Equal, Equal}]]];
         Return[{gr, le, grh, leh}]];
      If[Head[gr]===Max, gr=List@@gr, gr={gr}];
      If[Head[le]===Min, le=List@@le, le={le}];
      If[gr[[1]]-le[[1]]>0, Return[False]];
      gr=Select[lform, MemberQ[gr, #[[2]]]&];
      le=Select[lform, MemberQ[le, #[[2]]]&];
      And@@gr && And@@le]

listtoint[int_, x_] := 
   If[int[[3]]===Equal, 
      x==int[[1]],
      If[int[[1]]===-Infinity, 
         If[int[[2]]===Infinity, True, int[[4]][x, int[[2]]]],
         If[int[[2]]===Infinity, 
            If[int[[3]]===Less, x>int[[1]], x>=int[[1]]],
            Inequality[int[[1]], int[[3]], x, int[[4]], int[[2]]]]]]

rev[h_] := Switch[h,
                  Less, Greater,
                  LessEqual, GreaterEqual,
                  Greater, Less,
                  GreaterEqual, LessEqual,
                  _, h]
                  
(* form is constructed with x>a, x>=a, x<a, x<=a, x==a, x!=a, and *)
(* logical operations. Returns simplified form.                   *)

solsimp[form_, x_] := 
   Module[{eform, sol, int, nxt, flag, tmp},
      eform=LogicalExpand[form];
      eform=LogicalExpand[eform/.Unequal[a_, b_] :> (a<b || a>b)]; 
      eform=eform/.h_[a_, x]:>rev[h][x, a];
      If[Head[eform]===And,
         eform=simpand[eform, x];
         If[Head[eform]=!=List, Return[eform]];
         Return[listtoint[eform, x]]];
      If[Head[eform]=!=Or, Return[eform]];
      eform=simpand[#, x]&/@(List@@eform);
      sol=Select[eform, Head[#]=!=List&];
      eform=Complement[eform, sol];
      sol=Or@@sol;
      If[Length[eform]==0, Return[sol]];
      tmp=Select[eform, NumericQ[#[[1]]] || Head[#[[1]]]===DirectedInfinity&];
      eform=Complement[eform, tmp];
      tmp=Sort[tmp, #1[[1]]<=#2[[1]]&];
      eform=Sort[eform, #1[[1]]-#2[[1]]<=0&];
      eform=Join[tmp, eform];
      int=eform[[1]]; eform=Rest[eform];
      While[Length[eform]>0,
         nxt=eform[[1]]; 
         eform=Rest[eform];
         flag=True;
         If[Expand[nxt[[1]]-int[[1]]]==0 && int[[3]]=!=LessEqual && 
            nxt[[3]]=!=int[[3]], int[[3]]=LessEqual];
         If[NumericQ[nxt[[1]]] || Head[int[[2]]]=!=DirectedInfinity,
            tmp=Expand[nxt[[1]]-int[[2]]],
            tmp=1];
         If[tmp<0 || tmp==0 && (int[[4]]=!=Less || nxt[[3]]=!=Less),
            tmp=Expand[nxt[[2]]-int[[2]]];
            If[tmp>0,
               int[[2]]=nxt[[2]];
               int[[4]]=nxt[[4]];
               flag=False,
               If[tmp==0,
                  flag=False;
                  If[int[[4]]=!=LessEqual && nxt[[4]]=!=int[[4]],
                     int[[4]]=LessEqual],
                  If[tmp<0, flag=False]]]];
         If[flag, sol=sol || listtoint[int, x]; int=nxt]];
      sol=sol || listtoint[int, x];
      sol]
                
realQ[x_] := 
   Block[{nr},
      If[!NumericQ[x], Return[False]];
      If[Head[x]===Root && Length[x]==3, 
         FreeQ[x, Complex], 
         If[FreeQ[x, Root] && Im[x]===0,
            True,
            nr=N[x, 20];
            If[!NumberQ[nr],
               False,
               If[FreeQ[nr, Complex],
                  True,
                  If[Im[nr]!=0,
                     False,
                     If[Im[N[x, 20+$MaxExtraPrecision]]!=0,
                        False,
                        Message[InequalitySolve::real, x, $MaxExtraPrecision];
                        True]]]]]]]

positiveQ[x_] := 
   Block[{m=$MaxExtraPrecision, tmp, nr, pr},
         $MaxExtraPrecision=$MaxPrecision;
         tmp=If[x>0, True, False, False];                      
         $MaxExtraPrecision=m;
         tmp]

sortsol[sol_] :=
   Module[{n, ss, ns, pr=8, flag=True, rep, mult},
      n=Length[sol];
      While[Length[ss]!=n && flag,
            pr=2 pr;
            If[pr>16+$MaxExtraPrecision, flag=False; pr=16+$MaxExtraPrecision];
            ns={N[#, pr], #}&/@sol;
            ss=Union[ns, SameTest->(#1[[1]]==#2[[1]]&)];
            ];
      ss=Sort[ss];
      If[Length[ss]<n,
         rep=Function[t, Select[ns, #[[1]]==t[[1]]&]]/@ss;
         mult=Length/@rep;
         rep=Select[rep, Length[#]>1&];
         Message[InequalitySolve::mepr, #, $MaxExtraPrecision]&/@rep,
         mult=Table[1, {n}]];    
      {#[[2]]&/@ss, mult}]      

univariateQ[f_, x_] := 
   If[f===x || NumericQ[f],
      True,
      If[MemberQ[Attributes[Evaluate[Head[f]]], NumericFunction],
         And@@(univariateQ[#, x]&/@(List@@f)),
         False]]

(* rts are solutions of f==0. cond is a condition in the listtoint[] *)
(* form describing a part of the domain of f. Returns inequalities   *)
(* describing head[f, 0] && cond. head is one of >, >=, <, <=.       *)

condsol[f_, x_, head_, rts_, cond_] :=
   Module[{a, b, c, n, sol, tmp, pos, neg, zer, mflag},
      Switch[Head[cond],
         Symbol, 
            If[cond=!=True, Return[$Failed]];
            a=-Infinity; b=Infinity,
         Inequality,
            If[Length[cond]!=5 || cond[[3]]=!=x, Return[$Failed]];
            a=cond[[1]]; b=cond[[5]],
         Greater|GreaterEqual, 
            If[Length[cond]!=2 || cond[[1]]=!=x, Return[$Failed]];
            a=cond[[2]]; b=Infinity,
         Less|LessEqual, 
            If[Length[cond]!=2 || cond[[1]]=!=x, Return[$Failed]];
            a=-Infinity; b=cond[[2]],
         Equal,
            If[Length[cond]!=2 || cond[[1]]=!=x, Return[$Failed]];
            Return[head[f/.x->cond[[2]], 0] && cond],
         _, Return[$Failed]]; 
      sol=Select[rts, (cond/.x->#)=!=False&];                  
      n=Length[sol];      
      If[n==0,
         c=If[a===-Infinity, 
              If[b===Infinity, 0, b-1],
              If[b===Infinity, a+1, (a+b)/2]];
         tmp=f/.x->c;
         If[head===Greater || head===GreaterEqual,                
            If[positiveQ[tmp], Return[cond], Return[False]]];
         If[head===Less || head===LessEqual,                
            If[positiveQ[tmp], Return[False], Return[cond]]];
         If[positiveQ[tmp], 
            Return[{False, False, cond}], 
            Return[{cond, False, False}]]];
      sol=sortsol[sol][[1]];
      If[Length[sol]<n, n=Length[sol]; mflag=True, mflag=False];
      pos=False; 
      neg=False;
      If[(a==sol[[1]])===False,
         c=If[a===-Infinity, sol[[1]]-1, (a+sol[[1]])/2];
         tmp=f/.x->c;
         If[positiveQ[tmp], 
            Switch[head,
               Greater, pos=pos || x<sol[[1]],
               GreaterEqual, pos=pos || x<=sol[[1]],
               List, pos=pos || x<sol[[1]]],
            Switch[head,
               Less, pos=pos || x<sol[[1]],
               LessEqual, pos=pos || x<=sol[[1]],
               List, neg=neg || x<sol[[1]]]]];
      If[(b==sol[[n]])===False,
         c=If[b===Infinity, sol[[n]]+1, (sol[[n]]+b)/2];               
         tmp=f/.x->c;
         If[positiveQ[tmp], 
            Switch[head,
               Greater, pos=pos || x>sol[[n]],
               GreaterEqual, pos=pos || x>=sol[[n]],
               List, pos=pos || x>sol[[n]]],
            Switch[head,
               Less, pos=pos || x>sol[[n]],
               LessEqual, pos=pos || x>=sol[[n]],
               List, neg=neg || x>sol[[n]]]]];
      Do[tmp=f/.x->(sol[[i]]+sol[[i+1]])/2;
         If[positiveQ[tmp], 
            Switch[head,
               Greater, pos=pos || x>sol[[i]] && x<sol[[i+1]],
               GreaterEqual, pos=pos || x>=sol[[i]] && x<=sol[[i+1]],
               List, pos=pos || x>sol[[i]] && x<sol[[i+1]]],
            Switch[head,
               Less, pos=pos || x>sol[[i]] && x<sol[[i+1]],
               LessEqual, pos=pos || x>=sol[[i]] && x<=sol[[i+1]],
               List, neg=neg || x>sol[[i]] && x<sol[[i+1]]]],
         {i, n-1}];
      If[head===List,
         zer=Or@@((x==#)&/@sol);
         {neg && cond, zer, pos && cond},
         If[mflag && (head===LessEqual || head===GreaterEqual),
            pos=pos || Or@@((x==#)&/@sol)];
         pos && cond]]     

(* f is an expression with no denominator.  *)
(* Gives conditions on f to be real valued. *)

domain[f_] := 
   Switch[Head[f],
      Alternatives[Plus, Times, Sin, Cos, Tan, Cot, Csc, Sec, Sinh, Cosh, 
         Tanh, Coth, Csch, Sech, ArcTan, ArcCot, ArcSinh, ArcCsch],
         And@@(domain/@(List@@f)),
      Power, 
         With[{a=f[[1]], b=Drop[f, 1]},
            If[IntegerQ[b], 
               domain[a],
               If[Head[b]===Rational,
                  a>=0,
                  a>0 && domain[b]]]],
      Log, 
         And@@(#>0&/@(List@@f)),
      ArcSin|ArcCos,
         If[Length[f]==1, -1<=f[[1]]<=1, True],
      ArcSec|ArcCsc,
         If[Length[f]==1, f[[1]]<=-1 || f[[1]]>=1, True],
      ArcCosh,
         If[Length[f]==1, f[[1]]>=1, True],
      ArcSech,
         If[Length[f]==1, 0<f[[1]]<=1, True],
      ArcTanh,
         If[Length[f]==1, -1<f[[1]]<1, True],
      ArcCoth,
         If[Length[f]==1, f[[1]]<-1 || f[[1]]>1, True],         
      _, True]
                                  
oneineq[ineq_, x_] :=
   Module[{f, head, pos, neg, zer, sol, cond, mess, opt, tmp, 
           i, n, lc, even, odd, mflag},
          If[ineq===True || ineq===False, Return[ineq]];
          head=Head[ineq];
          If[head=!=List && head=!=Equal && head=!=Unequal && head=!=Less &&
             head=!=LessEqual && head=!=Greater && head=!=GreaterEqual,
             Return[$Failed]];
          f=Expand[ineq[[1]]-ineq[[2]]];
          If[NumericQ[f], 
             If[head===List, 
                Return[{f<0, f==0, f>0}],
                Return[head[f, 0]]]];
          If[FreeQ[f, x], Return[$Failed]];
          tmp=Select[Variables[f], 
                     Length[#]==1 && (Head[#]===Abs || Head[#]===Sign)&];
          If[Length[tmp]>0,
             tmp=tmp[[1]];
             If[Head[tmp]===Abs,
                pos=oneineq[head[f/.tmp->tmp[[1]], 0], x];
                neg=oneineq[head[f/.tmp->-tmp[[1]], 0], x];
                cond=oneineq[{tmp[[1]], 0}, x];
                If[pos===$Failed || neg===$Failed || cond===$Failed, 
                   Return[$Failed]];
                If[head===List,
                   Return[{neg[[1]] && cond[[1]] || pos[[1]] && cond[[2]] ||
                           pos[[1]] && cond[[3]], neg[[2]] && cond[[1]] || 
                           pos[[2]] && cond[[2]] || pos[[2]] && cond[[3]],
                           neg[[3]] && cond[[1]] || pos[[3]] && cond[[2]] ||
                           pos[[3]] && cond[[3]]}],
                   Return[neg && cond[[1]] || pos && cond[[2]] || 
                          pos && cond[[3]]]],
                pos=oneineq[head[f/.tmp->1, 0], x];
                neg=oneineq[head[f/.tmp->-1, 0], x];
                zer=oneineq[head[f/.tmp->0, 0], x];
                cond=oneineq[{tmp[[1]], 0}, x];
                If[pos===$Failed || neg===$Failed || cond===$Failed || 
                   zer===$Failed, Return[$Failed]];
                If[head===List,
                   Return[{neg[[1]] && cond[[1]] || zer[[1]] && cond[[2]] ||
                           pos[[1]] && cond[[3]], neg[[2]] && cond[[1]] || 
                           zer[[2]] && cond[[2]] || pos[[2]] && cond[[3]],
                           neg[[3]] && cond[[1]] || zer[[3]] && cond[[2]] ||
                           pos[[3]] && cond[[3]]}],
                   Return[neg && cond[[1]] || zer && cond[[2]] || 
                          pos && cond[[3]]]]]];
          If[!PolynomialQ[f, x],
             f=Together[f];
             If[!NumericQ[Denominator[f]],
                sol=oneineq[head[Numerator[f] Denominator[f], 0], x];
                cond=oneineq[Denominator[f]!=0, x];
                If[sol===$Failed || cond===$Failed, Return[$Failed]];
                If[head===List,
                   Return[(# && cond)&/@sol],
                   Return[sol && cond]]]];
          lc=FactorList[f];        
          tmp=Select[lc, !FreeQ[#, x]&];
          If[Max[#[[2]]&/@tmp]>1,
             even=Select[tmp, EvenQ[#[[2]]]&];
             odd=Times@@(#[[1]]&/@Complement[lc, even]);
             even=Times@@(#[[1]]&/@even);
             If[head===Equal || head===Unequal,
                Return[oneineq[head[even*odd, 0], x]]];
             odd=oneineq[head[odd, 0], x];
             even=oneineq[even==0, x];
             If[head===Greater || head===Less, Return[odd && Not[even]]];
             If[head===GreaterEqual || head===LessEqual, 
                Return[odd || even]];
             neg=odd[[1]] && Not[even];
             zer=odd[[2]] || even;
             pos=odd[[3]] && Not[even]; 
             Return[{neg, zer, pos}]];
          mess=$Messages;
          $Messages={};
          tmp=$MessageList;
          opt=Options[Roots];
          SetOptions[Roots, Cubics -> False, Eliminate -> False, 
                     EquatedTo -> Null, Modulus -> 0, Multiplicity -> 1, 
                     Quartics -> False, Using -> True];
          sol=Solve[f==0, x];
          SetOptions@@Prepend[opt, Roots];
          $Messages=mess;
          If[Head[sol]=!=List, Return[$Failed]];
          tmp=Complement[$MessageList, tmp];
          If[FreeQ[$MessageList, HoldPattern[InequalitySolve::npi]] &&
             (Length[tmp]>0 || !PolynomialQ[f, x]), 
             Message[InequalitySolve::npi]];
          If[sol==={{}},
             If[head===Equal || head===GreaterEqual || head===LessEqual, 
                Return[True], 
                If[head===List, Return[{False, True, False}], Return[False]]],
             If[sol=!={}, sol=x/.sol]];
          If[univariateQ[f, x],
             sol=Union[Select[sol, realQ]];
             If[head===Equal, Return[Or@@((x==#)&/@sol)]];
             If[head===Unequal, Return[And@@((x!=#)&/@sol)]];
             If[!PolynomialQ[f, x], 
                cond=univariateisolve[domain[f], x];
                If[cond===$Failed, Return[$Failed]],
                cond=True];
             Return[If[Head[cond]===Or, 
                       condsol[f, x, head, sol, #]&/@cond, 
                       condsol[f, x, head, sol, cond]]]];                
          If[!PolynomialQ[f, x], Return[$Failed]];
          lc=Coefficient[f, x, Exponent[f, x]];
          If[!realQ[lc] || (lc==0)=!=False, Return[$Failed]];
          n=Length[sol];
          sol=Union[sol];
          If[Length[sol]<n, Return[$Failed]];
          If[n==0,
             If[head===Greater || head===GreaterEqual,                
                If[positiveQ[lc], Return[True], Return[False]]];
             If[head===Less || head===LessEqual,                
                If[positiveQ[lc], Return[False], Return[True]]];
             If[head===List,
                If[positiveQ[lc], 
                   Return[{False, False, True}], 
                   Return[{True, False, False}]]];
             If[head===Equal, Return[False], Return[True]]];
          tmp=sol[[1]];
          sol=Expand/@(sol-tmp);
          If[!And@@(realQ/@sol), Return[$Failed]];
          If[head===Equal, 
             sol=Expand/@(sol+tmp); 
             Return[Or@@((x==#)&/@sol)]];
          If[head===Unequal, 
             sol=Expand/@(sol+tmp); 
             Return[And@@((x!=#)&/@sol)]];
          sol=sortsol[sol];
          zer=Or@@((x==Expand[#])&/@(sol[[1]]+tmp));
          If[Length[sol[[1]]]<n, 
             sol=Transpose[sol];
             sol=Select[sol, OddQ[#[[2]]]&];
             n=Length[sol]; 
             If[n==0, 
                If[head===Greater || head===GreaterEqual,                
                   If[positiveQ[lc], 
                      Return[True], 
                      If[head===Greater,
                         Return[False],
                         Return[zer]]]];
                If[head===Less || head===LessEqual,                
                   If[positiveQ[lc], 
                      If[head===Less,
                         Return[False],
                         Return[zer]],
                         Return[True]]];
                If[head===List,
                   If[positiveQ[lc], 
                      Return[{False, zer, True}], 
                      Return[{True, zer, False}]]]];
             sol=#[[1]]&/@sol;
             mflag=True,
             sol=sol[[1]];
             mflag=False];
          sol=Expand/@(sol+tmp);
          pos=False; 
          neg=False;
          If[positiveQ[lc], lc=1, lc=-1];
          If[lc>0, 
             Switch[head,
                    Greater, pos=pos || x>sol[[n]],
                    GreaterEqual, pos=pos || x>=sol[[n]],
                    List, pos=pos || x>sol[[n]]],
             Switch[head,
                    Less, pos=pos || x>sol[[n]],
                    LessEqual, pos=pos || x>=sol[[n]],
                    List, neg=neg || x>sol[[n]]]];
          Do[lc=-lc;
             If[lc>0, 
                Switch[head,
                       Greater, pos=pos || x>sol[[i]] && x<sol[[i+1]],
                       GreaterEqual, pos=pos || x>=sol[[i]] && x<=sol[[i+1]],
                       List, pos=pos || x>sol[[i]] && x<sol[[i+1]]],
                Switch[head,
                       Less, pos=pos || x>sol[[i]] && x<sol[[i+1]],
                       LessEqual, pos=pos || x>=sol[[i]] && x<=sol[[i+1]],
                       List, neg=neg || x>sol[[i]] && x<sol[[i+1]]]],
             {i, n-1, 1, -1}];
          lc=-lc;
          If[lc>0, 
             Switch[head,
                    Greater, pos=pos || x<sol[[1]],
                    GreaterEqual, pos=pos || x<=sol[[1]],
                    List, pos=pos || x<sol[[1]]],
             Switch[head,
                    Less, pos=pos || x<sol[[1]],
                    LessEqual, pos=pos || x<=sol[[1]],
                    List, neg=neg || x<sol[[1]]]];
          If[head===List,
             {neg, zer, pos},
             If[mflag && (head===LessEqual || head===GreaterEqual),
                pos || zer,
                pos]]]

(* Eliminates x from inequalities ineqs. *)

singleprojection[ineqs_, x_, vars_] :=
   Module[{fs, gs, gr, ge, ls, le, eq, sol},
          If[ineqs===True || ineqs===False, Return[{ineqs, ineqs}]];
          fs=If[Head[ineqs]===And, List@@ineqs, {ineqs}];
          fs=Head[#][Expand[#[[1]]-#[[2]]], 0]&/@fs;
          gs=Select[fs, FreeQ[#, x]&];
          fs=Complement[fs, gs];
          fs={CoefficientList[#[[1]], x], ineqtype[#]}&/@fs;
          fs={Expand[-#[[1, 1]]/#[[1, 2]]], 
              If[#[[1, 2]]>0, #[[2]], -#[[2]]]}&/@fs;
          eq=Select[fs, #[[2]]==0&];
          If[eq=!={},
             eq=eq[[1, 1]];
             fs=And@@Union[(typetoineq[#[[2]]][eq, #[[1]]])&/@fs, gs];
             Return[{fs, x==eq}]];
          gr=Select[fs, #[[2]]==-2&];
          ge=Select[fs, #[[2]]==-1&];
          ls=Select[fs, #[[2]]==2&];
          le=Select[fs, #[[2]]==1&];
          sol=And@@Join[(x>#[[1]]&/@ls), (x>=#[[1]]&/@le), 
              (x<=#[[1]]&/@ge), (x<#[[1]]&/@gr)];               
          fs=Union[gs, Flatten[{
             Table[le[[j, 1]]<=ge[[i, 1]], {i,Length[ge]}, {j,Length[le]}],
             Table[ls[[j, 1]]<ge[[i, 1]], {i,Length[ge]}, {j,Length[ls]}],
             Table[le[[j, 1]]<gr[[i, 1]], {i,Length[gr]}, {j,Length[le]}],
             Table[ls[[j, 1]]<gr[[i, 1]], {i,Length[gr]}, {j,Length[ls]}]}]];
          fs=reduceineqs[fs, vars];
          {fs, sol}]

(* ineqs is a list of equations inequalities. Gives a conjunction of  *)
(* ineqs with redundant inequalities eliminated.                      *)

reduceineqs[ineqs_, vars_] :=
    Catch[Module[{eq, fs, gs, hs, i, tmp, gr, ls},
          If[MemberQ[ineqs, False], Throw[False]];
          eq=Select[ineqs, Head[#]===Equal&];
          fs=Select[Complement[ineqs, eq], #=!=True&];
          fs={FactorTermsList[#[[1]]-#[[2]]], ineqtype[#]}&/@fs;
          fs=If[#[[1,1]]>0, {#[[1,2]],#[[2]]}, {#[[1,2]],-#[[2]]}]&/@fs;
          gs=Union[#[[1]]&/@fs];
          If[Length[gs]<Length[fs],
             hs={};
             Do[tmp=Select[fs, #[[1]]===gs[[i]]&];
                tmp=#[[2]]&/@tmp;
                gr=Max[tmp];
                ls=Min[tmp];
                If[gr>0 && ls<0,
                   If[Max[-ls, gr]==2, 
                      Throw[False], 
                      eq=Append[eq, gs[[i]]==0]],
                   If[gr>0, 
                      hs=Append[hs, typetoineq[-gr][-gs[[i]], 0]],
                      hs=Append[hs, typetoineq[ls][gs[[i]], 0]]]],
                {i,Length[gs]}],
             hs=If[#[[2]]>0, 
                   typetoineq[-#[[2]]][-#[[1]], 0],
                   typetoineq[#[[2]]][#[[1]], 0]]&/@fs];
          If[MemberQ[hs, False], Throw[False]];
          hs=Select[hs, #=!=True&];
          fs=And@@hs;  
          fs && (And@@eq)]]

(* inequality type number codes for singleprojection *)

ineqtype[expr_] := Switch[Head[expr],
                           Less, -2, 
                           LessEqual, -1,
                           Equal, 0,
                           GreaterEqual, 1, 
                           Greater, 2,
                           _, 3]

typetoineq[n_] := Switch[n,
                          -2, Less, 
                          -1, LessEqual,  
                          0, Equal, 
                          1, GreaterEqual,
                          2, Greater]

linearinequalityQ[ineq_, vars_] := ineqtype[ineq]=!=3 &&
    Length[ineq]==2 && linearQ[ineq[[1]], vars] && linearQ[ineq[[2]], vars]
    
linearQ[expr_, vars_] := 
  If[realQ[expr],
     True,
     If[Head[expr]===Plus, 
        And@@(linearQ[#, vars]&/@(List@@expr)),
        If[Head[expr]===Times,
           With[{v=Select[List@@expr, !realQ[#]&]},
                Length[v]==1 && MemberQ[vars, v[[1]], 1]],
           MemberQ[vars, expr, 1]]]]

(* l is a list of inequalities or conjunctions of inequalities of the form *)
(* x[i] (>, >=, ==, <=, or <) linear_function[x[j] : j<i]. The solution is *)
(* equal to the conjunction of elements of l. This function simplifies     *)
(* elements of l by using elements of l of the form x[i]==number to        *)
(* replace x[i], and using univariate InequalitySolve to simplify the      *)
(* univariate elements of l.                                               *)

simpsol[l_, x_] :=
   Module[{i, n=Length[l], rul, sol, elem},
      If[MemberQ[l, False, 1], Return[False]];
      sol=univariateisolve[l[[1]], x[[1]]];
      If[Head[sol]===Equal, rul={sol[[1]]->sol[[2]]}, rul={}];
      Do[elem=solsimp[l[[i]]/.rul, x[[i]]];
         If[Head[elem]===Equal, rul=Append[rul, elem[[1]]->elem[[2]]]];
         sol=sol && elem,
         {i, 2, n}];
      sol]
      
solvebyCAD[form_, vars_] :=
   Module[{mess, ans},
      mess=$Messages;
      $Messages=Null;
      CheckAbort[
         ans=Experimental`CylindricalAlgebraicDecomposition[form, vars];
         If[Head[ans]===Experimental`CylindricalAlgebraicDecomposition,
            ans=$Failed],
         $Messages=mess;
         Abort[]];
      $Messages=mess;
      ans]
         
multivariateisolve[form_, vars_] :=
   Module[{eform, ineqs, eqs, m, i},
      eform=solvebyCAD[form, vars];
      If[eform=!=$Failed, Return[eform]];
      m=Length[vars];
      eform=LogicalExpand[form];
      eform=LogicalExpand[eform/.Unequal[a_, b_] :> (a<b || a>b)];
      If[eform===True || eform===False, Return[eform]];
      If[Head[eform]===Or,
         eform=multivariateisolve[#, vars]&/@eform;
         If[!FreeQ[eform, $Failed], eform=$Failed];
         Return[eform]];
      If[Head[eform]===And, ineqs=List@@eform, ineqs={eform}];
      If[And@@(linearinequalityQ[#, vars]&/@ineqs)=!=True, Return[$Failed]];
      ineqs=reduceineqs[ineqs, vars];
      eform={};
      Do[ineqs=singleprojection[ineqs, vars[[m-i+1]], Take[vars, m-i]];
         eform=Prepend[eform, ineqs[[2]]];
         ineqs=ineqs[[1]], 
         {i, m-1}];
      eform=Prepend[eform, ineqs];
      eform=simpsol[eform, vars];
      If[FreeQ[eform, $Failed], eform, $Failed]]

univariateisolve[form_, x_] :=
   Module[{eform},
      eform=LogicalExpand[form];
      If[Head[eform]===Or,
         eform=univariateisolve[#, x]&/@eform,
         If[Head[eform]===And, 
            eform=oneineq[#, x]&/@eform,
            eform=oneineq[eform, x]]];
      If[FreeQ[eform, $Failed], solsimp[eform, x], $Failed]]

validvar[x_] := 
   Catch[Module[{vs, c},
      vs=If[Head[x]===List, vs=x, vs={x}];
      c=Complement[vs, Variables[vs]];
      If[Length[c]>0, 
         Message[InequalitySolve::ivar, c[[1]]]; 
         Throw[False]];
      Do[c=Select[Drop[vs, {i}], !FreeQ[#, vs[[i]]]&];
         If[Length[c]>0, 
            Message[InequalitySolve::nind, vs[[i]], c[[1]]]; 
            Throw[False]],
         {i, Length[vs]}];
      True]]
            
InequalitySolve[input_, x_ ? validvar] :=
   Module[{sol, form, flag=True},
      If[Head[input]===List, form=And@@input, form=input];
      If[x==={}, Return[form]];
      If[Head[$MaxExtraPrecision::meprec]===$Off,
         flag=False,
         Off[$MaxExtraPrecision::meprec]];
      If[Head[x]===List && Length[x]>1,
         sol=multivariateisolve[form, x];
         If[sol===$Failed, Message[InequalitySolve::nlin, form, x]],
         If[Head[x]===List && Length[x]==1,
            sol=univariateisolve[form, x[[1]]],
            sol=univariateisolve[form, x]];
         If[sol===$Failed, Message[InequalitySolve::nineq, form, x]]];
      If[flag, On[$MaxExtraPrecision::meprec]];
      sol/;(sol=!=$Failed)]

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

SetAttributes[InequalitySolve, {Protected, ReadProtected}]

EndPackage[]  (* Algebra`InequalitySolve` *)

            
