(* ::Package:: *)

(* ::Section::Closed:: *)
(*Ultra permissive "MIT copyright" notice ...*)


(* ::Text:: *)
(*Copyright (c) 2017, David R. Stoutemyer *)
(**)
(*Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:*)
(**)
(*The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.*)
(**)
(*THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.*)


(* ::Section:: *)
(*Summary ...*)


(* ::Text:: *)
(*This write-protected Mathematica package named  Integral.m  is meant to be used with a notebook file named  LoadAndUseIntegral.nb,  which should be stored in the same file directory.*)


(* ::Text:: *)
(*Both files are freely downloadable from   (to be determined).*)


(* ::Text:: *)
(*This version of this package is given by the value of the public package-level variable $IntegralVersionNumber.*)


(* ::Text:: *)
(*This package works for Mathematica version 11.0 or later.  Some of the capability is usable with earlier versions of Mathematica.*)


(* ::Subsection::Closed:: *)
(*History and purpose ...*)


(* ::Text:: *)
(*As of April 2017, this Mathematica package is under development by David R. Stoutemyer (dstout@hawaii.edu), in support of an article by Robert M. Corless, David J. Jeffrey and David R. Stoutemyer, titled "Integrals of functions containing parameters", for which we plan to post a preprint at  arXiv.org .*)
(**)
(*An extended abstract  is "Comprehensive anti-derivatives and parametric continuity", ACM Communications in Computer Algebra 52(2):32-33, DOI: 10.1145/3282678.3282680 .*)


(* ::Subsection::Closed:: *)
(*Some programming issues ...*)


(* ::Text:: *)
(*Naming conventions:*)
(*1.  All function names are capitalized.*)
(*2.  Global indeterminate read-only symbols that end users might see, such as optional-argument*)
(*        names, are merely capitalilzed.*)
(*     Otherwise global variables and dynamically-scoped variables intended for access outside*)
(*         the Block where they are defined begin with $.*)
(*3.  Such variables that are not expected to change value after "Run Package" are capitalized.*)


(* ::Section::Closed:: *)
(*Usage ...*)


BeginPackage["Integral`", {"NumericalCalculus`"}];


$IntegralVersionNumber::usage = "Please include the value of this variable in bug reports and citations."


Integral::usage =
"\[Bullet] The main benefit of  Integral [...]  is to return Comprehensive Piecewise results that include all cases depending on special values of the variables in the result.
\[Bullet] Another  optional  benefit is to combine the benefits of the Rubi antiderivative package freely downloadable from   https://rulebasedintegration.org/ with the built-in Mathematica Integrate function.
\[Bullet] Use Integral [...]  just like the built-in Mathematica  Integrate [...] function.
\[Bullet] Additional optional parameter names, their initial default values, and their alternative values are:
  Most important :
    ParameterTreatment \[Rule] Comprehensive,  {Generic, Conditional, Continuous, ApartContinuous};
    SubtaskTimeLimit \[Rule] 4.0,  {Infinity or any positive integer or floating-point real};
  Advanced :
    MostGeneralFirstQ \[Rule] True,  {False};
    TrueInPiecewise \[Rule] Automatic,  {None, All};
    NLimitUse \[Rule] None,  {Automatic, ExactNumberQ, InexactNumberQ};
    GenericIntegrationFunction \[Rule] IntThenIntegrate, {Integrate, Int};
    IntThenIntegrateUsedQ \[Rule] False,  {True};
    ShowIntegralStepsQ \[Rule] False,  {True}
Explanations with examples are in the associated file named  LoadAndUseIntegral.nb .
A download site for that file and this one is  (to be determined).";


IntThenIntegrate::usage = "If the Rubi integration package is loaded, then IntThenIntegrand [integrand, var]
first tries Rubi Int, then applies Integrate [...] to any parts that Rubi could not complete.
Otherwise IntThenIntegrate directly invokes Integrate. ";
PlotSurfaceAnd2DCurve::usage = "";


ParameterTreatment::usage = "";
SubtaskTimeLimit::usage = "";
MostGeneralFirstQ::usage = "";
TrueInPiecewiseMeans::usage = "";
NLimitUse::usage = "";
GenericIntegrationFunction::usage = "";
IntThenIntegrateUsedIntegrateQ::usage = "";
ShowIntegralStepsQ::usage = "";


Generic::usage =
"This alternative for the optional parameter  ParameterTreatment  in function  Integral [\[Ellipsis]]  returns only a generic integral.";
Conditional::usage =
"This alternative for the optional parameter  ParameterTreatment  in function  Integral [\[Ellipsis]]  returns a ConditionalExpression if needed.";
Comprehensive::usage = "";
"This alternative for the optional parameter  ParameterTreatment  in function  Integral [\[Ellipsis]]  returns a comprehensive result, using Piecewise [\[Ellipsis]] if necessary to do so.";
Continuous::usage = "";
"This alternative for the optional parameter  ParameterTreatment  in function  Integral [\[Ellipsis]]  attempts to return a Comprehensive result that is continuous with respect to the parameters therein.";
ApartContinuous::usage =
"This alternative for the optional parameter  ParameterTreatment  in function  Integral [\[Ellipsis]]  attempts to return a parametrically continuous result in which the constants of integration are separate terms.";

singularities::usage = "If [$VersionNumber < 12.0, Simplify`FunctionSingularities, FunctionProperties`Singularities]";


$DefaultGenericIntegrationFunction::usage =
 "Default value of option GenericIntegrationFunction, initially Automatic.";
$DefaultIntThenIntegrateUsedIntegrateQ::usage =
 "Default value of option IntThenIntegrateUsedIntegrateQ, initially False.";
$DefaultIntIfLoadedThenIntegrateUsedIntegrateQ::usage =
 "Default value of option IntIfLoadedThenIntegrateUsedIntegrateQ, initially False";
$DefaultMostGeneralFirstQ::usage = "Default value of option MostGeneralFirstQ, initially True.";
$DefaultNLimitUse::usage = "Default value of option NLimitUse, initially  None.";
$DefaultParameterTreatment::usage =
 "Default value of option ParameterTreatment, initially the string \"Comprehensive\".";
$DefaultShowIntegralStepsQ ::usage = "Default value of option ShowIntegralStepsQ, initially False.";
$DefaultSubtaskTimeLimit::usage = "Default value of option SubtaskTimeLimit, initially 4.0 seconds.";
$DefaultTrueInPiecewiseMeans::usage =
 "Default value of option TrueInPiecewiseMeans, initially Automatic.";


(* ::Text:: *)
(*Messages :*)


InvalidOptionValue::usage = "";
UnknownOptionValue::usage = "";

MaybeConsequentialTimeOut::usage =
 "Warning about a timeout that is likely to cause a less continuous, less comprehensive, or less concise result. "; 

MaybeDiscontinuous::usage =
 "Warning about a parent-child clause pair in a Piecewise result of  Integral [\[Ellipsis]]  that might be discontinuous for optional parameter value  ParameterTreatment \[Rule] Continuous  or  ParameterTreatment \[Rule] ApartContinuous.";

OnlyGenericAndConditionalAvailableForMathematica10::usage =
"Warning that the value of optional parameter named ParameterTreatment and $DefaultParameterTreatment to  Conditional  because you need Mathematica version 11.0 or later for values Comprehensive, Continuous, and ApartContinuous.";

OnlyGenericAvailableForMathematicaVersionsBefore10::usage =
"Warning that the value of optional parameter named ParameterTreatment and $DefaultParameterTreatment to  Generic  because you need Mathematica version 10.0 or later for  Conditional  and version 11.0 or later for values Comprehensive, Continuous, and ApartContinuous.";

PerhapsInconsequentialTimeOut::usage =
 "Warning that a timeOut occurred that MIGHT cause a less concise result than otherwise.";


(* ::Section::Closed:: *)
(*Initialize global variables and messages ...*)


Begin ["`Private`"];


On [Assert];


$IntegralVersionNumber = 1.1;


$DefaultGenericIntegrationFunction = IntThenIntegrate;
$DefaultIntThenIntegrateUsedIntegrateQ = False;
$DefaultMostGeneralFirstQ = True;
$DefaultNLimitUse = None;
$DefaultParameterTreatment = Comprehensive;
$DefaultShowIntegralStepsQ = False;
$DefaultSubtaskTimeLimit = 4.0;
$DefaultTrueInPiecewiseMeans = Automatic;

