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

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

(* :Context: 	PabloInequalities` *)

(* :Summary:
*)

(* :Alias:
*)

(* :Warnings: 
*)

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

Clear[SetMyInequalityFactBase];

SetMyInequalityFactBase::usage = 
     "SetMyInequalityFactBase[aList] allows for MyInequalityFactQ, \
      MyBoundedQ,MyLowerBound,and MyUpperBound to take one parameter \
      and aList is put in for the second parameter.";

Clear[MyInequalityFactQ];

MyInequalityFactQ::usage = 
     "MyInequalityFactQ[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[MyInequalitySearch];

MyInequalitySearch::usage = 
     "MyInequalitySearch[expr,aHead,aListOfFacts] searches \
      through the list aListOfFacts to find all expressions \
      gamma such that aHead[expr,gamma]===True. For example, \
      MyInequalitySearch[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 NonNumericalLeafs.";

Clear[MyInequalityToStandardForm];

MyInequalityToStandardForm::usage = 
     "MyInequalityToStandardForm[eqn] takes an inequality eqn \
      and changes it to a standardform. Here, standard form \
      is defined in terms of the Convert2 algorithm. For \
      example, MyInequalityToStandardForm[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[NonNumericalLeafs];

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

Clear[MyBoundedTriple];

MyBoundedTriple::usage = 
     "MyBoundedTriple[expr,aListOfFacts] is a more efficient version of \
      {MyBoundedQ[expr,aListOfFacts], \
       MyLowerBound[expr,aListOfFacts], \
       MyUpperBound[expr,aListOfFacts]}."; 

Clear[MyBoundedQ];

MyBoundedQ::usage = 
     "MyBoundedQ[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[MyLowerBound];

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

Clear[MyUpperBound];

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

Clear[SimplifyMyInequality];

SimplifyMyInequality::usage =
    "SimplifyMyInequality";

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[];
];

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

MyBoundedTriple[expr_] :=
Module[{low,high,bdd},
     low = MyLowerBound[expr];
     high = MyUpperBound[expr];
     bdd = TrueQ[And[Not[low===-Infinity],Not[high===Infinity]]];
     Return[{bdd,low,high}];
];

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

MyBoundedQ[x___] := MyBoundedTriple[x][[1]];
     
MyLowerBound[expr_,aListOfFacts_List] := 
     ( RecordFacts[aListOfFacts]; MyLowerBound[expr]);

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

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

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

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

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

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

SetMyInequalityFactBase[{}];

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

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

MyInequalityFactQ[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 = MyInequalityFactQ[Querys[[j]]];
     ];
     Return[result];
];

MyInequalityFactQ[True] := True;

MyInequalityFactQ[False] := False;

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

MyInequalityFactQ[x_] := MyInequalityFactQAux[x];

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

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

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

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

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

MyInequalitySearch[c_?NumberQ x_Symbol,aListOfHeads_List] := 
Module[{temp,result},
     If[c > 0, temp = MyInequalitySearch[x, aListOfHeads];
             , temp = MyInequalitySearch[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];

MyInequalitySearch[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];
MyInequalitySearch[x___] := BadCall["MyInequalitySearch",x];
 
MyInequalityToStandardForm[x_List] := Map[MyInequalityToStandardForm,x];

MyInequalityToStandardForm[True] := True;

MyInequalityToStandardForm[False] := False;

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

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

MyInequalityToStandardForm[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];
];  

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

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

NonNumericalLeafs[expr_,aListOfHeads_List,NumberOfIterations_Integer?Positive] :=
Module[{temp,temp2,j},
     temp = Union[Flatten[{expr}]];
     temp2 = {};
     For[j=1,j<=NumberOfIterations,j++,
         temp2 = Map[MyInequalitySearch[#,aListOfHeads]&,temp];
         If[j===1,
            temp = Union[Flatten[temp2]],
            temp = Union[Flatten[temp2],temp]
           ];
     ];
     Return[temp]
];

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

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

ListOfNumbersQ[_] := False;

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

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

SimplifyMyInequality[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];
];

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

simplifyInequality[ineq_List] :=
Module[{aList,withEqual,withoutEqual,rules,result,m,n,k,y,z,x},
     withoutEqual = MyInequalityToStandardForm[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 = MyInequalityToStandardForm[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[]
