(* :Title: 	Inequalites // Mathematica 1.2 and 2.0 *)

(* :Author: 	Mark Stankus (mstankus). *)

(* :Context: 	Inequalities` *)

(* :Summary:
*)

(* :Alias:
*)

(* :Warnings: 
*)

(* :History: 
   03/07/93 Changed Convert2 to Convert2`Convert2 for clarity.
   10/18/93 Changed Converts to linearConvert
*)
BeginPackage["Inequalities`",
       "Errors`"];

Clear[SetInequalityFactBase];

SetInequalityFactBase::usage = 
     "SetInequalityFactBase[aList] allows for InequalityFactQ, \
      BoundedQ,LowerBound,and UpperBound to take one parameter \
      and aList is put in for the second parameter.";

Clear[InequalityFactQ];

InequalityFactQ::usage = 
     "InequalityFactQ[eqn,aListOfFacts] gives True if the equation \
      eqn can easily be deduced from the list of assumptions \
      aListOfFacts eqn is a True inequality, False if it is False \
      and unknown if the program algorithm cannot decide if it is \
      True of False. Right now the code only works for GreaterEqual \
      equations.";

Clear[InequalitySearch];

InequalitySearch::usage = 
     "InequalitySearch[expr,aHead,aListOfFacts] searches \
      through the list aListOfFacts to find all expressions \
      gamma such that aHead[expr,gamma]===True. For example, \
      InequalitySearch[a-b,GreaterEqual, \
      {a>=2, a<= 7, b<=9, b>= -10,b<=10}] \
      would give {-7,-8} (since 2-9 == -7 and 2-10 == -8). \
      list. See also NumericalLeafs.";

Clear[InequalityToStandardForm];

InequalityToStandardForm::usage = 
     "InequalityToStandardForm[eqn] takes an inequality eqn \
      and changes it to a standardform. Here, standard form \
      is defined in terms of the Convert2 algorithm. For \
      example, InequalityToStandardForm[x-y <=0] is y>=x \
      (becuase Convert2[x-y==0] is y->x and so the inequality \
      is rearranged so that the left hand side is y).";

Clear[NumericalLeafs];

NumericalLeafs::usage = 
     "NumericalLeafs[expr,aHead,aListOfFacts] evaluates \
      NumericalLeafs[expr,aHead,aListOfFacts,40]; \
      NumericalLeafs[expr,>=,aListOfFacts,n] tries to find \
      all expressions >= 0. Since >= is transitive, this \
      may require many calls to InequalitySearch. n gives \
      the maximum number of times NumercalLeafs calls the \
      InequalitySearch module.";

Clear[BoundedTriple];

BoundedTriple::usage = 
     "BoundedTriple[expr,aListOfFacts] is a more efficient version of \
      {BoundedQ[expr,aListOfFacts], \
       LowerBound[expr,aListOfFacts], \
       UpperBound[expr,aListOfFacts]}."; 

Clear[BoundedQ];

BoundedQ::usage = 
     "BoundedQ[expr,aListOfFacts] returns True if the program \
      can easily deduce from the list of assumptions given in \
      aListOfFacts two inequalities m<= expr and expr<=n where m and \
      n are numbers.";

Clear[LowerBound];

LowerBound::usage = 
     "LowerBound[expr,aListOfFacts] returns the maximum number n such \
      that InequalitySearch can deduce expr>=n. If no such \
      inequality is generated by InequalitySearch, \
      LowerBound[expr,aListOfFacts] is -Infinity.";

Clear[UpperBound];

UpperBound::usage = 
     "UpperBound[expr,aListOfFacts] returns the minimum number n such \
      that InequalitySearch can deduce n>=expr. If no such \
      inequality is generated by InequalitySearch, \
      UpperBound[expr,aListOfFacts] is Infinity.";

Clear[SimplifyInequality];

SimplifyInequality::usage =
    "SimplifyInequality";

Begin["`Private`"];

aListOfHeads = {Greater, GreaterEqual, Equal, Less, LessEqual};

RecordFacts[aListOfFacts_List] := 
Module[{aHead,j},
     Do[ aHead = aListOfHeads[[j]];
         factdatabase[aHead] = Select[aListOfFacts,(Head[#]===aHead)&];
     ,{j,1,Length[aListOfHeads]}];
     Return[];
];

BoundedTriple[expr_,aListOfFacts_List] := 
     ( RecordFacts[aListOfFacts]; BoundedTriple[expr]);

BoundedTriple[expr_] :=
Module[{low,high,bdd},
     low = LowerBound[expr];
     high = UpperBound[expr];
     bdd = TrueQ[And[Not[low===-Infinity],Not[high===Infinity]]];
     Return[{bdd,low,high}];
];

BoundedTriple[x___] := BadCall["BoundedTriple",x];

BoundedQ[x___] := BoundedTriple[x][[1]];
     
LowerBound[expr_,aListOfFacts_List] := 
     ( RecordFacts[aListOfFacts]; LowerBound[expr]);

LowerBound[expr_] := 
Module[{lowerlist,aListOfNumbers,result},
    lowerlist = NumericalLeafs[expr,{GreaterEqual,Greater,Equal}];
    aListOfNumbers = Select[lowerlist,NumberQ];
    If[aListOfNumbers==={}, result = -Infinity;
                          , result = Apply[Max,aListOfNumbers];
    ];
    Return[result]
];

LowerBound[x___] := BadCall["LowerBound",x];

UpperBound[expr_,aListOfFacts_List] := 
     ( RecordFacts[aListOfFacts]; UpperBound[expr]);

UpperBound[expr_] := 
Module[{upperlist,aListOfNumbers,result},
    upperlist = NumericalLeafs[expr,{LessEqual,Less,Equal}];
    aListOfNumbers = Select[upperlist,NumberQ];
    If[aListOfNumbers ==={}, result = Infinity;
                           , result = Apply[Min,aListOfNumbers];
    ];
    Return[result]
];

UpperBound[x___] := BadCall["UpperBound",x];

InequalityFactQ[aListOfQuerys_List,aListOfFacts_List] :=
     ( RecordFacts[aListOfFacts]; InequalityFactQ[aListOfQuerys];);

InequalityFactQ[aListOfQuerys_List] :=
Module[{Querys,len,j,result},
     result = True;
     If[MemberQ[aListOfQuerys,False] 
              , result = False;
                Querys = {};
              , Querys = Complement[Union[aListOfQuerys],{True}];
     ];
     len = Length[Querys];
     For[j=1,j<=len && result,j++,
         result = InequalityFactQ[Querys[[j]]];
     ];
     Return[result];
];

SetInequalityFactBase[aList_List] :=
Module[{},
     InequalitiesDefault = InequalityToStandardForm[Flatten[aList]];
     RecordFacts[InequalitiesDefault]; 
];

SetInequalityFactBase[{}];

SetInequalityFactBase[x__] := BadCall["SetInequalityFactBase",x];

InequalityFactQ[True] := True;

InequalityFactQ[False] := False;

InequalityFactQ[x_Equal] := 
Module[{first,second,result},
     first = x[[1]];
     second = x[[2]];
     result = And[InequalityFactQAux[first >= second],
                  InequalityFactQAux[first <= second]
                 ];
     Return[result]
];

InequalityFactQ[x_] := InequalityFactQAux[x];

InequalityFactQ[x___] := BadCall["InequalityFactQ",x];

InequalityFactQAux[anInequality_] := 
Module[{theInequality,difference,list,temp,result},
    theInequality = InequalityToStandardForm[anInequality];
    If[Or[theInequality===True,
          theInequality===False]
         , result = theInequality;
         , difference = theInequality[[1]] - theInequality[[2]];
           list = InequalitySearch[difference,{Head[theInequality]}];
           temp = Map[(InequalityFactQ[Head[theInequality][#,0]])&,
                      list];
           If[MemberQ[temp,True], result = True;
                                , result = unknown;
           ]; 
    ];
    Return[result]
];

InequalityFactQAux[x___] := BadCall["InequalityFactQAux in Inequalities.m",x];

InequalitySearch[aSymbol_Symbol,aListOfHeads_List] := 
Module[{correctHead,correctSymbol},
    correctHead = Map[factdatabase,aListOfHeads];
    correctHead = Flatten[correctHead];
    correctSymbol = Select[correctHead,(#[[1]]===aSymbol)&];
    Return[Map[#[[2]]&,correctSymbol]];
];

InequalitySearch[x_ + y_, aListOfHeads_List] := 
Module[{temp,temp2,result},
     temp = InequalitySearch[x,aListOfHeads];
     temp2 = InequalitySearch[y,aListOfHeads];
     result = Table[ temp[[j]] + temp2[[k]]
                    ,{j,1,Length[temp]},{k,1,Length[temp2]}];
     result = Union[Flatten[result]];
     Return[result]
];

InequalitySearch[c_?NumberQ x_Symbol,aListOfHeads_List] := 
Module[{temp,result},
     If[c > 0, temp = InequalitySearch[x, aListOfHeads];
             , temp = InequalitySearch[x,Alternate[aListOfHeads]];
             , Abort[];
     ];
     result = c temp;
     result = Union[Flatten[result]];
     Return[result]
]; 

Alternate[GreaterEqual] := LessEqual;
Alternate[LessEqual] := GreaterEqual;
Alternate[Less] := Greater;
Alternate[Greater] := Less;
Alternate[Equal] := Equal;
(*
Alternate[_] := "nothing";
*)
Alternate[aList_List] := Map[Alternate,aList];
Alternate[x___] := BadCall["Alternate",x];

InequalitySearch[c_?NumberQ,aListOfHeads_List] := 
     Flatten[Map[inequalitySearch[c,#]&,aListOfHeads]];

inequalitySearch[c_?NumberQ,GreaterEqual] := {c};
inequalitySearch[c_?NumberQ,LessEqual] := {c};
inequalitySearch[c_?NumberQ,Equal] := {c};
inequalitySearch[c_?NumberQ,Less] := {c+$MachineEpsilon};
inequalitySearch[c_?NumberQ,Greater] := {c-$MachineEpsilon};
inequalitySearch[Infinity,_] := {};
inequalitySearch[-Infinity,_] := {};

inequalitySearch[x___] := BadCall["inequalitySearch",x];
InequalitySearch[x___] := BadCall["InequalitySearch",x];
 
InequalityToStandardForm[x_List] := Map[InequalityToStandardForm,x];

InequalityToStandardForm[True] := True;

InequalityToStandardForm[False] := False;

InequalityToStandardForm[eq:aHead_[_Symbol,_?NumberQ]] := eq;

InequalityToStandardForm[-aHead_[s_Symbol,c_?NumberQ]] := 
          Alternate[aHead][s,-c];

InequalityToStandardForm[x_] :=  
Module[{temp,result},
     temp = x[[1]] - x[[2]];
     result = Head[x][temp,0];
     If[And[Not[result===True],Not[result===False]]
          , result = inequalityToStandardForm[result];
     ];
     Return[result];
];

inequalityToStandardForm[x_] := 
Module[{result,head,expr,rule,top,rest,leadcoeff},
   head = Head[x];
   expr = x[[1]] - x[[2]];
   rule = Global`linearConvert[expr==0][[1]];
   top = rule[[1]];
   rest = expr/.{top->0};
   leadcoeff = Expand[(expr - rest)/top];
(* 
     Now expr ==something top + rest
*)   
   If[Not[FreeQ[leadcoeff,top]], Abort[]];
   Which[Positive[leadcoeff], 
             result = Apply[head,{rule[[1]],rule[[2]]}];
       , Negative[leadcoeff], 
             result = Apply[Alternate[head],{rule[[1]],rule[[2]]}]
       , True, result = Apply[head,{leadcoeff top, -rest}]
   ];
   Return[result];
];  

InequalityToStandardForm[x___] := BadCall["InequalityToStandardForm",x];

NumericalLeafs[expr_,aListOfHeads_List] :=
             NumericalLeafs[expr,aListOfHeads,10];

NumericalLeafs[expr_,aListOfHeads_List,NumberOfIterations_Integer?Positive] :=
Module[{temp,j},
     temp = Union[Flatten[{expr}]];
     For[j=1,j<=NumberOfIterations && Not[ListOfNumbersQ[temp]],j++,
         temp = Map[InequalitySearch[#,aListOfHeads]&,temp];
         temp = Union[Flatten[temp]];
     ];
     Return[temp]
];

NumericalLeafs[x___] := BadCall["NumericalLeafs",x];

ListOfNumbersQ[{___?NumberQ}] := True;

ListOfNumbersQ[_] := False;

ListOfNumbersQ[x___] := BadCall["ListOfNumbersQ",x];

SimplifyInequality[aList_List] := Map[SimplifyInequality,aList];

SimplifyInequality[Global`RuleTuple[aRule_,ineq_,var_]] :=
Module[{newineq,result},
     If[False (* aRule===(0->0) *)
          , result = {};
          , newineq = simplifyInequality[ineq];
            If[MemberQ[newineq,False]
                 , result = Global`RuleTuple[0->0,{},{}];
                 , result = Global`RuleTuple[aRule,newineq,var];
            ];
     ];
     Return[result];
];

SimplifyInequality[x___] := BadCall["SimplifyInequality",x];

simplifyInequality[ineq_List] :=
Module[{aList,withEqual,withoutEqual,rules,result,m,n,k,y,z,x},
     withoutEqual = InequalityToStandardForm[Union[Flatten[ineq]]];
     withEqual = Select[withoutEqual,(Head[#]===Equal)&];
     If[Not[withEqual==={}],
        rules = Map[Global`linearConvert,withEqual];
        rules = Union[Flatten[rules]];
        withoutEqual = Union[withoutEqual//.rules];
        If[MemberQ[withoutEqual,False], Return[{False}];];
        withoutEqual = InequalityToStandardForm[withoutEqual];
        withoutEqual = withoutEqual//.{
                 {x___,n_>=m_?NumberQ,y___,n_>=k_?NumberQ,z___}->
                      {x,n>=Max[m,k],y,z},
                 {x___,m_?NumberQ>=n_,y___,k_?NumberQ>=n_,z___}->
                      {x,Min[m,k]>=n,y,z}
                               };
     ];
     If[MemberQ[withoutEqual,False], result = {False};
                                   , result = Union[withoutEqual,withEqual];
     ];
     result = Complement[result,{True}];
     Return[result];
];
     
simplifyInequality[x___] := BadCall["simplifyInequality",x];

End[];
EndPackage[]