$DeprecatedInResults = {
  Integrate, DirectedInfinity,
  Rubi`Int, Rubi`Unintegrable. Rubi`CannonIntegrate, Rubi`Subst, Rubi`Dist};

$printedDefaultGenericIntegrationFunctionMessage = False;

$MaxAllowedNKahanianEntropy10ForRational = 1.7;
$MaxAllowedNKahanianEntropy10ForIrrational = 4.0;

$NLimitChopAbsTolerance = 1.0*^-5;
$NLimitWorkingPrecision = 64;

$RationalizeTolerance = 1.0*^-5;
$RubiWasLoadedBeforeIntegralQ = MemberQ [$Packages, "Rubi`"];


Protect [Generic, Conditional, Comprehensive, Continuous, ApartContinuous];


InvalidOptionValue::SubtaskTimeLimit =
"Optional parameter SubtaskTimeLimit in function  Integral [\[Ellipsis]]  must be Infinity or a positive integer or floating-point number rather than  `1` .
\[Bullet] The implicit time unit is seconds.
\[Bullet] The default is the value of the reassignable global variable  $DefaultSubtaskTimeLimit, currently  `2`."; 

InvalidOptionValue::GenericIntegrationFunctionInt =
"The optional free Rubi antiderivative package from   https://rulebasedintegration.org/  must be loaded _before_ package  Integral.m  if you want to set the value of optional parameter  GenericIntegrationFunction  to  Int  or  IntThenIntegrate.
\[Bullet] I have just made the assignment  $DefaultGenericIntegrationFunction = Integrate  so that deleting your use of  GenericIntegrationFunction \[Rule] \[Ellipsis]  will now automatically use  Integrate [\[Ellipsis]].  
\[Bullet] If Rubi is installed on your computer and you want to exploit it, then use the Mathematica  Evaluation  menu to Quit the kernel, then Load Rubi, then re-load  Integral.m .";

UnknownOptionValue::GenericIntegrationFunction =
"`1` is an unknown option value for optional parameter  GenericIntegrationFunction  in function  Integral [...].
\[Bullet] Known alternatives are the symbols in  {IntThenIntegrate, Integrate, Int}.
\[Bullet] If the optional Rubi antiderivative package is not loaded before package  Integral.mAlternative, then  Int  is invalid and  IntThenIntegrate  is equivalent to  Integrate.
\[Bullet] The default is the value of the reassignable global variable  $DefaultGenericIntegrationFunction,  currently  `2`.";

UnknownOptionValue::ParameterTreatment =
"`1` is an unknown option value for optional parameter  ParameterTreatment  in function  Integral [...].
\[Bullet] Known successively more ambitious but slower alternatives are the Symbols in  `2`.
\[Bullet] The default is the value of the reassignable global variable  $DefaultParameterTreatment,  currently  `3`."; 

UnknownOptionValue::MostGeneralFirstQ =
"`1` is an unknown option value for optional parameter  MostGeneralFirstQ  in function  Integral [...].
\[Bullet] Known alternatives are the Symbols in  `2`.
\[Bullet] The default is the value of the reassignable global variable  $DefaultMostGeneralFirstQ,  currently  `3`."; 

UnknownOptionValue::TrueInPiecewiseMeans =
"`1` is an unknown option value for optional parameter  TrueInPiecewiseMeans  in function  Integral [...].
\[Bullet] Known alternatives are the Symbols in  `2`.
\[Bullet] The default is the value of the reassignable global variable  $DefaultTrueInPiecewiseMeans,  currently  `3`."; 

UnknownOptionValue::NLimitUse =
"`1` is an unknown option value for optional parameter  NLimitUse  in function  Integral [...].
\[Bullet] Known alternatives are the Symbols in  `2`.
\[Bullet] The default is the value of the reassignable global variable  $DefaultNLimitUse,  currently  `3`."; 

UnknownOptionValue::PrincipalValue =
"`1` is an unknown option value for optional parameter  PrincipalValue  in function  Integral [...].
\[Bullet] Known values are the symbols in  `2`.";

UnknownOptionValue::GenerateConditions =
"`1` is an unknown option value for optional parameter  GenerateConditions  in function  Integral [...].
\[Bullet] Known values are the symbols in  `2`.";

UnknownOptionValue::ShowIntegralStepsQ =
"`1` is an unknown option value for optional parameter  ShowIntegralStepsQ  in function  Integral [...].
\[Bullet] Known values are the symbols in  `2`.";


MaybeConsequentialTimeOut::subtask =
 "A timeOut occurred during `1`[\[Ellipsis], `2`, `3`] that is likely to cause a MaybeDiscontinuous message.
To request more subtask time than the value of the global variable $DefaultSubtaskTimeLimit, either use optional Integral [\[Ellipsis]] keyword argument  SubtaskTimeLimit \[Rule] \[Ellipsis], or else assign a larger value to $DefaultSubtaskTimeLimit."; 

MaybeDiscontinuous::mutualContinuity =
 "Case `1` might be parametrically discontinuous with case `2` in the Piecewise result of Integral [\[Ellipsis]].";
MaybeDiscontinuous::selfContinuity = "Subexpression `1` might be parametrically discontinuous.";
MaybeDiscontinuous::NLimitPropose =
"Used NLimit then the AskConstants Propose [\[Ellipsis]] function to determine `1` as a limit of `2`.";
MaybeDiscontinuous::NLimitRationalize =
"Used NLimit then the Rationalize [\[Ellipsis]] function to determine `1` as a limit of `2`.
Downloading then loading the AskConstants package from AskConstants.org enables identification of many IRRATIONAL floating-point constants too.";

OnlyGenericAndConditionalAvailableForMathematica10::sorry =
"Changed the value of optional parameter named ParameterTreatment and $DefaultParameterTreatment to Conditional because you need Mathematica version 11.0 or later for values Comprehensive, Continuous, and ApartContinuous.";
OnlyGenericAvailableForMathematicaVersionsBefore10::sorry =
"Changed the value of optional parameter named ParameterTreatment and $DefaultParameterTreatment to Generic because you need Mathematica version 10.0 or later for value Conditional, or version 11.0 or later for the additional alternatives Comprehensive, Continuous, and ApartContinuous";

PerhapsInconsequentialTimeOut::subtask =
 "A timeOut occurred during `1` that MIGHT cause a less concise result than otherwise.
To request more subtask time than the value of the global variable $DefaultSubtaskTimeLimit, either use optional Integral [\[Ellipsis]] keyword argument  SubtaskTimeLimit \[Rule] \[Ellipsis], or else assign a larger value to that variable.";

TrueInPiecewiseMeans::trueMeans = "True means `1` in the Piecewise result of Integral [\[Ellipsis]].";


(* ::Section::Closed:: *)
(*Function definitions ...*)


(* ::Subsubsection::Closed:: *)
(*0 ^ exponent_ := Piecewise [...] and fix bug in Simplify [...Piecewise..., Element...] ...*)


With [{test = 0 ^ \[FormalCapitalO]},
  If [Head [test] === Power && First [test] === 0,
       Unprotect [Power];
       0 ^ exponent_ := Piecewise [{{0, Re[exponent] > 0},
                                     {Indeterminate, Re[exponent] == 0}},
                                    ComplexInfinity];
       Protect [Power]]];

Unprotect [Element];
Element [Reals, Reals] = True;
Element [Integers, Reals] = True;
Element [Integers, Integers] = True;
Protect [Element];


(* ::Subsubsection::Closed:: *)
(*AllVariables [expression_] ...*)


ClearAll [AllVariables];

AllVariables [atom_? AtomQ] :=
  If [Head[atom] === Symbol && ! MemberQ [Attributes  [atom], Constant] &&
         ! MemberQ [{True, Indeterminate, Undefined, Less, Greater, LessEqual,
                      GreaterEqual, Equal, Unequal, Integers, Reals, Complexes,$TimedOut}, atom],
      {atom},
      {}];

AllVariables [f_[zeroOrMoreArguments___]] :=
   Apply [Union, Map [AllVariables, {zeroOrMoreArguments}]];


(* ::Subsubsection::Closed:: *)
(*ApartKahanian [antiderivative_, x_, subtaskTimeLimit_:$DefaultSubtaskTimeLimit] ...*)


ApartKahanian [antiderivative_, x_, subtaskTimeLimit_:$DefaultSubtaskTimeLimit] :=
  Module [{pieces = TermsFreeOfVarAndWithVar [antiderivative, x]},
    Apply [Plus, If [First [pieces] === 0 || Last [pieces] === 0,
                     pieces,
                     QuietTimeConstrainedSimp [pieces, subtaskTimeLimit]]]];


(* ::Subsubsection::Closed:: *)
(*BestSolution [solutions_] ...*)


BestSolution [solutions_] :=
  If [solutions === {},
      solutions,
      Module [{bestSolution = First [solutions],
                bestRightSide, result},
        bestRightSide = Last [bestSolution];
        Scan [If  [NumberQ [Last [#]] && ! NumberQ [bestRightSide] ||
                     NumericQ [Last [#]] && ! NumericQ [bestRightSide] ||
                     PolynomialOverGaussianRationalsQ [Last [#]] &&
                       ! PolynomialOverGaussianRationalsQ [bestRightSide] ||
                       RationalFunctionOverGaussianRationalsQ [Last [#]] &&
                         ! RationalFunctionOverGaussianRationalsQ [bestRightSide] ||
                       Entropy10 [Last [#]] < Entropy10 [bestRightSide],
                   bestSolution = #;
                   bestRightSide = Last [#]] &,  Rest [solutions]];
        {bestSolution}]];


(* ::Subsubsection::Closed:: *)
(*CandidateRulesForDiscontinuities [expression_,  varsOfInterest_,  excludedVars_, showStepsQ_:False] ...*)


CandidateRulesForDiscontinuities [expression_, varsOfInterest_, excludedVars_, showStepsQ_:False] :=
  Module [{allVariables = AllVariables [expression],
           domain = True,
           domainSubstitutions = {},
           singularitiesOfInterest, singularitySubstititons},
    If [$VersionNumber >= 11.0,
        singularitiesOfInterest = singularities [expression, allVariables,
          If [AllTrue [varsOfInterest, TrueQ [Refine [Element [#, Reals]]] &],
              {"DEFCUTS", "POLES", "ESSENTIAL", "BRANCHPOINTS", "IGNORE", "REAL"},
              {"POLES", "ESSENTIAL", "BRANCHPOINTS", "IGNORE"}]];
        singularitySubstititons = SingularitiesToReplacementRules [
                                    singularitiesOfInterest, varsOfInterest, excludedVars];
        If [showStepsQ,
            Print ["Singularites of interest = ", singularitiesOfInterest,
                   ", Derived substitutions = ", singularitySubstititons]];
        domainSubstitutions = Union [domainSubstitutions, singularitySubstititons]];
    SortBy [domainSubstitutions, Entropy10]];


(* ::Subsubsection::Closed:: *)
(*ComprehensiveAntiderivative [u_, x_, options:OptionsPattern[]] ...*)


(* ::Text:: *)
(*ComprehensiveAntiderivative:*)
(*1. Determines a generic antiderivative.*)
(*2. If that antiderivative is parametric, then ComprehensiveAntiderivative uses FunctionDomain [...] and singularities [...] to determine whether or not that antiderivative might have special cases.*)
(*3. If so and OptionValue [ParameterTreatment] === Comprehensive, then ComprehensiveAntiderivative invokes PiecewiseAntiderivative, which recursively identifies and assembles special cases by using independent antiderivatives for each case.*)
(*4. Otherwise, ComprehensiveAntiderivative:*)
(*    a) Determines the lowest complexity generic-case Kahanian that it can -- perhaps 0.*)
(*    b) Recursively invokes  ComprehensiveAntiderivativeAllCases to return a perhaps Piecewise result containing the generic and special cases, each with a special-case Kahanian some or all of which might be 0.*)
(*5. At recursion level 0, if the result of ComprehensiveAntiderivativeAllCases is piecewise and OptionValue [ParameterTreatment] === ApartContinuous, then ComprehensiveAntiderivative embellishes the returned result with an indication of the Boolean condition for the piecewise True clause.*)


ClearAll [ComprehensiveAntiderivative];

Options[ComprehensiveAntiderivative] = {
 Assumptions -> $Assumptions,
 GenerateConditions -> True,
 GenericIntegrationFunction :> $DefaultGenericIntegrationFunction,
 MostGeneralFirstQ :> $DefaultMostGeneralFirstQ,
 NLimitUse :> $DefaultNLimitUse,
 ParameterTreatment :> $DefaultParameterTreatment,
 PrincipalValue -> False,
 RecursionLevel -> 0, (* :Internal use only. *)  
 ShowIntegralStepsQ :> $DefaultShowIntegralStepsQ,
 SubtaskTimeLimit :> $DefaultSubtaskTimeLimit,
 TrueInPiecewiseMeans :> $DefaultTrueInPiecewiseMeans};

ComprehensiveAntiderivative [integrandOfx_, x_, options:OptionsPattern[]] :=
  Block [{genericIntegrationFunction = OptionValue [GenericIntegrationFunction],
          recursionLevel = OptionValue [RecursionLevel],
          showStepsQ = OptionValue [ShowIntegralStepsQ],
          subtaskTimeLimit  = OptionValue [SubtaskTimeLimit],
          trueInPiecewiseMeans = OptionValue [TrueInPiecewiseMeans],
          parameterTreatment = OptionValue [ParameterTreatment],
          rulesExceptRecursionLevel = FilterRules [{options}, Except [RecursionLevel]],
          realDomainQ = TrueQ [Refine [Element [x, Reals]]],
          allCases, allVariables, anchorCandidates, anchorSubstCandidates,
          antiderivative, bestAnchor, bestKahanian, complexityOfBestKahanian,
          complexityOfKahanian, integrandWithout0ToExponent, kahanian,
          kahanAntiderivative, parameters, parametricDiscontinuities,
          probablyDiscontinuousAnchorCandidates, solveAnchors, solveIntegrandAnchors, together},
    If [showStepsQ ,
        Print [""];
        If [recursionLevel === 0,
            Print ["Using ", $Version];
            If [$RubiWasLoadedBeforeIntegralQ,
                Print ["Rubi is loaded"]]];  
        Print ["Entered ComprehensiveAntiderivative with integrand = ",
               integrandOfx, ", and options = ", options]];
               
    If [OptionValue [ParameterTreatment] === Comprehensive,
        antiderivative = PiecewiseAntiderivative [integrandOfx, x,
          FilterRules [Prepend [rulesExceptRecursionLevel, RecursionLevel -> recursionLevel + 1],
                       Options[PiecewiseAntiderivative]]];
        Return [EmbellishComprehensiveExpression [antiderivative, x,
                  AllVariables [antiderivative], False,
                  OptionValue [MostGeneralFirstQ], subtaskTimeLimit, trueInPiecewiseMeans]]];
                  
    antiderivative = Together [RemoveSomeRemovableSingularities [Quiet [
                       If [genericIntegrationFunction === Integrate,
                           Integrate [integrandOfx, x, FilterRules[rulesExceptRecursionLevel,
                                      Options[Integrate]]],
                           genericIntegrationFunction [integrandOfx, x]]]],
                               Extension ->  Automatic, Trig -> True];                             
    (* First convert 0 ^ ... to Piecewise ?? : *)
    antiderivative = QuietTimeConstrainedSimp [antiderivative, subtaskTimeLimit];
    If [! FreeQAll [antiderivative, {Integrate, Rubi`Int}],
        Return [antiderivative]];
        
    allVariables = AllVariables [antiderivative];
    parameters = Complement [allVariables, {x}];
    If [Length [allVariables] <= 1 || UndeterminedQ [antiderivative] ||
          ! ProbablyHasDiscontinuitiesInVarsQ [antiderivative, parameters],
        If [showStepsQ,
            Print ["Returning ", antiderivative]];
        Return [antiderivative]];
        
    If [showStepsQ,
        Print ["Traditional antiderivative = ", antiderivative]];
    probablyDiscontinuousAnchorCandidates = CandidateRulesForDiscontinuities[
            antiderivative, {x}, {}(*parameters*), showStepsQ];
    parametricDiscontinuities = CandidateRulesForDiscontinuities [
            antiderivative, parameters, {x}, showStepsQ];
    solveAnchors = QuietTimeConstrainedSolveToReplacementRules [
                     antiderivative == 0, {x}, {}, subtaskTimeLimit];
    solveIntegrandAnchors = QuietTimeConstrainedSolveToReplacementRules [
                              integrandOfx == 0, {x}, {}, subtaskTimeLimit];
    anchorSubstCandidates = FunctionPeriod [antiderivative, x, Reals];
    If [showStepsQ,
        Print ["Antiderivative FunctionPeriod = ", anchorSubstCandidates]];
    anchorSubstCandidates = Join [{x -> 0, x -> 1, x -> -1, x -> 1/2, x -> -1/2},
                                  If [NumericQ  [anchorSubstCandidates],
                                      {x -> 1/2*anchorSubstCandidates, x -> 1/4*anchorSubstCandidates},
                                      {}],
                                  If [! FreeQ [antiderivative, Complex],
                                      {x -> I, x -> -I},
                                      {}]];
    anchorSubstCandidates = Select [DeleteDuplicates [Join [solveAnchors,
                              solveIntegrandAnchors, anchorSubstCandidates]],
                                    ! MemberQ [probablyDiscontinuousAnchorCandidates, #] &];
    anchorSubstCandidates = CullSubstitutions [anchorSubstCandidates, parametricDiscontinuities];
    anchorCandidates = DeleteDuplicates [Join [anchorSubstCandidates(*,
                         probablyDiscontinuousAnchorCandidates*)]];
    If [QuietTimeConstrainedLimit [integrandOfx, x -> -Infinity, -1, subtaskTimeLimit, False] === 0,
        AppendTo [anchorCandidates, x -> -Infinity]];
    If [QuietTimeConstrainedLimit [integrandOfx, x -> Infinity, 1, subtaskTimeLimit, False] === 0,
        AppendTo [anchorCandidates, x -> Infinity]];
    If [showStepsQ,
        Print ["solveAnchors = ", solveAnchors,
               ", solveIntegrandAnchors = ", solveIntegrandAnchors,
               ", probablyDiscontinuousAnchorCandidates = ", probablyDiscontinuousAnchorCandidates,
               ", anchorSubstCandidates = ", anchorSubstCandidates,
               ", anchorCandidates = ", anchorCandidates]];

    (* Compute generic bestAnchor and bestKahanian: *)
    bestAnchor = First [anchorCandidates];
    bestKahanian = If [MemberQ [anchorSubstCandidates, bestAnchor],
                       Quiet [antiderivative /. bestAnchor],
                       QuietTimeConstrainedLimit [antiderivative, bestAnchor, -1, subtaskTimeLimit, False]];
    bestKahanian = FailedOrTimeLimitedFullSimplifyToFinite [bestKahanian,realDomainQ,subtaskTimeLimit];
    Assert [Head [bestKahanian] =!= ConditionalExpression, ComprehensiveAntiderivative];
    If [showStepsQ,
        Print ["Initial anchor = ", bestAnchor, ", Initial Kahanian = ", bestKahanian]];
    If [bestKahanian =!= 0,
        complexityOfBestKahanian = If [bestKahanian === $Failed,
                                       Infinity,
                                       Entropy10 [bestKahanian]];
        Scan [(kahanian =
                If [MemberQ [anchorSubstCandidates, #],
                    Quiet [antiderivative /. #],
                    QuietTimeConstrainedLimit [antiderivative, #, -1, subtaskTimeLimit, False]];
               If [showStepsQ,
                   Print ["Finished internal Replace or Limit for anchor ", #, ", kahanian = ", kahanian]];
               kahanian = FailedOrTimeLimitedFullSimplifyToFinite [kahanian, realDomainQ, subtaskTimeLimit];
               If [kahanian === 0,
                   bestKahanian = kahanian;
                   Return []]; (* :Merely quit Scan *)
               Assert [Head [bestKahanian] =!= ConditionalExpression, ComprehensiveAntiderivative];
               If [kahanian =!= $Failed,
                   complexityOfKahanian = Entropy10 [kahanian];
                   If [complexityOfKahanian < complexityOfBestKahanian,
                       bestAnchor = #;
                       complexityOfBestKahanian = complexityOfKahanian;
                       bestKahanian = kahanian]]) &,
              Rest [anchorCandidates]]]; (* :end of If [bestKahanian =!= 0, ...] *)

    (* Compute generic kahanAntiderivative: *) 
    If [bestKahanian === $Failed,
        kahanAntiderivative = antiderivative,
    (* else bestKahanian =!= $Failed: *)
        (* Together might transformationally remove some removable singularities: *)
        kahanAntiderivative = antiderivative - bestKahanian;
        If [showStepsQ,
            Print ["kahanAntiderivative before Together = ", kahanAntiderivative,
            ", bestKahanian = ", bestKahanian], ", bestAnchor = ", bestAnchor];
        together = Together [kahanAntiderivative, Extension -> Automatic, Trig -> True];
        If [showStepsQ,
            Print ["kahanAntiderivative after Together = ", together]];
        If [Length [FunctionNames [together]] <= Length [FunctionNames [kahanAntiderivative]],
            (* Do not use FullSimplify here: *)
            kahanAntiderivative = QuietTimeConstrainedSimp [together, subtaskTimeLimit];
            If [showStepsQ,
                Print ["kahanAntiderivative after Simplify = ", kahanAntiderivative]]]];

    (* Combine generic kahanAntiderivative with special cases: *)
    anchorSubstCandidates = Select [anchorSubstCandidates, NumericQ [Last [#]] &];
    kahanAntiderivative = ComprehensiveAntiderivativeAllCases [
      integrandOfx, x, kahanAntiderivative, parameters, anchorSubstCandidates,
      FilterRules [Prepend [rulesExceptRecursionLevel, RecursionLevel -> recursionLevel + 1],
                   Options [ComprehensiveAntiderivativeAllCases]]];
    If [recursionLevel =!= 0,
        If [bestKahanian === $Failed,
            Message [MaybeDiscontinuous::selfContinuity, kahanAntiderivative]];
        kahanAntiderivative,  (* :?? *)
        If [Head [kahanAntiderivative] =!= Piecewise &&
              StrictEquivalentQ [antiderivative, kahanAntiderivative] &&
              Entropy10 [antiderivative] < Entropy10 [kahanAntiderivative],
            If [showStepsQ,
                Print ["Replaced ", kahanAntiderivative, " with lower Entropy10 equivalent ",
                       antiderivative]];
            kahanAntiderivative = antiderivative];
        EmbellishComprehensiveExpression [kahanAntiderivative, x,
           AllVariables [kahanAntiderivative], OptionValue [ParameterTreatment] === ApartContinuous, 
           OptionValue [MostGeneralFirstQ], subtaskTimeLimit, trueInPiecewiseMeans]]];


(* ::Subsubsection::Closed:: *)
(*ComprehensiveAntiderivativeAllCases [integrandOfx_, x_, antiderivative_,  parameters_, anchorSubstCandidates_, options:OptionsPattern[]] ...*)


(* ::Text:: *)
(*For each parameter value where the generic antiderivative might be inapplicable;*)
(*1. ComprehensiveAntiderivativeAllCases [...] makes that substitution.*)
(*2. If that substitution does not immediately reveal a discontinuity, then ComprehensiveAntiderivativeAllCases recursively applies itself to the substituted antiderivative in case one of its subcases causes a discontinuity.*)
(*3. If that substitution reveals a discontinuity, then ComprehensiveAntiderivativeAllCases invokes SpecialCaseComprehensiveAntiderivOrFailed [...] to compute the special case, which will be piecewise if it has special cases.*)
(*4. If there are one or more special cases, then ComprehensiveAntiderivativeAllCases creates a corresponding Piecewise result.*)


ClearAll [ComprehensiveAntiderivativeAllCases];

Options[ComprehensiveAntiderivativeAllCases] = {
 Assumptions -> $Assumptions,
 GenerateConditions -> True,
 GenericIntegrationFunction :> $DefaultGenericIntegrationFunction,
 MostGeneralFirstQ :> $DefaultMostGeneralFirstQ,
 NLimitUse :> $DefaultNLimitUse,
 ParameterTreatment :> $DefaultParameterTreatment,
 PrincipalValue -> False,
 RecursionLevel -> 0, (* :Internal use only. *)  
 ShowIntegralStepsQ :> $DefaultShowIntegralStepsQ,
 SubtaskTimeLimit :> $DefaultSubtaskTimeLimit,
 TrueInPiecewiseMeans :> $DefaultTrueInPiecewiseMeans};

ComprehensiveAntiderivativeAllCases [integrandOfx_, x_, antiderivative_,
    parameters_, anchorSubstCandidates_, options:OptionsPattern[]] :=
  Block [{(*$Assumptions = True,*)
          genericIntegrationFunction = OptionValue [GenericIntegrationFunction],
          recursionLevel = OptionValue [RecursionLevel],
          showStepsQ = OptionValue [ShowIntegralStepsQ],
          subtaskTimeLimit  = OptionValue [SubtaskTimeLimit],
          rulesExceptRecursionLevel = FilterRules [{options}, Except [RecursionLevel]],
          (* The following local var avoids a Mathematica \[LessEqual] 11.1.1.0 Piecewise bug: *)
          antideriv = antiderivative, 
          antiderivOfSubstIntegrand, bestAntideriv,exponent, limFromLeft, limFromRight,
          limitsOfKahanIntegral, nKahanian, numericAnchorCandidates,
          parameterSubstitutions,sowedNKahanianQ,
          specialCase, specialCases, substitutedAntiderivative, substitutedIntegrand},
    If [parameters === {},
        Return [antideriv]];
    If [showStepsQ,
        Print [{options}];
        Print ["Entered ComprehensiveAntiderivativeAllCases with antiderivative = ", antideriv,
               ", parameters = ", parameters,
               ", anchorSubstCandidates = ",anchorSubstCandidates,
               ", options = ", options]];
    parameterSubstitutions = CandidateRulesForDiscontinuities[antideriv, parameters, {x}, showStepsQ];
    If [showStepsQ,
        Print ["candidateParameterSubstitutions = ",  parameterSubstitutions]];
    If [parameterSubstitutions === {},
        Return [antideriv]];
    numericAnchorCandidates = Select [anchorSubstCandidates, First[#] === x && NumericQ [Last[#]] &];
    specialCases = Last [Reap [Scan [(
      substitutedIntegrand = Quiet [integrandOfx /. #];
      substitutedIntegrand = QuietTimeConstrainedSimp [substitutedIntegrand, subtaskTimeLimit];
      antiderivOfSubstIntegrand = Together [RemoveSomeRemovableSingularities [Quiet [
        If [genericIntegrationFunction === Integrate,
            Integrate [substitutedIntegrand, x,
              FilterRules[rulesExceptRecursionLevel, Options[Integrate]]],
            genericIntegrationFunction [substitutedIntegrand, x]]]],
                                            Extension ->  Automatic, Trig -> True];
      antiderivOfSubstIntegrand = QuietTimeConstrainedSimp [antiderivOfSubstIntegrand, subtaskTimeLimit];
      If [showStepsQ,
          Print ["Parameter substitution is ", #,
                 ", substitutedIntegrand = ", substitutedIntegrand,
                 ", antiderivOfSubstIntegrand = ", antiderivOfSubstIntegrand ]];
      Which [
  (**) Head [substitutedIntegrand] === Power && First [substitutedIntegrand] === 0,
        exponent = Last [substitutedIntegrand];
        Sow [{-Infinity, First[#] == Last[#] && Im [exponent] == 0 && Re [exponent] < 0}];
        Sow [{0, First[#] == Last[#] && (Im [exponent] =!= 0 || Re [exponent] >= 0)}],
  (**) DeprecatedInResultsQ [antiderivOfSubstIntegrand],
        Sow [{antiderivOfSubstIntegrand ,  SubstitutionToEquality [#]}],
  (**) substitutedAntiderivative = Quiet [antideriv /. #];
        bestAntideriv = substitutedAntiderivative =
          UndeterminedOrTimeLimitedFullSimplify [substitutedAntiderivative,subtaskTimeLimit]; 
        If [showStepsQ,
            Print [", substitutedAntiderivative = ", substitutedAntiderivative]];
        ! DeprecatedInResultsQ [substitutedAntiderivative] &&
           AllVariables [substitutedAntiderivative] === AllVariables [antiderivOfSubstIntegrand],
          (* :Not a special case, but some of its subcases might be. *)
          If [StrictEquivalentQ [antiderivOfSubstIntegrand, bestAntideriv] &&
                Entropy10 [antiderivOfSubstIntegrand] < Entropy10 [substitutedAntiderivative],
              (* antiderivOfSubstIntegrand is a simpler equivalent: *)
              bestAntideriv = antiderivOfSubstIntegrand];
          bestAntideriv = ComprehensiveAntiderivativeAllCases  [
            substitutedIntegrand, x, bestAntideriv,
            Complement [AllVariables [bestAntideriv], {x}], anchorSubstCandidates,
            FilterRules[Prepend [rulesExceptRecursionLevel, RecursionLevel -> recursionLevel + 1],
                        Options [ComprehensiveAntiderivativeAllCases]]];
            If [Head [bestAntideriv] === Piecewise, (* :There was a special SUBcase. *)
                 Sow [{bestAntideriv, SubstitutionToEquality [#]}]],
  (**) True, (* Rubi`Subst in antideriv \[Rule] deprecated. This case is immediately special: *)
         specialCase = SpecialCaseComprehensiveAntiderivOrFailed [
           integrandOfx, x, antideriv, #, anchorSubstCandidates,
           numericAnchorCandidates, parameterSubstitutions,
           FilterRules[Prepend [rulesExceptRecursionLevel, RecursionLevel -> recursionLevel + 1],
                       Options [SpecialCaseComprehensiveAntiderivOrFailed]]];
         If [specialCase =!= $Failed,
             Sow [{specialCase, SubstitutionToEquality [#]}]]]) &, parameterSubstitutions]]];
    If [specialCases =!= {},
        specialCases = First [specialCases]];
    If [showStepsQ,
        Print ["specialCases in ComprehensiveAntiderivativeAllCases = ", specialCases]];
    If [specialCases =!= {},
        $Assumptions = True;
        antideriv = Piecewise [specialCases, antideriv];
        If [showStepsQ,
            Print ["antiderivative = ", antideriv]];
        antideriv = SortedPiecewiseExpand [antideriv, x,
                                           OptionValue [MostGeneralFirstQ], subtaskTimeLimit]];
    If [showStepsQ,
        Print ["Antiderivative returned from ComprehensiveAntiderivativeAllCases = ", antideriv]];
    antideriv];


(* ::Subsubsection::Closed:: *)
(*CullSubstitutions  [substitutions_, singularities_] ...*)


CullSubstitutions  [substitutions_, singularities_]:=
  Select [substitutions, NumericQ [Last[#]] || FiniteForAllSubstitutions [Last[#], singularities] &];


(* ::Subsubsection::Closed:: *)
(*DeprecatedInResultsQ  [expression_] ...*)


(* Listed in order of decreasing desirability: *)
DeprecatedInResultsQ [expression_] := UndeterminedQ [expression] || ! FreeQAll [expression, $DeprecatedInResults];


(* ::Subsubsection::Closed:: *)
(*DiscardNegligibleImOrRe  [number_] ...*)


DiscardNegligibleImOrRe [number_] :=
  With [{re = Re [number],
         im = Im [number]},
    (* Exploit "fuzzy ==": *)     
    Which [re == re + im,   re,
           im == re + im,   im,
           True,            number]];


(* ::Subsubsection::Closed:: *)
(*EmbellishComprehensiveExpression [expression_, integrationVar_, allVars_,  apartKahanianQ_, mostGeneralFirstQ_, subtaskTimeLimit_,trueInPiecewiseMeans_] ...*)


EmbellishComprehensiveExpression [expression_, integrationVar_, allVars_,
       apartKahanianQ_, mostGeneralFirstQ_, subTaskTimeLimit_, trueInPiecewiseMeans_] :=
  Module [{result = expression,
           condition = $Assumptions,
           workingDomain},
    If [! TrueQ [condition],
        workingDomain = If [TrueQ [Refine [Element [integrationVar, Reals]]],
                            Reals,
                            Complexes];
        (*If [workingDomain === Reals && ! AllTrue [allVars, TrueQ [Refine [Element [#, Reals]]]&],
              Message [PartialReal::trueMeans, $Assumptions, allVars, integrationVar]];*)
        condition = condition && Which [
          FreeQ [$Assumptions, Element],
            True,
          Head [result] =!= Piecewise,
            FunctionDomain [result, allVars, workingDomain],
          True,
            Apply [Or, Append [
              Map [FunctionDomain [First[#], AllVariables [First[#]], workingDomain] && Last[#] &,
                   First [result]],
              FunctionDomain [Last [result], AllVariables [Last [result]], workingDomain] &&
                InPiecewiseTrueMeans [result]]]];
        result = QuietTimeConstrainedSimp [ConditionalExpression [result, condition], subtaskTimeLimit];
        If [Head [result] =!= ConditionalExpression,
            result = ConditionalExpression [result, condition]]];
    result = PossiblySeparateIndependentExpressions[result, integrationVar, apartKahanianQ, subtaskTimeLimit];
    result = PossibleReSortPiecewise [result, integrationVar, mostGeneralFirstQ, subtaskTimeLimit];
    PossibleTrueMeansMessage [result, trueInPiecewiseMeans];
    result]


(* ::Subsubsection::Closed:: *)
(*Entropy10 [expression_] ...*)


If [DownValues [Entropy10]=== {},  (* AskConstants is not already loaded *)

    $PrecisionOfMachineNon0 = $MachinePrecision - 2.4;

    $OperationEntropy10 = 0.5;
    $PureFunctionEntropy10Overhead = 2.0;

    $CommonFunctionList = {(*Power, *)Log, Sin, Cos, Tan, ArcTan, ArcSin, ArcCos};
    $CommonFunctionEntropy10 = 0.75;

    $LessCommonFunctionList = {Cot, Csc, Sec, ArcCot, ArcCsc, ArcSec,
      Sinh, Cosh, Tanh, Csch, Sech, Coth, ArcSinh, ArcCosh, ArcTanh, ArcCoth, ArcCsch, ArcSech,
      Erf, Erfc, Factorial, Gamma, Zeta, ProductLog, PolyLog, PolyGamma};
    $LessCommonFunctionEntropy10 = 1.0;

    $RareFunctionEntropy10 = 1.5;
    
    Unprotect [Entropy10];

    ClearAll [Entropy10];

    Entropy10 [0.0] := 0.1;

    Entropy10 [integer_Integer] :=
      Which [
      (**) Abs [integer] === 1,  0.2,
      (**) integer =!= 0,        RealExponent [integer],
      (**) True,                 0.1] + If [integer < 0,  0.25,  0.0];

    Entropy10 [rat_Rational] := 0.5 + Entropy10 [Numerator [rat] * Denominator [rat]];

    (* Should not occur for Propose[...], but included for completeness: *)
    Entropy10 [Complex [re_, im_]] := 1.5 + Entropy10 [re] + Entropy10 [im];

    (* Should not occur for Propose[...], but included for completeness: *)
    Entropy10 [$TimedOut] := 1000.0;

    (* Should not occur for Propose[...], but included for completeness: *)
    Entropy10 [Undefined] := 1000.0;

    (* Should not occur for Propose[...], but included for completeness: *)
    Entropy10 [Indeterminate] := 900.0;

    Entropy10 [ComplexInfinity] := 800000;

    (* Should not occur for Propose[...], but included for completeness: *)
    Entropy10 [directedInfinity_DirectedInfinity] := 700.0 + Entropy10 [First[directedInfinity]];

    (* Should not occur for Propose[...], but included for completeness: *)
    Entropy10 [real_Real] :=
      Which [
      (**) real === Underflow [],  Infinity,
      (**) real === Overflow [],   Infinity,
      (**) MachineNumberQ [real],  $PrecisionOfMachineNon0,
      (**) True,                    Precision [real]];

    Entropy10 [Pi] := 0.5;

    Entropy10 [E] := 0.625;

    Entropy10 [EulerGamma | Catalan | Glaisher | Khinchin] := 1.0;

    Entropy10 [\[FormalT] | \[FormalS]] = 0.0;

    Entropy10 [_? AtomQ] := 2.0;

    Entropy10 [Sqrt [base_]] := $CommonFunctionEntropy10 + Entropy10 [base];

    Entropy10 [Root [pureFunction_, ___]] :=
    Entropy10 [pureFunction] - $PureFunctionEntropy10Overhead;

    Entropy10 [functionalForm_] :=
      With [{head = Head [functionalForm]},
       Total [Map [Entropy10, Apply [List, functionalForm]]]
       + Max [1.0, Length [functionalForm] - 1] 
       * Which [
(**) head === Plus,            0.5,
              (**) head === Times,           0.5,
              (**) head === Power,           0.5,
              (**) MemberQ [$CommonFunctionList, head],         $CommonFunctionEntropy10,
              (**) MemberQ [$LessCommonFunctionList, head], $LessCommonFunctionEntropy10,
              (**) True,                                  $RareFunctionEntropy10]];

Protect [Entropy10]

];  (* end If [DownValues [Entropy10]=!= {}, ... *)


(* ::Subsubsection::Closed:: *)
(*ExtraBooleanSimplify [boolean_] ...*)


ClearAll [ExtraBooleanSimplify, ExtraBooleanSimplifyAux];

ExtraBooleanSimplify [u_] := FullSimplify [ExtraBooleanSimplifyAux [u]];

ExtraBooleanSimplifyAux [disjunction_Or /; ! FreeQAll [disjunction, {Re,Im}]]:=
  Not [Apply [And, Map [Not [ExtraBooleanSimplifyAux [#]] &, Apply [List, disjunction]]]];

ExtraBooleanSimplifyAux [conjunction_And /; ! FreeQAll [conjunction, {Re,Im}]]:=
  Not [Apply [Or, Map [Not [ExtraBooleanSimplifyAux [#]] &, Apply [List, conjunction]]]]

ExtraBooleanSimplifyAux [! u_] := ! ExtraBooleanSimplifyAux [u];

ExtraBooleanSimplifyAux [anythingElse_] := anythingElse;


(* ::Subsubsection::Closed:: *)
(*FailedOrTimeLimitedFullSimplifyToFinite [expression_, realDomainQ_, timeLimit_:$DefaultSubtaskTimeLimit] ...*)


FailedOrTimeLimitedFullSimplifyToFinite [
    expression_, realDomainQ_, subtaskTimeLimit_:$DefaultSubtaskTimeLimit]:=
  With [{simplified = QuietTimeConstrainedSimp [expression, subtaskTimeLimit]},
    If [UndeterminedQ [simplified] || ! FreeQ [simplified, DirectedInfinity ||
          Head [simplified] === ConditionalExpression (*??*) ||
          realDomainQ  && !FreeQ [simplified, Complex]], (*!TrueQ [Refine [Element [expression, Reals]]] *)
        $Failed,
        simplified]];






















(* ::Subsubsection::Closed:: *)
(*FiniteCases [piecewise_] ...*)


ClearAll [FiniteCases];

FiniteCases [piecewise_Piecewise]:=
  With [{finiteClauses = Select [Append [First [piecewise],
                                         {Last [piecewise], InPiecewiseTrueMeans[piecewise]}],
                                 FreeQAll [#, {DirectedInfinity, Indeterminate, Undefined}]&]},
    Switch  [Length [finiteClauses],
      0,  None,
      1,  finiteClauses[[1,1]],
      _,  Piecewise [Most [finiteClauses], First [Last [finiteClauses]]]]];
   
FiniteCases [anythingElse_] :=
  If [FreeQAll [anythingElse, {DirectedInfinity, Indeterminate, Undefined}],
      anythingElse,
      None];


(* ::Subsubsection::Closed:: *)
(*FiniteForAllSubstitutions  [expression_, singularities_] ...*)


FiniteForAllSubstitutions[expression_, singularities_]:=
  (* Should really consider every combination of one alternative for each parameter: *)
  Quiet [AllTrue  [singularities,
                   FreeQAll [expression /. #, {DirectedInfinity, Indeterminate}]&]];


(* ::Subsubsection::Closed:: *)
(*FreeQAll [expression_, items_] ...*)


FreeQAll[expression_,items_] :=
  Scan [If [!FreeQ [expression,#],
            Return [False]] &,
        items] === Null;


(* ::Subsubsection::Closed:: *)
(*FunctionNames [expression_] ...*)


ClearAll [FunctionNames];

FunctionNames [atom_? AtomQ] := {};

FunctionNames [(Times | Plus | Power) [args__]] := Apply [Union, Map [FunctionNames, {args}]];

FunctionNames [f_[args___]] :=  Union [{f}, Apply [Union, Map [FunctionNames, {args}]]];

FunctionNames [_] := {};


(* ::Subsubsection::Closed:: *)
(*InPiecewiseTrueMeans [piecewise_] ...*)


ClearAll [InPiecewiseTrueMeans];

InPiecewiseTrueMeans [piecewise_Piecewise] :=
  FullSimplify [! Apply [Or, Map [Last, First [piecewise]]]];

InPiecewiseTrueMeans [_] := "";


(* ::Subsubsection::Closed:: *)
(*Integral [integrand_, ...] ...*)


Unprotect [Integral];
ClearAll [Integral];

Options[Integral] = {
 Assumptions :> $Assumptions,
 GenerateConditions -> True,
 GenericIntegrationFunction :> $DefaultGenericIntegrationFunction,
 IntThenIntegrateUsedIntegrateQ :> $DefaultIntThenIntegrateUsedIntegrateQ,
 MostGeneralFirstQ :> $DefaultMostGeneralFirstQ,
 NLimitUse :> $DefaultNLimitUse,
 ParameterTreatment :> $DefaultParameterTreatment,
 PrincipalValue -> False,
 RecursionLevel -> 0, (* :Internal use only. *)  
 ShowIntegralStepsQ :> $DefaultShowIntegralStepsQ,
 SubtaskTimeLimit :> $DefaultSubtaskTimeLimit,
 TrueInPiecewiseMeans :> $DefaultTrueInPiecewiseMeans};

Integral [integrand_, Element [vars_, region_],  options:OptionsPattern[]] :=
  If [! ValidIntegralOptionsQ [options],
      $Failed,
      Integrate  [integrand,  Element [vars, region], FilterRules[{options}, Options[Integrate]]]];

Integral [integrand_, {var_, lowerLimit_, upperLimit_}, options:OptionsPattern[]] :=
  If [! ValidIntegralOptionsQ [options],
      $Failed,
      If [OptionValue [PrincipalValue] === True,
          Integrate [integrand, {var, lowerLimit, upperLimit},
                     FilterRules[{options}, Options[Integrate]]],
          Module [{antiderivative = Integral [integrand, var,
                                               FilterRules[{options}, Options[Integral]]],
                   valueAtLowerLim, valueAtUpperLim},
            If [FreeQ [antiderivative, Integrate] && TrueQ [Refine [
                  FunctionDomain [antiderivative, AllVariables [antiderivative], Complexes],
                                                                    lowerLimit < var < upperLimit]],
                valueAtUpperLim = Limit [antiderivative, var -> upperLimit, Direction -> 1];
                If [FreeQ [valueAtUpperLim, Limit],
                    valueAtLowerLim = Limit [antiderivative, var -> lowerLimit, Direction -> -1];
                    If [FreeQ [valueAtLowerLim, Limit],
                        Return [valueAtUpperLim - valueAtLowerLim]]]];
            Integrate [integrand, {var, lowerLimit, upperLimit},
                       FilterRules[{options}, Options[Integrate]]]]]];
                   
Integral [integrand_, var_, options:OptionsPattern[]] :=
  Block [{genericIntegrationFunction = OptionValue [GenericIntegrationFunction],
          parameterTreatment = OptionValue [ParameterTreatment],
          $neededIntegrateQ = False,
          antiderivative},
    If [! ValidIntegralOptionsQ [options],
        Return [$Failed]];
    If [$VersionNumber < 10.0 && parameterTreatment =!= Generic,
        parameterTreatment = $DefaultParameterTreatment = Generic;
        Message [OnlyGenericAvailableForMathematicaVersionsBefore10::sorry],
        If [$VersionNumber < 11.0 &&
               MemberQ [{Comprehensive, Continuous, ApartContinuous}, parameterTreatment],
            parameterTreatment = $DefaultParameterTreatment = Conditional;
            Message [OnlyGenericAndConditionalAvailableForMathematica10::sorry]]];
    
    antiderivative = Switch [parameterTreatment,
    Comprehensive,
       ComprehensiveAntiderivative [integrand, var,
         FilterRules [Join [{options}, {RecursionLevel -> 0}],
                      Options[ComprehensiveAntiderivative]]],
    Continuous,
       ComprehensiveAntiderivative [integrand, var,
         FilterRules [Join [{options}, {RecursionLevel -> 0}],
                      Options[ComprehensiveAntiderivative]]],
    ApartContinuous,
       ComprehensiveAntiderivative [integrand, var,
         FilterRules [Join [{options}, {RecursionLevel -> 0}],
                      Options[ComprehensiveAntiderivative]]],
    antiderivative =
      Switch [genericIntegrationFunction,
        Rubi`Int,
          With [{$Assumptions = OptionValue [Assumptions]},
                Rubi`Int [integrand, var]],
        IntThenIntegrate,
          With [{$Assumptions = OptionValue [Assumptions]},
                IntThenIntegrate [integrand, var(*,
                                          FilterRules [{options}, Options[Integrate]]*)]],
        _,
          Integrate [integrand, var(*, FilterRules[{options}, Options[Integrate]]*)]];
                            
    Conditional,   
       ConditionalExpression [antiderivative, FunctionDomain [
            antiderivative, AllVariables[antiderivative], Complexes]],
    _,  (* : Generic *)
       antiderivative];
    If [OptionValue [IntThenIntegrateUsedIntegrateQ] && $neededIntegrateQ,
        Print ["Function  IntThenIntegrate [\[Ellipsis]]  needed Integrate [\[Ellipsis]], which succeeded."]];
    antiderivative];
    
Protect [Integral];    


(* ::Subsubsection::Closed:: *)
(*IntThenIntegrate [integrand_, var_, options:OptionsPattern[]] ...*)


Unprotect [IntThenIntegrate];
ClearAll [IntThenIntegrate];

Options[IntThenIntegrate] =
 {Assumptions -> $Assumptions,
  GenerateConditions -> True,
  PrincipalValue -> False};

IntThenIntegrate [integrand_, var_, options:OptionsPattern[]] :=
  If [! $RubiWasLoadedBeforeIntegralQ,
      Integrate [integrand, var, options],
      Module [{result = Rubi`Int [integrand, var]},
        If [FreeQAll [ans, {Rubi`Int, Rubi`Subst}],
            result,
            result = ReplaceIntWithIntegrate [result,
                       FilterRules[{options}, Options[ReplaceIntWithIntegrate]]];
            If [FreeQAll [result, Integrate],
                $neededIntegrateQ = True];
            result]]];
            
Protect [IntThenIntegrate];


(* ::Subsubsection::Closed:: *)
(*Make Rubi do   Int [Piecewise [...] ...*)


Unprotect [Rubi`Int];
    Rubi`Int [piecewise_Piecewise, var_] := Integrate [piecewise, var];
Protect [Rubi`Int];


(* ::Subsubsection::Closed:: *)
(*OrderedClauses [clause1_, clause2_] ...*)


ClearAll [OrderedClauses];

OrderedClauses [Undefined, _] := True;

OrderedClauses [_, Undefined] := False;

OrderedClauses [Indeterminate, _] := True;

OrderedClauses [_, Indeterminate] := False;

OrderedClauses [ComplexInfinity, _] := True;

OrderedClauses [_, ComplexInfinity] := False;

OrderedClauses [infiniteMagnitude1_DirectedInfinity, infiniteMagnitude2_DirectedInfinity] := Entropy10 [First [infiniteMagnitude1]] < Entropy10 [First [infiniteMagnitude2]];

OrderedClauses [infiniteMagnitude_DirectedInfinity, _] := True;

OrderedClauses [_, infiniteMagnitude_DirectedInfinity] := False;

OrderedClauses [clause1_, clause2_] :=
  Module [{consequent1 = First [clause1],
           consequent2 = First [clause2],
           lenVars1, lenVars2},
    lenVars1 = Length [AllVariables [consequent1]];
    lenVars2 = Length [AllVariables [consequent2]];
    ! FreeQ [consequent1, DirectedInfinity] && FreeQ [consequent2, DirectedInfinity] ||
       lenVars1 < lenVars2 || lenVars1===lenVars2 && Entropy10 [clause1] < Entropy10[clause2]];


(* ::Subsubsection::Closed:: *)
(*PiecewiseAntiderivative [u_, x_, OptionsPattern[]] ...*)


ClearAll[PiecewiseAntiderivative];

Options[PiecewiseAntiderivative] = {
 Assumptions -> $Assumptions,
 GenerateConditions -> True, 
 GenericIntegrationFunction :> $DefaultGenericIntegrationFunction,
 MostGeneralFirstQ :> $DefaultMostGeneralFirstQ,
 PrincipalValue -> False,
 RecursionLevel -> 0, (* :Internal use only. *)
 ShowIntegralStepsQ :> $DefaultShowIntegralStepsQ,
 SubtaskTimeLimit :> $DefaultSubtaskTimeLimit};

PiecewiseAntiderivative [integrandOfx_, x_, options:OptionsPattern[]] :=
  Block [{allVariables = AllVariables [integrandOfx],
          genericIntegrationFunction = OptionValue [GenericIntegrationFunction],
          recursionLevel = OptionValue [RecursionLevel],
          showStepsQ = OptionValue [ShowIntegralStepsQ],
          subtaskTimeLimit  = OptionValue [SubtaskTimeLimit],
          rulesExceptRecursionLevel = FilterRules [{options}, Except [RecursionLevel]],
          (*$Assumptions = True,*)
          antiderivative, parameters, specialCases},
    If [showStepsQ,
        Print [""];
        Print ["Entering PiecewiseAntiderivative with Integrand = ", integrandOfx,
               ", options = ", options]];
    (* First convert for 0 ^ ... to Piecewise ?? : *)
    antiderivative = Together [RemoveSomeRemovableSingularities [Quiet [
                       If [genericIntegrationFunction === Integrate,
                           Integrate [integrandOfx, x,
                                      FilterRules[rulesExceptRecursionLevel, Options[Integrate]]], 
                           genericIntegrationFunction [integrandOfx, x]]]],
                               Extension ->  Automatic, Trig -> True];
    (*antiderivative = TimeConstrained [Simplify [
        antiderivative], subtaskTimeLimit, antiderivative]*);
    antiderivative = QuietTimeConstrainedSimp [antiderivative, subtaskTimeLimit];
    If [showStepsQ,
        Print ["antiderivative = ",antiderivative]];
    Which [
    (**) UndeterminedQ [antiderivative] (*|| !FreeQ [antiderivative, Integrate]*),
           antiderivative,
    (**) parameters = Complement [allVariables, {x}];
         parameters === {},
         (*If [recursionLevel === 1 || keepNonSpecialIntermediateCasesQ,*)
           antiderivative,
    (**) True,
            specialCases = PiecewiseAntiderivativeSpecialCases[integrandOfx, x, antiderivative,
              FilterRules[{Prepend [rulesExceptRecursionLevel, RecursionLevel -> recursionLevel + 1]},
                           Options [PiecewiseAntiderivativeSpecialCases]]];
            If [showStepsQ,
                Print ["specialCases in PiecewiseAntiderivative = ",specialCases]];
            If [specialCases === {},
                antiderivative,
                (*$Assumptions = True;*)
                (*If [recursionLevel > 1 && ! keepNonSpecialIntermediateCasesQ,*)
                antiderivative = Piecewise [specialCases,antiderivative];
                If [showStepsQ,
                    Print ["unsorted unexpanded all cases = ",specialCases]];
                antiderivative = SortedPiecewiseExpand [antiderivative, x, 
                                   OptionValue [MostGeneralFirstQ], subtaskTimeLimit];
                If [showStepsQ,
                    Print ["SortedPiecewiseExpand = ", antiderivative]];
                antiderivative]]];


(* ::Subsubsection::Closed:: *)
(*PiecewiseAntiderivativeSpecialCases [u_, x_, genericAntiderivative_, recursionLevel_:0, OptionsPattern[]] ...*)


ClearAll[PiecewiseAntiderivativeSpecialCases];

Options[PiecewiseAntiderivativeSpecialCases] = {
 Assumptions -> $Assumptions,
 GenerateConditions -> True, 
 GenericIntegrationFunction :> $DefaultGenericIntegrationFunction,
 MostGeneralFirstQ :> $DefaultMostGeneralFirstQ,
 PrincipalValue -> False,
 RecursionLevel -> 0, (* :Internal use only. *)
 ShowIntegralStepsQ :> $DefaultShowIntegralStepsQ,
 SubtaskTimeLimit :> $DefaultSubtaskTimeLimit};

PiecewiseAntiderivativeSpecialCases [
    integrandOfx_, x_, genericAntiderivative_, options:OptionsPattern[]]:=
  Block [{recursionLevel = OptionValue [RecursionLevel],
          showStepsQ = OptionValue [ShowIntegralStepsQ],
          rulesExceptRecursionLevel = FilterRules [{options}, Except [RecursionLevel]],
          (*$Assumptions = True,*)
          allVariables = AllVariables [integrandOfx],
          candidateParameterSubstitutions, exponent, integralOfSubstitutedIntegrand,
          parameters, specialCases, substitutedAntiderivative, substitutedIntegrand},
    If [showStepsQ,
        Print ["Entering PiecewiseAntiderivativeSpecialCases with Integrand = ", integrandOfx,
               ", genericAntiderivative = ", genericAntiderivative,
               ", options = ", options]];
    parameters = Complement [allVariables, {x}];
    If [parameters === {},
        {},
        candidateParameterSubstitutions = CandidateRulesForDiscontinuities[
                                            genericAntiderivative, parameters, {x}, showStepsQ];
        If [showStepsQ,
            Print ["candidateParameterSubstitutions = ", candidateParameterSubstitutions]];
        specialCases = Last [Reap [Scan [
          (substitutedAntiderivative = Quiet [genericAntiderivative /. #];
           substitutedIntegrand = Quiet [integrandOfx /. #];
           integralOfSubstitutedIntegrand = PiecewiseAntiderivative [substitutedIntegrand, x,
              FilterRules[{Prepend [rulesExceptRecursionLevel, RecursionLevel -> recursionLevel + 1]},
                          Options[PiecewiseAntiderivative]]];
           If [Head [integralOfSubstitutedIntegrand] === Piecewise || 
                 ! FreeQ [substitutedAntiderivative, DirectedInfinity] ||
                 UndeterminedQ [substitutedAntiderivative],
               Sow [{integralOfSubstitutedIntegrand, SubstitutionToEquality[#]}]];
           If [showStepsQ,
               Print ["integralOfSubstitutedIntegrand = ", integralOfSubstitutedIntegrand]])&,
                                         candidateParameterSubstitutions]]];
        If [specialCases === {},
            specialCases,
            First [specialCases]]]];     


(* ::Subsubsection::Closed:: *)
(*PlotSurfaceAnd2DCurve [surface,  curve_,  variableRange_,  parameterRange_,  resultRange_, parameterValue_] ...*)


Unprotect [PlotSurfaceAnd2DCurve];
ClearAll [PlotSurfaceAnd2DCurve];

Options [PlotSurfaceAnd2DCurve] = {Exclusions -> All};

PlotSurfaceAnd2DCurve [surface_, curve_, integrationVarRange_,
    parameterRange_, surfaceRange_, specialParameterRule_, options:OptionsPattern[]] :=
  Show [
    Plot3D [surface, integrationVarRange, parameterRange, PlotRange -> surfaceRange,
            AxesLabel ->{First[integrationVarRange], First [parameterRange]}, PlotPoints -> 100,
            MeshStyle -> Gray,
            Exclusions -> OptionValue [Exclusions]],
    ParametricPlot3D[{First[integrationVarRange], Last [specialParameterRule], curve},
                     integrationVarRange,
                     PlotRange -> surfaceRange,  PlotStyle -> Black,
                     Exclusions -> OptionValue [Exclusions]]];
                     
Protect [PlotSurfaceAnd2DCurve];


(* ::Subsubsection::Closed:: *)
(*PolynomialOverGaussianRationalsQ [expression_] ...*)


ClearAll [PolynomialOverGaussianRationalsQ];

PolynomialOverGaussianRationalsQ [number_? NumberQ] := True;

PolynomialOverGaussianRationalsQ [symbol_Symbol /;
   ! MemberQ [{Indeterminate, Undefined, E, Pi, EulerGamma, Khinchin, Glaisher},
              symbol]] := True;

PolynomialOverGaussianRationalsQ [
  (base_? PolynomialOverGaussianRationalsQ) ^ (posInt_Integer /; posInt > 0)] := True;

PolynomialOverGaussianRationalsQ [plus_Plus] :=
   AllTrue [plus, PolynomialOverGaussianRationalsQ];

PolynomialOverGaussianRationalsQ [times_Times] :=
   AllTrue [times, PolynomialOverGaussianRationalsQ];

PolynomialOverGaussianRationalsQ [_] := False;


(* ::Subsubsection::Closed:: *)
(*PossibleReSortPiecewise [expression_, integrationVar_, mostGeneralFirstQ_, subtaskTimeLimit_] ...*)


ClearAll [PossibleReSortPiecewise];

PossibleReSortPiecewise [ConditionalExpression [
    consequent_, condition_], integrationVar_, mostGeneralFirstQ_, subtaskTimeLimit_] :=
  ConditionalExpression [If [Head [consequent] === Piecewise,
                             PossibleReSortPiecewise [consequent, integrationVar, mostGeneralFirstQ, subtaskTimeLimit],
                               Last [SyntacticTermsFreeOfVarAndWithVar [consequent, integrationVar]]],
                         condition];

PossibleReSortPiecewise [
    expression_Piecewise, integrationVar_, mostGeneralFirstQ_, subtaskTimeLimit_] :=
  SortedPiecewiseExpand [expression, integrationVar, mostGeneralFirstQ, subtaskTimeLimit];

PossibleReSortPiecewise [expression_, __] := expression;


(* ::Subsubsection::Closed:: *)
(*PossiblySeparateIndependentExpressions [expression_, x_, apartKahanianQ_, subtaskTimeLimit_:$DefaultSubtaskTimeLimit] ...*)


ClearAll [PossiblySeparateIndependentExpressions];

PossiblySeparateIndependentExpressions [expression_, _, False, _] := expression;

PossiblySeparateIndependentExpressions [piecewise_Piecewise, x_, _,
    subtaskTimeLimit_:$DefaultSubtaskTimeLimit] :=
  Piecewise [Map [{PossiblySeparateIndependentExpressions [First[#],x, True, subtaskTimeLimit], Last[#]} &,
                  First[piecewise]],
             PossiblySeparateIndependentExpressions [Last [piecewise], x, True, subtaskTimeLimit]];

PossiblySeparateIndependentExpressions [ConditionalExpression [consequent_, condition_], x_, _,
    subtaskTimeLimit_:$DefaultSubtaskTimeLimit] :=
  ConditionalExpression [PossiblySeparateIndependentExpressions [consequent, x, _, subtaskTimeLimit],
                         condition];

PossiblySeparateIndependentExpressions [atom_? AtomQ, __] := atom;

PossiblySeparateIndependentExpressions [expression_, x_, _, subtaskTimeLimit_:$DefaultSubtaskTimeLimit] /;
   FreeQ [expression, Piecewise] := ApartKahanian [expression, x, subtaskTimeLimit];

(*PossiblySeparateIndependentExpressions [f_[args___], x_, _, subtaskTimeLimit_:$DefaultSubtaskTimeLimit] :=
  Apply [f, Map [PossiblySeparateIndependentExpressions [#, x, True, subtaskTimeLimit] &, {args}]];*)


(* ::Subsubsection::Closed:: *)
(*PossibleTrueMeansMessage [expression_, trueInPiecewiseMeans_] ...*)


ClearAll [PossibleTrueMeansMessage];

PossibleTrueMeansMessage [_, None] := Null;

PossibleTrueMeansMessage [conditionalExpression_ConditionalExpression, Automatic] :=
  PossibleTrueMeansMessage [First [conditionalExpression], Automatic];

PossibleTrueMeansMessage [conditionalExpression_ConditionalExpression, All] :=
  PossibleTrueMeansMessage [First [conditionalExpression], All];

PossibleTrueMeansMessage[piecewise_Piecewise, trueInPiecewiseMeans_]:=
  If [trueInPiecewiseMeans === All ||
          trueInPiecewiseMeans == Automatic && Length [First [piecewise]] > 1,
      Print ["In the Piecewise [...] result of Integral [...], True means ",
             InPiecewiseTrueMeans [piecewise], " :"]];

PossibleTrueMeansMessage [__] := Null;


(* ::Subsubsection::Closed:: *)
(*ProbablyHasDiscontinuitiesInVarsQ [expression_, vars_] ...*)


ProbablyHasDiscontinuitiesInVarsQ [expression_, vars_] :=
  Flatten [singularities [expression, vars,
             If  [AllTrue [vars, TrueQ [Refine [Element [#, Reals]]] &],
                  {"REALS", "POLES", "ESSENTIAL", "BRANCHPOINTS", "IGNORE"},
                  {"POLES", "ESSENTIAL", "BRANCHPOINTS", "IGNORE"}]]] =!= {};


(* ::Subsubsection::Closed:: *)
(*QuietLimitViaNLimitOrFailed [univariateLimitand_, varToValue_, direction_, scale_, showStepsQ_:False] ...*)


 (* Moderate magnitudes, Unlikely to be at discontinuities. Most likely to give at least one real result: *)

QuietLimitViaNLimitOrFailed[univariateLimitand_, varToValue_, direction_, scale_,
                            nLimitUse_, showStepsQ_:False] :=
  Module [{eulerSumLimit = Quiet [NLimit [univariateLimitand, varToValue,
                                           Direction -> direction, Scale -> scale,
                                           WorkingPrecision -> $NLimitWorkingPrecision,
                                           (*Terms \[Rule] 7, *)Method -> EulerSum]],
           sequenceLimit = Quiet [NLimit [univariateLimitand,varToValue,
                                           Direction -> direction, Scale -> scale,
                                           WorkingPrecision -> $NLimitWorkingPrecision,
                                           (*Terms \[Rule] 7,*)Method -> SequenceLimit]],
           floatLimits, limits},
    floatLimits = Sort [Map [DiscardNegligibleImOrRe,
                             Select [{eulerSumLimit, sequenceLimit}, NumberQ]], Precision];
    If  [True(*Q [showStepsQ]*),
         Print ["univariate limitand = ", univariateLimitand,
                ", varToValue = ", varToValue,
                ", scale \[Rule] ", scale,
                ", NLimits = ", floatLimits]];
    If [floatLimits === {},
        Return [$Failed]];            
    If [nLimitUse === InexactNumberQ,
        Return [First [floatLimits]]];
        
    floatLimits = Chop [SetPrecision [floatLimits, MachinePrecision], $NLimitChopAbsTolerance];  (* ?? *)
    If [True(*Q [showStepsQ]*),
        Print ["Chopped successful MachinePrecision Limits = ", floatLimits]];
    limits = If [MemberQ [$Packages, "Propose`"],
                 Map [First, Select [Map [Propose [#, MaxAcceptableCount -> 1] &, floatLimits],
                                     # =!= {} &]],
                 Select [Rationalize [floatLimits, $RationalizeTolerance],
                         ExactNumberQ [#] && Entropy10 [#] < $MaxAllowedNKahanianEntropy10ForRational &]];
    If [True(*Q [showStepsQ]*),
        Print ["limits after Propose or Rationalize = ", limits]];
    If [Length [limits] === 1 || Length[limits] === 2 && First [limits] === Last [limits],
        First [limits],  (* :Choose one with greatest Excellent Merit ?? *)
        If [nLimitUse === Automatic,
            Last [floatLimits],
            $Failed]]];






















(* ::Subsubsection::Closed:: *)
(*QuietTimeConstrainedSimp [expression_, timeLimit_:$DefaultSubtaskTimeLimit] ...*)


Unprotect  [QuietTimeConstrainedSimp];

QuietTimeConstrainedSimp[expression_,timeLimit_:$DefaultSubtaskTimeLimit] :=
  (* TO DO: Check to see if FullSimplify introduces removable singularities: *)
  Module [{simplified = Quiet [TimeConstrained [FullSimplify [expression], timeLimit]]},
    If [simplified === $Aborted,
        Message [PerhapsInconsequentialTimeOut::subtask, FullSimplify];
        Quiet [simplified = TimeConstrained [Simplify [expression], timeLimit];
               If [simplified === $Aborted,
                   Together [expression, Extension -> Automatic, Trig -> True],
                   simplified]],
        simplified]];

Protect [QuietTimeConstrainedSimp];


(* ::Subsubsection::Closed:: *)
(*QuietTimeConstrainedLimit [limitand_, substitution_, direction_, timeLimit_:$DefaultSubtaskTimeLimit] ...*)


QuietTimeConstrainedLimit [limitand_, substitution_, direction_:0,
    timeLimit_:$DefaultSubtaskTimeLimit, consequentialQ_:True] :=
  With [{limit = Quiet [TimeConstrained [Limit [
                    limitand, substitution, Direction -> direction], timeLimit, $TimedOut]]},
    If [limit === $TimedOut,
        If [consequentialQ,
            Message [MaybeConsequentialTimeOut::subtask, Limit,
                              substitution, Direction -> direction],
            Message [PerhapsInconsequentialTimeOut::subtask, Limit]]];
    (*If [Head [limit] === ConditionalExpression,
          Print [{limitand, substitution, direction, limit}]];*)
    limit];


(* ::Subsubsection::Closed:: *)
(*QuietTimeConstrainedSolveToReplacementRules [boolean_, varsOfInterest_, disallowedVars_, timeLimit_] ...*)


QuietTimeConstrainedSolveToReplacementRules [
    boolean_, varsOfInterest_, disallowedVars_, timeLimit_] :=
  Module [{solutions = Quiet [TimeConstrained [Solve [boolean, varsOfInterest], timeLimit]]},
    If [solutions === $Aborted || Head [solutions] === Solve,
        {},
        solutions = Map [If [Head [Last[#]]=== ConditionalExpression,
                             First [#] -> #[[2,1]],
                             #] &,  Flatten [solutions]];
        Select [Union [solutions(*/.{C[1] \[Rule]0,C[2]\[Rule]1}*)], (*:Too hard otherwise*)
                MemberQ [varsOfInterest, First [#]] && FreeQAll [#, disallowedVars] &]]];


(* ::Subsubsection::Closed:: *)
(*RationalFunctionOverGaussianRationalsQ [expression_] ...*)


RationalFunctionOverGaussianRationalsQ [expression_] :=
  With [{together = Together [expression]},
    PolynomialOverGaussianRationalsQ [Numerator [together]] &&
      PolynomialOverGaussianRationalsQ [Denominator [together]]];


(* ::Subsubsection::Closed:: *)
(*RemoveSomeRemovableSingularities [expression_] ...*)


ClearAll [RemoveSomeRemovableSingularities];

RemoveSomeRemovableSingularities [HeavisideTheta [u_]]:= UnitStep [u];

RemoveSomeRemovableSingularities [HeavisidePi [u_]] := UnitBox [u];

(*RemoveSomeRemovableSingularities [(Sin[u_]/u_)]:= Sinc [u];*)
RemoveSomeRemovableSingularities [(Sin[u_]/u_)*w_.]:= Sinc [u]*w;

RemoveSomeRemovableSingularities [Csc[u_]*u_] := Sinc [u] ^ -1;

(* Why Does this not work?? : *)
(*RemoveSomeRemovableSingularities [Sin [u_]^v_ / u_^v_]:= Sinc [u] ^ v;*)

RemoveSomeRemovableSingularities [f_[args__]] :=
  Apply [f, Map [RemoveSomeRemovableSingularities, {args}]];

RemoveSomeRemovableSingularities [anythingElse_] := anythingElse;


(* ::Subsubsection::Closed:: *)
(*ReplaceIntWithIntegrate [expression_] ...*)


ClearAll [ReplaceIntWithIntegrate];

Options[ReplaceIntWithIntegrate] =
  {Assumptions -> $Assumptions,
   GenerateConditions -> True,
   PrincipalValue -> False};

ReplaceIntWithIntegrate [Rubi`Int [expression_, var_], options:OptionsPattern[]] :=
  Integrate [expression, var, options];

ReplaceIntWithIntegrate [Defer [Rubi`Int] [expression_, var_], options:OptionsPattern[]] :=
  Integrate [expression, var, options]

ReplaceIntWithIntegrate [f_ [args__], options:OptionsPattern[]] :=
  Apply [f, Map [ReplaceIntWithIntegrate [#, options] &, {args}]];

ReplaceIntWithIntegrate [anythingElse_, ___] := anythingElse;


(* ::Subsubsection::Closed:: *)
(*ScaleForNLimit [var_, substitutions_] ...*)


ScaleForNLimit [var_, substitutions_] :=
  With [{minDifference = Min [Abs [Differences [Map [Last,Sort [Select [
           substitutions, First [#] === var && NumericQ [Last[#]] &]]]]]]},
    If [minDifference === Infinity,
        1,
        minDifference] / 8];


(* ::Subsubsection::Closed:: *)
(*singularities [expression, vars_, types_] ...*)


singularities = If [$VersionNumber < 12.0,
                    Simplify`FunctionSingularities,
                    FunctionProperties`Singularities];


(* ::Subsubsection::Closed:: *)
(*SingularitiesToReplacementRules [singularities_, varsOfInterest_, disallowedVars_] ...*)


SingularitiesToReplacementRules [singularities_, varsOfInterest_, disallowedVars_] :=
  Union [Flatten [Map [QuietTimeConstrainedSolveToReplacementRules [
                          Apply [And, #], varsOfInterest, disallowedVars, 0.2] &,
                       Flatten [singularities, 1]]]];


(* ::Subsubsection::Closed:: *)
(*SortedPiecewiseExpand [piecewise_Piecewise, x_, mostGeneralFirstQ_:True,  timeLimit_:$DefaultSubtaskTimeLimit] ...*)


ClearAll [SortedPiecewiseExpand];

SortedPiecewiseExpand [piecewise_Piecewise, x_,
    mostGeneralFirstQ_:True,  timeLimit_:$DefaultSubtaskTimeLimit]:=
  Block [{newPiecewise = piecewise,
          $Assumptions = True, (* :This & Block make PiecewiseExpand work correctly. *)
          flattenedClauses, explicitTrueClause, sortedClauses},
    If [!FreeQ [Apply [List, piecewise], Piecewise],
        newPiecewise =
          PiecewiseExpand [piecewise, Method -> {"OrderlessConditions" -> True,
                                                 "ValueSimplifier" -> Identity,
                                                 "ConditionSimplifier" -> FullSimplify}]];
    flattenedClauses= First [newPiecewise];
    flattenedClauses = Map [{First [#], ExtraBooleanSimplify [Last [#]]} &, flattenedClauses];
    explicitTrueClause = {Last [newPiecewise],
      ExtraBooleanSimplify [FullSimplify [!Apply [Or, Map [Last, flattenedClauses]]]]};
    AppendTo [flattenedClauses, explicitTrueClause];
    sortedClauses = Sort [flattenedClauses, OrderedClauses];
    (*Print ["sortedClauses = ", sortedClauses];*)
    If [mostGeneralFirstQ,
        sortedClauses = Reverse [sortedClauses]];
    If [sortedClauses === flattenedClauses,
        newPiecewise,
        newPiecewise =Piecewise [Most [sortedClauses], First [Last[sortedClauses]]];
        If [InPiecewiseTrueMeans [newPiecewise] === False,
Print [sortedClauses]; (* A debugging output *)
            Piecewise [Most [First [newPiecewise]], First [Last [First [newPiecewise]]]],
            newPiecewise]]];

SortedPiecewiseExpand [anythingElse_, ___]:= anythingElse;


(* ::Subsubsection::Closed:: *)
(*SpecialCaseComprehensiveAntiderivOrFailed [integrand_, integrateVar_, kahanAntiderivative_, parameterToValue_, anchorSubstCandidates_, numericAnchorCandidates_,  parameterSubstitutions_, options:OptionsPattern[]] ...*)


(* ::Text:: *)
(*For one parameter substitution that makes a generic antiderivative discontinuous, SpecialCaseComprehensiveAntiderivOrFailed [...] attempts to compute a special-case Kahanian as a limit of the Difference between the Kahan generic antiderivative and the antiderivative of the substituted integrand as the parameter value approaches its special value.*)
(**)
(*The reason for using a limit to compute the Kahanian rather than to compute the entire special-case Kahan antiderivative is that the Kahanian is independent of the integration variable, making the resulting limit faster and more likely to be computed.*)
(**)
(*Moreover, since the Kahanian is independent of the integration variable but the above difference generally depends on the integration variable, a sequence of non-special numeric values are substituted for the integration variable in this difference so that the limitand will have one less variable, making it much easier to compute the limit of the simplest of these substituted differences from the right and the left.*)
(**)
(*If either of these limits cannot be computed, and the substituted limitand is univariate, then NLimit is used to attempt computing a floating-point approximation. If it succeeds and the AskConstants program freely downloadable from AskConstants.org is loaded, then SpecialCaseComprehensiveAntiderivOrFailed invokes the AskConstants Propose function to attempt proposing a concise nearby symbolic constant such as Pi/7 + Sqrt[2]  for the exact limit.;  Otherwise the Rationalize function is used to attempt proposing a simple exact rational number for the exact limit.*)
(**)
(*If by any of these means a limit is determined, then it is used as the Kahanian for further recursion.*)
(*If only one 1-sided limit is computed or the two 1-sided limits differ, then the simpler one is used, achieving only one-sided continuity.*)


ClearAll [SpecialCaseComprehensiveAntiderivOrFailed];

Options[SpecialCaseComprehensiveAntiderivOrFailed] = {
 Assumptions -> $Assumptions,
 GenerateConditions -> True,
 GenericIntegrationFunction :> $DefaultGenericIntegrationFunction,
 MostGeneralFirstQ :> $DefaultMostGeneralFirstQ,
 NLimitUse :> $DefaultNLimitUse,
 ParameterTreatment :> $DefaultParameterTreatment,
 PrincipalValue -> False,
 RecursionLevel -> 0, (* :Internal use only. *)  
 ShowIntegralStepsQ :> $DefaultShowIntegralStepsQ,
 SubtaskTimeLimit :> $DefaultSubtaskTimeLimit,
 TrueInPiecewiseMeans :> $DefaultTrueInPiecewiseMeans};

SpecialCaseComprehensiveAntiderivOrFailed [integrand_, integrateVar_, kahanAntiderivative_,
    parameterToValue_, anchorSubstCandidates_, numericAnchorCandidates_,
    parameterSubstitutions_, options:OptionsPattern[]] :=
  Module [{(*$Assumptions = True,*)
           genericIntegrationFunction = OptionValue [GenericIntegrationFunction],
           nLimitUse = OptionValue [NLimitUse],
           recursionLevel = OptionValue [RecursionLevel],
           showStepsQ = OptionValue [ShowIntegralStepsQ],
           subtaskTimeLimit  = OptionValue [SubtaskTimeLimit],
           rulesExceptRecursionLevel = FilterRules [{options}, Except [RecursionLevel]],
           realDomainQ = TrueQ [Refine [Element [x, Reals]]],
           parameter = First [parameterToValue],
           NLimitFromLeftQ = False,
           NLimitFromRightQ = False,
           specialIntegrand = integrand /. parameterToValue,
           difference, finiteCases, limFromLeft, limFromRight, limit, limitands,
           parameters, sameLimitsQ, scale, simplestLimitand, specialAntiderivative,
           specialCase, minusSpecialKahanian, useLimFromLeft},
    specialIntegrand = QuietTimeConstrainedSimp [specialIntegrand, subtaskTimeLimit];
    (*If [UndeterminedQ [specialIntegrand],
          Return [Indeterminate]];*)
    specialAntiderivative = Together [RemoveSomeRemovableSingularities [Quiet [
                              If  [genericIntegrationFunction === Integrate,
                                   Integrate [specialIntegrand, integrateVar,
                                              FilterRules[{options}, Options[Integrate]]], 
                                   genericIntegrationFunction [specialIntegrand, integrateVar]]]],
                                      Extension -> Automatic, Trig -> True];
    specialAntiderivative = QuietTimeConstrainedSimp [specialAntiderivative, subtaskTimeLimit];
    If [showStepsQ,
        Print ["Entered SpecialCaseComprehensiveAntiderivOrFailed with integrand = ", integrand,
               ", kahanAntiderivative = ", kahanAntiderivative,
               ", parameterToValue = ", parameterToValue,
               ", anchorSubstCandidates = ", anchorSubstCandidates,
               ", numericAnchorCandidates = ", numericAnchorCandidates,
               ", parameterSubstitutions = ", parameterSubstitutions,
               ", options = ", options];
        Print ["specialIntegrand = ", specialIntegrand,
               ", specialAntiderivative = ", specialAntiderivative]];
    If [UndeterminedQ [specialAntiderivative],
        Return [specialAntiderivative]];

    (* First try to compute the consistent Kahanian for this special case: *)   
    difference = Together [kahanAntiderivative - specialAntiderivative,
                           Extension -> Automatic, Trig -> True];
    limitands = Map [QuietTimeConstrainedSimp [Quiet [difference /. # ], subtaskTimeLimit]&,
                     numericAnchorCandidates];
    If [showStepsQ,
        Print ["difference = ", difference, ", limitands = ", limitands]];
    limitands = SortBy [Select [limitands, !UndeterminedQ [#] && FreeQ [#, DirectedInfinity] &],
                        Entropy10];
    If [showStepsQ, Print ["limitands = ", limitands]];
    simplestLimitand = Cases [limitands, limitand_/; FreeQ [limitand, parameter], 1, 1];
    If [simplestLimitand =!= {}, (* We got lucky: simplestLimitand is the kahanian: *)
        limFromLeft = limFromRight = First [simplestLimitand];
        If [showStepsQ,
            Print ["Lucky: simplestLimitand is the special-case Kahanian = ", limFromLeft]],
    (* else simplestLimitand === {}: *)
        If [limitands === {},
            If [showStepsQ,
                Print ["Unlucky: All substitutions for integration var failed."]];
            (* We got unlucky. Try computing limit with integration var not substituted away: *)
            limFromLeft = QuietTimeConstrainedLimit [difference, parameterToValue, 1,
                                                     subtaskTimeLimit, True];
            limFromLeft = FailedOrTimeLimitedFullSimplifyToFinite [
                                              limFromLeft, realDomainQ, subtaskTimeLimit];
            Assert [Head [limFromLeft] =!= ConditionalExpression || "limDifferenceFromLeft",
                    SpecialCaseComprehensiveAntiderivOrFailed];
            limFromRight = QuietTimeConstrainedLimit [difference, parameterToValue, -1,
                                                      subtaskTimeLimit, True];
            limFromRight = FailedOrTimeLimitedFullSimplifyToFinite [
                                              limFromRight, realDomainQ, subtaskTimeLimit];
            Assert [Head [limFromLeft] =!= ConditionalExpression || "limDifferenceFromRight",
                    SpecialCaseComprehensiveAntiderivOrFailed],

        (* else limitands =!= {}. Medium luck, 1st try simplest limitand freeof integrate var: *)
            simplestLimitand = First [limitands];
            limFromLeft = QuietTimeConstrainedLimit [simplestLimitand, parameterToValue, 1,
                                                     subtaskTimeLimit, True];
            (*limFromLeft = UndeterminedOrTimeLimitedFullSimplify [limFromLeft, subtaskTimeLimit];*)
            limFromLeft = FailedOrTimeLimitedFullSimplifyToFinite [
                                              limFromLeft, realDomainQ, subtaskTimeLimit];
            If [! FreeQ [limFromLeft, DirectedInfinity],
                limFromLeft = $Failed,
                Assert [Head [limFromLeft] =!= ConditionalExpression || "limSubstFromLeft",
                        SpecialCaseComprehensiveAntiderivOrFailed];
                If [nLimitUse =!= None,
                    If [limFromLeft === Indeterminate, (* Try NLimit on any univariate substituted limitands: *)
                        scale = ScaleForNLimit[parameter,parameterSubstitutions];
                        limFromLeft = Scan [
                          If [AllVariables [#] === {parameter},
                              limit = QuietLimitViaNLimitOrFailed[#, parameterToValue, 1,
                                              scale, nLimitUse, showStepsQ];
                              If [limit =!= $Failed,
                                  NLimitFromLeftQ = True,
                                  Return [limit]]] &, limitands];  (* :Return only from Scan *)
                        If [limFromLeft === Null,
                            limFromLeft = $Failed]]]];
  
            (* Try the same sequence of limits for limFromRight: *)
            limFromRight = QuietTimeConstrainedLimit [
                             simplestLimitand, parameterToValue, -1, subtaskTimeLimit, True];          
            limFromRight = FailedOrTimeLimitedFullSimplifyToFinite [
                                               limFromRight, realDomainQ, subtaskTimeLimit];
            (*limFromRight = UndeterminedOrTimeLimitedFullSimplify [limFromRight, subtaskTimeLimit];*)

            Assert [Head [limFromRight] =!= ConditionalExpression || "limSubstFromRight",
                    SpecialCaseComprehensiveAntiderivOrFailed];
            If [! FreeQ [limFromRight, DirectedInfinity],
                limFromRight = $Failed,
                If [nLimitUse =!= None,
                    If [limFromLeft === Indeterminate,
(*Print ["scale = ", scale];*)
                        If [Head [scale]=== Symbol,
                            scale = ScaleForNLimit[parameter,parameterSubstitutions]];
                            limFromRight = Scan [
                              If [AllVariables [#] === {parameter},
                                  limit = QuietLimitViaNLimitOrFailed[#, parameterToValue, -1,
                                            scale, nLimitUse, showStepsQ];
                                  If [limit =!= $Failed,
                                      NLimitFromRightQ = True,
                                      Return [limit]]] &, limitands];  (* :Return only from Scan *)
                            If [limFromRight === Null,
                                limFromRight = $Failed]]]]]];

   (* Decide which limit to use as the Kahanian, if any: *)
   If [showStepsQ,
       Print ["limFromLeft = ", limFromLeft, ",  limFromRight = ", limFromRight]];
   If [limFromLeft =!= $Failed,
       If [limFromRight =!= $Failed,
           sameLimitsQ = StrictEquivalentQ [limFromLeft, limFromRight];
           useLimFromLeft = Entropy10 [limFromLeft] < Entropy10 [limFromRight],            
       (* else limFromRight === $Failed: *)
           sameLimitsQ = False;
           useLimFromLeft = True],
    (* else limFromLeft === $Failed: *)
       useLimFromLeft = False;
       If [limFromRight =!= $Failed,
           sameLimitsQ = False,
       (* else both limits === $Failed: *)
           sameLimitsQ = $Failed]];
    minusSpecialKahanian = If [useLimFromLeft,
                               limFromLeft,
                               limFromRight];
    If [showStepsQ,
        Print ["sameLimitsQ = ", sameLimitsQ, ", useLimFromLeft = ", useLimFromLeft,
               ", minusSpecialKahanian = ", minusSpecialKahanian]];

   (* Recursively compute raw special case, then optionally embellish it : *)
   If [UndeterminedQ [minusSpecialKahanian] || ! FreeQ [minusSpecialKahanian, DirectedInfinity],
        (* Continuity not guaranteed, so use antiderivative of substituted integrand: *)
       specialCase = ComprehensiveAntiderivative  [specialIntegrand, integrateVar,     
                       FilterRules[Prepend [rulesExceptRecursionLevel, RecursionLevel -> recursionLevel + 1],
                                   Options[ComprehensiveAntiderivative]]];
       finiteCases = FiniteCases [specialCase];
(*Print ["finiteCases = ", finiteCases];*)
       If [finiteCases =!= None && Head [finiteCases] =!= Piecewise,
           Message [MaybeDiscontinuous::mutualContinuity, finiteCases, kahanAntiderivative]];
       If [showStepsQ,
           Print ["specialCase = ", specialCase]];
       specialCase,
    (* else use best Kahanian for at least one-sided continuity: *)
       specialAntiderivative = Together [specialAntiderivative + minusSpecialKahanian,
                                         Extension -> Automatic, Trig -> True];
       parameters = Complement [AllVariables [specialAntiderivative], {integrateVar}];
       If [showStepsQ,
           Print ["minusSpecialKahanian = ",  minusSpecialKahanian,
                  ", specialAntiderivative = ", specialAntiderivative]];
       specialCase =
         If [parameters === {},
             specialAntiderivative,
             ComprehensiveAntiderivativeAllCases [specialIntegrand, integrateVar,
               Together [specialAntiderivative, Extension -> Automatic, Trig -> True],
               parameters, anchorSubstCandidates,
               FilterRules[Prepend [rulesExceptRecursionLevel, RecursionLevel -> recursionLevel + 1],
               Options [ComprehensiveAntiderivativeAllCases]]]];
       If [sameLimitsQ,
           If [(NLimitFromLeftQ || NLimitFromRightQ) &&
                  None =!= (finiteCases = FiniteCases [specialCase]),
               If [MemberQ [$Packages, "Propose`"],
                   Message [MaybeDiscontinuous::NLimitPropose, finiteCases, kahanAntiderivative],
                   Message [MaybeDiscontinuous::NLimitRationalize, finiteCases, kahanAntiderivative]]];
           If [showStepsQ,
               Print ["specialCase = ", specialCase]];
           specialCase,
       (* else ! sameLimitsQ: *)
           finiteCases = FiniteCases [specialCase];
(*Print ["finiteCases = ", finiteCases];*)
           If [finiteCases =!= None && Head [finiteCases] =!= Piecewise,
               Message [MaybeDiscontinuous::mutualContinuity, finiteCases, kahanAntiderivative]];
           If [showStepsQ,
               Print ["specialCase = ", specialCase]];
           specialCase]]];


(* ::Subsubsection::Closed:: *)
(*StrictEquivalentQ [expression1_, expression2_, timeLimit_:$DefaultSubtaskTimeLimit] ...*)


StrictEquivalentQ [expression1_, expression2_, subtaskTimeLimit_:$DefaultSubtaskTimeLimit] :=
    (* For example, StrictEquivalentQ [Infinity, Infinity] returns False *)
  ! (UndeterminedQ [expression1] || UndeterminedQ [expression2]) &&
  Quiet [Module [{difference = expression1 - expression2},
           difference = QuietTimeConstrainedSimp [difference, subtaskTimeLimit];
           TimeConstrained [PossibleZeroQ [difference], subtaskTimeLimit, False]]];


(* ::Subsubsection::Closed:: *)
(*SubstitutionsFromQuietTimeConstrainedSolve [expression_, vars_, timeLimit_, timeOutMessagesQ_] ...*)


SubstitutionsFromQuietTimeConstrainedSolve [relation_, vars_, timeLimit_] :=
  Union [Flatten [Map [
    SubstitutionsFromQuietTimeConstrainedSolveAux [relation, #, vars, timeLimit] &,
                       vars]]];

SubstitutionsFromQuietTimeConstrainedSolveAux [relation_, solutionVar_, allowedVars_, timeLimit_] :=
  With [{ans = Quiet [TimeConstrained [Solve [relation, solutionVar], timeLimit]]},
    Which [
(**) ans === $Aborted,
        Message [PerhapsInconsequentialTimeOut::subtask, Solve];
        {},  
(**) ! FreeQ [ans, Solve],
        {},
(**) True,
        Map [SubstitutionsFromSolveRule [#, allowedVars]&, Flatten [ans]]]];


(* ::Subsubsection::Closed:: *)
(*SubstitutionsFromSolveRule [expression_, solutionVars_] ...*)


ClearAll [SubstitutionsFromSolveRule];

(*SubstitutionsFromSolveRule [{}, solutionVars_] := {};

SubstitutionsFromSolveRule [{{}}, solutionVars_] := {};

SubstitutionsFromSolveRule [List [rules__], solutionVars_] :=
Apply [Union,Map [SubstitutionsFromSolveRule [#,solutionVars] &,
                                   rules]];*)

SubstitutionsFromSolveRule [
    solutionVar_ -> condExpr_ConditionalExpression,solutionVars_] /; MemberQ[solutionVars, solutionVar] :=
  SubstitutionsFromSolveRule [solutionVar -> First [condExpr], solutionVars];

SubstitutionsFromSolveRule [solutionVar_ -> numericExpr_? NumericQ, solutionVars_] /;
  MemberQ [solutionVars, solutionVar]:= {solutionVar -> numericExpr};

SubstitutionsFromSolveRule [solutionVar_ -> expr_/; !FreeQ[expr, C[1]], solutionVars_]/;
    MemberQ [solutionVars, solutionVar] :=
  SubstitutionsFromSolveRule [solutionVar -> Quiet [expr /. C[1] -> 0], solutionVars];

SubstitutionsFromSolveRule [solutionVar_ -> expr_/;!FreeQ[expr, C[2]], solutionVars_]/;
  MemberQ [solutionVars,solutionVar] :=
SubstitutionsFromSolveRule [solutionVar -> Quiet [expr /. C[2] -> 0], solutionVars];

SubstitutionsFromSolveRule [solutionVar_ -> 0^_,solutionVars_] /; 
  MemberQ [solutionVars,solutionVar] := {solutionVar -> 0, solutionVar -> 1};

SubstitutionsFromSolveRule [solutionVar_ -> expr_, solutionVars_] /;
    SubsetQ  [solutionVars, AllVariables [expr]] && FreeQ [expr, solutionVar] :=
  {solutionVar -> expr};

SubstitutionsFromSolveRule [_, _] := {};


(* ::Subsubsection::Closed:: *)
(*SubstitutionToEquality [substitution_] ...*)


SubstitutionToEquality [substitution_] :=
  With [{left = First [substitution],
          right = Last [substitution]},
    (FreeQ [right, C[1]] || Element [C[1],Integers]) &&
      (FreeQ [right, C[2]] || Element [C[2],Integers]) &&
      If [Denominator [right] =!= 1,
          Denominator [right] * left == Numerator [right],
          left == right]];


(* ::Subsubsection::Closed:: *)
(*SyntacticTermsFreeOfVarAndWithVar [expression_, var_] ...*)


SyntacticTermsFreeOfVarAndWithVar [expression_, var_] :=
  Which [
(**) Head [expression] === Plus,
          {Select [expression, FreeQ [#, var] &],  Select [expression, !FreeQ [#, var] &]},
(**) FreeQ [expression, var],
          {expression, 0},
(**) True,
          {0, expression}];       


(* ::Subsubsection::Closed:: *)
(*TermsFreeOfVarAndWithVar [termsAlreadyFreeOfVar_, expression_, var_] ...*)


TermsFreeOfVarAndWithVar [expression_,var_] :=
  (* I make an effort here to avoid morphing the two results more than necessary. *)
 Module [{termsFreeOfVar, termsWithVar, moreTermsFreeOfVar},
   {termsFreeOfVar, termsWithVar} = SyntacticTermsFreeOfVarAndWithVar [expression, var];

    (* Avoiding because could introduce removable singularities: *)
   (*{moreTermsFreeOfVar, termsWithVar} = 
         SyntacticTermsFreeOfVarAndWithVar [Apart [termsWithVar], var];
    termsFreeOfVar += moreTermsFreeOfVar;*)
   {moreTermsFreeOfVar, termsWithVar} = 
         SyntacticTermsFreeOfVarAndWithVar [ExpandAll [termsWithVar], var];
    termsFreeOfVar += moreTermsFreeOfVar;

   {moreTermsFreeOfVar, termsWithVar} = 
        SyntacticTermsFreeOfVarAndWithVar [FunctionExpand [termsWithVar], var];
    termsFreeOfVar += moreTermsFreeOfVar;

   {moreTermsFreeOfVar, termsWithVar} =
        SyntacticTermsFreeOfVarAndWithVar [ExpandAll [termsWithVar], var];
   Simplify [{termsFreeOfVar + moreTermsFreeOfVar, termsWithVar}]];


(* ::Subsubsection::Closed:: *)
(*UndeterminedOrTimeLimitedFullSimplify [expression_, timeLimit_:$DefaultSubtaskTimeLimit] ...*)


UndeterminedOrTimeLimitedFullSimplify [expression_, subtaskTimeLimit_:$DefaultSubtaskTimeLimit]:=
  If [UndeterminedQ [expression],
      Indeterminate,
      QuietTimeConstrainedSimp [expression, subtaskTimeLimit]];


(* ::Subsubsection::Closed:: *)
(*UndeterminedQ [expression_] ...*)


UndeterminedQ [expression_] :=
 ! FreeQAll [expression, {Undefined,Indeterminate,(*DirectedInfinity,*) Interval,
               $TimedOut, $Aborted, $Failed, Limit, Solve(*, Integrate, Rubi`Int, Series, Sum*)}];


(* ::Subsubsection::Closed:: *)
(*ValidIntegralOptionsQ [options:OptionsPattern[ ]] ...*)


$ValidNLimitUseOptions = {None, InexactNumberQ, ExactNumberQ, Automatic};
$ValidParameterTreatmentOptions =
  {Generic, Conditional, Comprehensive, Continuous, ApartContinuous};
$ValidMostGeneralFirstQOptions = $ValidPrincipalValueOptions =
  $ValidGenerateConditions = $ValidShowIntegralStepsQOptions = {True, False};
$ValidTrueInPiecewiseMeansOptions = {None, Automatic, All};


Options [ValidIntegralOptionsQ] = {
 ParameterTreatment :> $DefaultParameterTreatment,
 Assumptions :> $Assumptions,
 GenerateConditions -> True,
 GenericIntegrationFunction :> $DefaultGenericIntegrationFunction,
 MostGeneralFirstQ :> $DefaultMostGeneralFirstQ,
 NLimitUse :> $DefaultNLimitUse, 
 PrincipalValue -> False,
 RecursionLevel -> 0, (* :Internal use only. *)  
 ShowIntegralStepsQ :> $DefaultShowIntegralStepsQ,
 SubtaskTimeLimit :> $DefaultSubtaskTimeLimit,
 TrueInPiecewiseMeans :> $DefaultTrueInPiecewiseMeans};


ValidIntegralOptionsQ [options:OptionsPattern[]] :=
  Module [{allValidQ = True,
           optionValue = OptionValue [GenericIntegrationFunction],
           validOptionQ},
    If [! $printedDefaultGenericIntegrationFunctionMessage,
        Print ["Initialization report for package  Integral.m,  version ", $IntegralVersionNumber];      
        If [$RubiWasLoadedBeforeIntegralQ,
            $DefaultGenericIntegrationFunction = IntThenIntegrate;             
            Print ["\[Bullet] The optional Rubi antiderivative package was loaded before Integral.
\[Bullet] Therefore the default value of the optional parameter  GenericIntegrationFunction  was set to  IntThenIntegrate."],
            $DefaultGenericIntegrationFunction = Integrate;
            Print ["\[Bullet] The optional Rubi antiderivative package was not loaded before Integral.
\[Bullet] Therefore the default value of the optional parameter  GenericIntegrationFunction  was set to  Integrate."]]; 
        $printedDefaultGenericIntegrationFunctionMessage = True];

    If [! $RubiWasLoadedBeforeIntegralQ && ToString [optionValue] === "Int",
        $DefaultGenericIntegrationFunction = Integrate;
        Message [InvalidOptionValue::GenericIntegrationFunctionInt, optionValue];
        allValidQ = False,
    (* else: *)
        If [optionValue =!= Integrate && optionValue =!= IntThenIntegrate,
            Message [UnknownOptionValue::GenericIntegrationFunction, optionValue,
                     $DefaultGenericIntegrationFunction];
            allValidQ = False]];
                           
    optionValue = OptionValue [ParameterTreatment];
    If [! MemberQ [$ValidParameterTreatmentOptions, optionValue],
        Message [UnknownOptionValue::ParameterTreatment, optionValue,
                 $ValidParameterTreatmentOptions, $DefaultParameterTreatment];
        allValidQ = False];
        
    optionValue = OptionValue [SubtaskTimeLimit];    
    validOptionQ = (IntegerQ [optionValue] || Head[optionValue] === Real) && optionValue > 0.0 ||
                     optionValue === Infinity;
    If [! validOptionQ,
        Message [InvalidOptionValue::SubtaskTimeLimit, optionValue, $DefaultSubtaskTimeLimit];
        allValidQ = False];
        
    optionValue = OptionValue [MostGeneralFirstQ];    
    validOptionQ = MemberQ [$ValidMostGeneralFirstQOptions, optionValue];
    If [! validOptionQ,
        Message [UnknownOptionValue::MostGeneralFirstQ, optionValue,
                 $ValidMostGeneralFirstQOptions, $DefaultMostGeneralFirstQ];
        allValidQ = False];

    optionValue = OptionValue [TrueInPiecewiseMeans];    
    validOptionQ = MemberQ [$ValidTrueInPiecewiseMeansOptions, optionValue];
    If [! validOptionQ,
        Message [UnknownOptionValue::TrueInPiecewiseMeans, optionValue,
                 $ValidTrueInPiecewiseMeansOptions, $DefaultTrueInPiecewiseMeans];
        allValidQ = False];

    optionValue = OptionValue [NLimitUse];           
    validOptionQ = MemberQ [$ValidNLimitUseOptions, optionValue];
    If [! validOptionQ,
        Message [UnknownOptionValue::NLimitUse, optionValue,
                 $ValidNLimitUseOptions, $DefaultNLimitUse];
        allValidQ = False];
               
    optionValue = OptionValue [PrincipalValue];    
    validOptionQ = MemberQ [$ValidPrincipalValueOptions, optionValue];
    If [! validOptionQ,
        Message [UnknownOptionValue::PrincipalValue, optionValue,
                 $ValidPrincipalValueOptions];
        allValidQ = False];

    optionValue = OptionValue [GenerateConditions];    
    validOptionQ = MemberQ [$ValidGenerateConditions, optionValue];
    If [! validOptionQ,
        Message [UnknownOptionValue::GenerateConditions, optionValue,
                 $ValidGenerateConditions];
        allValidQ = False];

    (* Not validating Assumptions. *)
    
    optionValue = OptionValue [ShowIntegralStepsQ];    
    validOptionQ = MemberQ [$ValidShowIntegralStepsQOptions, optionValue];
    If [! validOptionQ,
        Message [UnknownOptionValue::ShowIntegralStepsQ, optionValue,
                 $ValidShowIntegralStepsQOptions];
        allValidQ = False];
    allValidQ];



(* ::Section::Closed:: *)
(*End private context then End the package ...*)


End [];
EndPackage [];
