Integration and Summation
of Piecewise Functions

Discussion

PiecewiseIntegrate

PiecewiseIntegrate[
        f, {x, xmin, xmax}, {y, ymin, ymax}, …] gives the definite integral of function f. It is intended for integrating piecewise continuous functions, and also generalized functions. It handles integrands and integration bounds involving the following expressions:

◆ UnitStep,  Sign,  Abs,  Min,  Max
◆ Floor,  Ceiling,  Round,  IntegerPart,  
      FractionalPart,  Quotient,  Mod
◆ DiracDelta and its derivatives,  
      DiscreteDelta,  KroneckerDelta
◆ If,  Which,  Element,  NotElement
◆ Piecewise,  Boole,  Clip

The basic idea is to decompose the integration region into components where the piecewise function can be simplified to one of its 'branches':

PiecewiseIntegrate[
      f[x, y, z], {x, 0, 1}, {y, 0, 1}, {z, 0, Max[x, y]}]

The arguments of the piecewise functions can also be nonlinear (and non-algebraic) and contain symbolic parameters:

PiecewiseIntegrate[
      UnitStep[-(x^2) + a x + 
          a], {x, -∞, ∞}]

PiecewiseIntegrate[Floor[a Log[x] + b], {x, 0, 1}, 
      Assumptions -> 
        a ≥ 0 ∧ 
          0 ≤ b ≤ 1]

Similar methods can be used to compute integrals involving DiracDelta:

PiecewiseIntegrate[
      DiracDelta[
        r - (Sqrt[(x^2) + (y^2)])], {x, -∞, ∞}, {y, -∞, ∞}]

When integrating expressions which contain the delta function δ(x), there is a question as to how (∫_( 0)^∞) δ(
            x) dx should be interpreted. The convention used by PiecewiseIntegrate is that boundary points are always treated the same way as interior points, that is, \(TraditionalForm\`∫\_\(\(\ \)\(0\)\)\%∞\( δ(
              x)\) dx = \(∫\_\(-∞\)\%\(\(\ \)\(0\)\)\(δ(x)\) dx = 1\)\). When adding integrals with intersecting boundaries you can always explicitly subtract the value of the integral over the boundary, because \(TraditionalForm\`∫\_\(\(\ \)\(0\)\)\%\(\(\ \)\(0\)\)\(δ(x)\) dx\) is also assumed by PiecewiseIntegrate to be equal to 1.

In most cases you should try to keep the integrand as simple as possible, e.g. integrating sums of piecewise functions term by term whenever it can be done. Otherwise the complexity of the decomposition increases very fast, particularly for Floor and other arithmetic piecewise functions. However, for complex integrands Integrate may be the most time-consuming part of the routine rather than Reduce, in which case it is advantageous to perform the decomposition of the whole expression first in order to minimize the number of calls to Integrate.

PiecewiseIntegrate doesn't support functions that are defined using conditions, e.g. as f[x_] := 1 /; 0 ≤ x ≤ 1. PiecewiseIntegrate needs to analyze the structure of the integrand, and even if it is possible to search through DownValues/UpValues/SubValues of all the symbols in the integrand to extract the conditions, it doesn't seem to be worth the effort. The above definition of f[x] can be easily rewritten as Which[0 ≤ x ≤ 1, 
        1, …].

PiecewiseSum

PiecewiseSum[
        f,  {i,  imin,  imax}, {j,  jmin,  jmax}, …] evaluates the sum of f. It uses the standard Mathematica iteration specification. A sum from -∞ to a with steps d > 0 is understood as the sum over the points a,  a - d,  a - 2  d,  ... , and similarly for other cases involving ±∞:

PiecewiseSum[(1/(k^2)), {k, ∞, n, -2}]

PiecewiseSum doesn't emulate the behaviour of Sum exactly:

PiecewiseSum[k, {k, n}]

This is different from the result given by Sum but consistent with the result obtained by the numerical evaluation of the sum for any value of the parameter n, including non-integer or negative values. To avoid getting unnecessary complicated results, you can explicitly specify that n is integer using the option Assumptions (see the Options section below):

PiecewiseSum[j, {i, 2  k + 1}, {j, (i/2)}, 
      Assumptions -> k ∈ Integers]

Again, this is in agreement with the numerical evaluation of the sum, because PiecewiseSum doesn't assume i/2 to be integer. PiecewiseSum can handle the same set of piecewise functions as PiecewiseIntegrate:

PiecewiseSum[Floor[2  Sin[((π k)/6)]] (1/(k^3)), {k, ∞}] // 
      FullSimplify

Naturally DiracDelta is more commonly used with PiecewiseIntegrate, and DiscreteDelta,  KroneckerDelta,  Element,  
      NotElement are more useful for PiecewiseSum:

PiecewiseSum[DiscreteDelta[k - a], {k, n}]

PiecewiseSum evaluates the sum term by term without performing the decomposition whenever possible. However, in some cases the best strategy is to evaluate the sum symbolically first and then substitute specific values for the parameters. A multiple sum such as the one below would be impossible to compute directly for a large n, while the decomposition effectively reduces the dimension of the sum due to the presence of DiscreteDelta:

PiecewiseSum[
      DiscreteDelta[i + j - k], {i, -n, n}, {j, -n, n}, {k, -n, 
        n}, 
Assumptions -> 
        n ∈ Integers]

Options

PiecewiseIntegrate and PiecewiseSum take the following options:

({{name, default, description}, {Assumptions, :>$Assumptions, assumptions about the parameters}, {ImplicitConditions, ->False, \(whether\ to\ generate\ conditions\ in\ the\ form\ of\ explicit\ ranges\)}, {Parameters, ->Automatic, \(the\ independent\ parameters\)}})

The meaning of the Assumptions option is the same as for Integrate. Notice the difference between the function Assuming and using Assumptions: Assuming adds new assumptions to $Assumptions while setting Assumptions overrides the current $Assumptions value.

Assuming[
      a > 0, 
Assuming[
        a < 1, 
PiecewiseIntegrate[
          UnitStep[a - x], {x, 0, 1}]]]
Assuming[
      a > 0, 
PiecewiseIntegrate[
        UnitStep[a - x], {x, 0, 1}, 
Assumptions -> 
          a < 1]]

ImplicitConditions determines whether the generated conditions may contain unsolved equations or inequalites describing the domain of the parameter values implicitly, e.g. n == Floor[n] which may often appear in the output of PiecewiseSum. This is particularly useful for complex sums:

PiecewiseSum[
      Max[i + j + k, a], {i, n}, {j, n}, {k, 
        n}, 
ImplicitConditions -> 
        True, 
Assumptions -> 
        5 < a < n ∧ n ∈ Integers]

The option Parameters can be used to explicitly specify the independent parameters. Its value can be a single parameter, a list of parameters or Automatic, the latter being equivalent to an empty list. In the following example PiecewiseIntegrate by default chooses UnitStep[Sin[a]] as the free parameter:

PiecewiseIntegrate[
      Abs[x], {x, 0, UnitStep[Sin[a]]}]
PiecewiseIntegrate[Abs[x], {x, 0, UnitStep[Sin[a]]}, 
      Parameters -> Sin[a]]
PiecewiseIntegrate[Abs[x], {x, 0, UnitStep[Sin[a]]}, 
      Parameters -> a]

NPiecewiseIntegrate/NPiecewiseSum

NPiecewiseIntegrate and NPiecewiseSum use a hybrid approach: first they decompose the integration/summation region into subregions and then compute numerical approximations to the integral or the sum over each subregion. This can significantly speed up the computation of integrals, at the cost of performing Reduce first. Also this method avoids the roundoff errors which may affect the numerical evaluation of discontinuous piecewise functions (UnitStep,  Floor and others) and make purely numerical methods unsuitable for such sums.

PiecewiseSum and NPiecewiseSum do not have the HoldAll attribute, because their primary purpose is to transform the summand to several simpler expressions, and in order to do that the summand needs to be evaluated and analyzed. For the same reason all Piecewise*/NPiecewise* functions treat If and Which (and Piecewise in Mathematica 5.1 and higher) as non-holding functions. Generally speaking, it is hard to expect to have complete control over the evaluation of the expression when performing integration or symbolic summation.

NPiecewiseIntegrate has the following options:

({{name, default, description}, {NIntegrateOptions, :>Options[
                NIntegrate], the options to pass to NIntegrate}, {NSumOptions, :>Options[
                NSum], the options to pass to NSum}})

NPiecewiseSum has the option NSumOptions only. NSumOptions is needed for NPiecewiseIntegrate because it may be necessary to sum the integrals over the subregions. The value for each of these options is in turn an option or a list of options that should be passed to NIntegrate or NSum:

NPiecewiseIntegrate[
      Mod[(E^(-x)), (E^(-2 x))], {x, 
        0, ∞}, 
NIntegrateOptions -> {WorkingPrecision -> 
            20}, 
NSumOptions -> {NSumTerms -> 
            1000, WorkingPrecision -> 20}]

Even though NPiecewiseSum doesn't have any other options, it uses the same format:

NPiecewiseSum[
      DiscreteDelta[2  i - 3  j - 10] 
        (E^(-Abs[i + j])), {i, -∞, ∞}, {j, -∞, ∞}, 
NSumOptions -> {Method -> Fit, WorkingPrecision -> 20}]

Another method is to set the options for NIntegrate and NSum directly using SetOptions. Incidentally, this is the only way to control the behaviour of NIntegrate when it is called internally by NSum.

Limitations

The Piecewise*/NPiecewise* functions assume that all operations are performed over the domain Reals. Thus in the following example PiecewiseIntegrate will look for subsets of the integration range where Sign is equal to either -1, 0 or 1. All the subsets are empty (except for one point x = 0) and so PiecewiseIntegrate will return 0. In principle it is possible to employ certain heuristics based on FindInstance to detect such cases:

PiecewiseIntegrate[Sign[(Sqrt[x]) + 1], {x, -1, 0}]

Secondly, since the idea of the algorithm is to break the integration/summation range into components, which are then processed separately, it may involve changing the order of the terms in the series (or computing improper integrals as the limit of the integrals over some specific exhausting sequence of subregions). Therefore, this method is guaranteed to work only for absolutely convergent integrals and series. In the next example the value of the sum is ln(2), but rearranging the terms leads to the indeterminate result ∞ - ∞:

PiecewiseSum[((2  Mod[k, 2] - 1)/k), {k, ∞}]

Next, in Mathematica 5 the function Reduce doesn't always convert the original set to the union of disjoint subsets. In such cases summing over those subsets will yield an incorrect result, since some points will be counted more than once. Here is the Reduce output:

Reduce[(k ≠ 0 ∨ n ≠ 0) ∧ 
        0 ≤ k ≤ n ∧ 
        k ∈ Integers]

The values k = 1, 
      2, …, ⌊n⌋ are counted twice, and consequently the corresponding sum comes out wrong:

PiecewiseSum[
      If[k ≠ 0 ∨ n ≠ 0, 1, 0], {k, 0, n}]

Another issue is that currently Sum cannot use assumptions and may not give a generic result. Therefore, the output of PiecewiseSum will not be generically valid either. Here the result given by PiecewiseSum is incorrect for -1 ≤ n < 0:

(∑_(k = -2)^n) (1/k)
PiecewiseSum[(1/k), {k, -2, n}]

Examples

Further examples

More complicated variations of the two of the previous examples:

PiecewiseIntegrate[
      f[x, y, z], {x, 0, 1}, {y, 0, 1}, {z, 0, 
        Max[x + y, a x y]}]

PiecewiseIntegrate[
      UnitStep[-(x^4) + a x + 
          a], {x, -∞, ∞}]

Exampes where the decomposition effectively returns an infinite number of regions:

PiecewiseIntegrate[
      Mod[(E^(-x)), (E^(-2 x))], {x, 0, ∞}]

PiecewiseIntegrate[DiracDelta''[
        x Cos[x]], {x, -∞, ∞}]

PiecewiseSum[
      DiscreteDelta[2  i - 3  j - 10] 
        (E^(-Abs[i + j])), {i, -∞, ∞}, {j, -∞, ∞}]

A more complicated example involving a symbolic parameter:

int = PiecewiseIntegrate[UnitStep[Sin[a x]], {x, 0, π}]

Here's a numerical check; evaluate the next cell to see the graphic. Dashing allows us to see one curve superimposed on another.

Plot[
{int, 
NIntegrate[
          UnitStep[
            Sin[a x]], 
Evaluate@{x, 
              Sequence @@ 
                Range[0, π, π/
                    Abs[a]], π}, 
Method -> 
            GaussKronrod, 
          AccuracyGoal -> 2]}, 
{a, -10, 
        10}, 
PlotStyle -> {Hue[0], {Hue[2/3], 
            Dashing[{ .05,  .05}]}}]
Clear[int]

Notice that we specified the potentially problematic points in the iterator of NIntegrate, namely the points π k/a where the integrand is discontinuous. In that case it is preferable to specify the integration method as well, otherwise Mathematica 5.1 will still use the new methods for numerical integration of piecewise functions, thus negating the speed increase. Additionally, we used AccuracyGoal rather than PrecisionGoal, because it is impossible to attain any level or precision when the exact value of the integral is zero (in other words, we cannot obtain any significant digits). Value of AccuracyGoal -> 2 is sufficient for the plot.

Antiderivative

The antiderivative of the Floor[(x^2)] function:

indef = PiecewiseIntegrate[Floor[(ξ^2)], {ξ, 0, x}]

Since the antiderivative is continuous by construction, it can be used to evaluate definite integrals:

def1 = (indef /. x -> 3) - (indef /. 
            x -> -3)

The direct computation of the same definite integral:

def2 = 
      PiecewiseIntegrate[
        Floor[(ξ^2)], {ξ, -3, 3}]
def1 - def2 // FullSimplify
Clear[indef, def1, def2]

Delta function

If \(TraditionalForm\`f(x\&_) = f(x\_1, …, x\_n)\) and \(TraditionalForm\`ϕ(x\&_) = ϕ(x\_1, …, x\_n)\) are scalar-valued functions, then

\(TraditionalForm\`∫\+\(ℝ\^n\)\(δ(\(f(x\&_)\)\/\(\(|\)\(∇\(f(x\&_)\)\)\(|\)\))\) \(ϕ(
              x\&_)\) dx\&_ = ∫\+\(f(x\&_) = 0\)\(ϕ(x\&_)\) dS\),

where the integral on the right is the surface integral over the surface \(TraditionalForm\`f(x\&_) = 0\). As an example, let's take the cone given by the equation z = k (Sqrt[ (x^2) + (y^2)]), with z running from 0 to h. First we compute the norm of the gradient; notice the 'trick' way of using a pattern in the assumptions:

\(Norm[\(∂\_#\((z - k \@\( x\^2 + y\^2\))\) &\) /@ {x, y, z}] // 
      Simplify[#, _ ∈ Reals] &\)

And then find the area of the lateral surface of the cone:

PiecewiseIntegrate[
      DiracDelta[((z - k (Sqrt[ (x^2) + (y^2)]))/(Sqrt[1 + (k^2)]))], {x, -∞, ∞}, {y, -∞, ∞}, {z, 0, 
        h}, 
Assumptions -> {h > 0, k > 0}]

Probability integrals can also be conveniently expressed in terms of the delta function. For example, if two random variables (Subscript[ξ, 1]),  (Subscript[ξ, 2]) have probability density functions (PDF) (Subscript[f, (Subscript[ξ, 1])]) and (Subscript[f, (Subscript[ξ, 2])]), then the PDF of g((Subscript[ξ, 1]),  (Subscript[ξ, 2])) can be written as

\(TraditionalForm\`\(f\_\(g(ξ\_1, \ ξ\_2)\)\)(
        x) = ∫\+\(-∞\)\%∞\(∫\+\(-∞\)\%∞\( δ(
              x - g(τ\_1, τ\_2))\) \(\(f\_\(ξ\_1\)\)(τ\_1)\) \(\(f\_\(ξ\_2\)\)(τ\_2)\) dτ\_1 dτ\_2\)\)

Thus the following integral gives the PDF for the square of the distance from the origin to a point uniformly distributed in the square [0, 1]×[0, 1]:

f[x_] = 
      PiecewiseIntegrate[
        DiracDelta[x - (Subscript[τ, 1]^2) - (Subscript[τ, 2]^2)], {(Subscript[τ, 1]), 0, 
          1}, {(Subscript[τ, 2]), 0, 1}]

The expectation of the distance:

PiecewiseIntegrate[(Sqrt[x]) f[x], {x, -∞, ∞}] // 
        TrigToExp // Simplify
Clear[f]

Sometimes the representation of the result as a sum of UnitStep terms may be preferable to the conditional form. This is the case if the result is to be differentiated as a generalized function; for example, in this way one can obtain the correct expression for the PDF, which will involve the delta function, from a discontinuous cumulative distribution function (CDF). Integration of such PDF will be effectively equivalent to taking the Stieltjes integral with respect to the corresponding CDF.

Consider another example:

PiecewiseIntegrate[DiracDelta'[2  x + y] ϕ[x, 
          y], {x, -∞, ∞}, {y, -∞, ∞}]
PiecewiseIntegrate[DiracDelta'[2  x + y] ϕ[x, 
          y], {y, -∞, ∞}, {x, -∞, ∞}]

The results do not seem to be equal. However, we should remember that we are working with the set of test functions ϕ that have finite support. Therefore, since

\(TraditionalForm\`∫\+\(-∞\)\%∞\((\(d\/dx\) \(ϕ(
                  x, \(-2\) 
                    x)\))\) dx = \(∫\+\(-∞\)\%∞\((\(ϕ\^\((1, 0)\)\)(x, \(-2\) x) - 
                  2 \(\( ϕ\^\((0, 1)\)\)(
                      x, \(-2\) x)\))\) dx = \(\(\(ϕ(
                x, \(-2\) x)\)\( | \_\(-∞\)\%∞\)\) = 
            0\)\)\),

the change of variables (ϕ^((1, 0)))(
          x, -2 x) -> (ϕ^((1, 0)))(-(y/2), y) shows that the two outputs above are indeed equal to each other.

The integration can also be extended to multidimensional Dirac delta functions which are not a direct product of univariate delta functions:

PiecewiseIntegrate[
Derivative[1, 1][
            DiracDelta][(x^2) + (y^2) - 1, 2  x - y - 1] ϕ[x, 
          y], 
{x, -∞, ∞}, {y, -∞, ∞}]

Square triangle picking

We can apply PiecewiseIntegrate to solving the square triangle picking problem. The problem is as follows: given three points independently and uniformly distributed in a unit square, find the expected value of the area of the triangle formed by these points. The area is considered to be nonnegative; the expectation of the signed area is, of course, 0. We will use two simple considerations: firstly, we integrate over the region where the signed area is positive and then multiply by 2, and secondly, we select only the region where for the three points ((Subscript[x, i]), (Subscript[y, i])) (Subscript[x, 1]) > (Subscript[x, 2]) > (Subscript[x, 3]) and then multiply by 3! permutations of (Subscript[x, i]). This significantly simplifies the decomposition of the region of integration and is already enough to obtain the result sufficiently fast. The anwer is 11/144.

Module[{area}, 
area = (1/2) 
            Last[{x2 - x1, y2 - y1, 0}\[Cross]{x3 - x1, y3 - y1, 
                  0}]; 
12 
           PiecewiseIntegrate[
            area UnitStep[area], 
{x1, 0, 1}, {x2, 0, 
              x1}, {x3, 0, x2}, {y1, 0, 1}, {y2, 0, 1}, {y3, 0, 
              1}]
] // Timing

Here's a numerical check. Notice that the default value of the precision goal for MonteCarlo and QuasiMonteCarlo methods is 2.

Module[{area}, 
area = (1/2) 
            Last[{x2 - x1, y2 - y1, 0}\[Cross]{x3 - x1, y3 - y1, 
                  0}]; 
(11/144) - 
          NIntegrate[
            Abs[area], 
{x1, 0, 1}, {x2, 0, 1}, {x3, 0, 
              1}, {y1, 0, 1}, {y2, 0, 1}, {y3, 0, 
              1}, 
Method -> QuasiMonteCarlo, 
            MaxPoints -> (10^6), 
            PrecisionGoal -> 3]
] // Timing

Birds on wire

Here's another interesting probability problem (author: Prof. Marcin Kuczma): N birds (random points) sit on a wire (segment [0, 1]). For each bird find its nearest neighbour and paint the segment of wire between them. Then find the total length L of all painted segments. What is the limit of the expectation of L as N -> ∞?

A few comments are needed here: points are again assumed to be independent uniformly distributed random variables; secondly, some intervals can be painted twice (from left to right and from right to left) — those intervals are counted only once; and finally, the segments to be painted are determined only after all N birds land. I will not give a complete solution here, just the technical computation to obtain an approximation to the PDF for the length of one segment that allows us to find the asymptotics for the moments:

Module[{f, 
        F}, 
F[x_] = (1 - (E^(-N x))) 
          UnitStep[x]; 
f[
          x_] = (∂_x) F[x]^3; 
Assuming[N > 0, 
        1 - (N/3) 
            PiecewiseIntegrate[
              x f[x], {x, -∞, ∞}]]
]

The answer is 7/18. Finding asymptotics for the variance is slightly more complicated (and depends on the model; here it is assumed that for N points there will be (N - 1) painted intervals, some of them having length zero):

Module[{f, 
        F, μ, σ}, 
F[
          x_] = (1 - (E^(-N x))) UnitStep[x]; 
f[
          x_] = (∂_x) F[x]; 
Assuming[
        N > 0, 
F[x_] = 
          2  PiecewiseIntegrate[
                f[(Subscript[τ, 1])] f[(Subscript[τ, 2])] 
                  f[(Subscript[τ, 3])], 
{(Subscript[τ, 1]), -∞, x}, {(Subscript[τ, 2]), (Subscript[τ, 1]), ∞}, {(Subscript[τ, 3]), -∞, (Subscript[τ, 2])}] + 
(1/3) UnitStep[x]; 
f[
            x_] = (∂_x) F[x]; 
μ = 
          PiecewiseIntegrate[
            x f[x], {x, -∞, ∞}]; 
σ = 
          PiecewiseIntegrate[(((x - μ))^2) 
              f[x], {x, -∞, ∞}]]; 
{μ, σ}
]

Here's a numerical simulation. We need to declare the type of the value returned by Mean and Variance because they can operate on structures with arbitrary dimensions.

{μ, σ} = Compile[{n}, 
Module[{L, 
              Llen}, 
L = 
              Sort@Table[Random[], {n}]; 
Llen = L = 
                ListCorrelate[{-1, 1}, 
                  L]; 
Do[
If[
                L[[i]] > 
                    L[[i - 1]] && 
                  L[[i]] > 
                    L[[i + 
                        1]], 
Llen[[i]] = 0], 
{i, 2, 
                n - 2}
]; 
{n*Mean[Llen], 
              n^2*Variance[
                  Llen]}
], 
{{_Mean, _Real}, {_Variance, _Real}}
][10^6]

{μ, σ} - {(7/18), (89/324)}
Clear[μ, σ]

Hypercube slicing

Suppose we take a 4-dimensional hypercube and cut it with a 3-dimensional hyperplane. We want to find the 3-dimensional volume of the resulting cross-section. First we would like to visualize this cross-section. We will do it by rotating the xyzw space to make the normal to the cutting hyperplane collinear to the vector {0, 0, 0, 1}. This means that the cross-section will lie on a hyperplane parallel to the Oxyz subspace. Therefore, CAD (cylindrical algebraic decomposition) will transform the inequalities describing the hypercube together with the equation of the hyperplane to inequalities for x, y, z and an equation w = (Subscript[w, 0]). Then the inequalities for the first three coordinates can be passed to InequalityPlot3D to visualize the solid.

Here is the cross-section of the hypercube 0 ≤ x ≤ 1, 
      0 ≤ y ≤ 1, 0 ≤ z ≤ 1, 
      0 ≤ w ≤ 1 on the hyperplane x + y + z + w = 2. We obtain a regular octahedron.

<< graphics`
Module[{rotm, slice}, 
rotm[{x_, y_, z_, w_}] := 
        Module[{u, v, s, t}, 
{u, v, s, t} = 
            Norm /@ {{x, y}, {x, y, z}, {x, y, z, w}, {z, 
                  w}}; 
If[
            TrueQ[x == y == 0], 
If[
              TrueQ[z == 0], 
IdentityMatrix[4]*
                Sign[w], 
{{1, 0, 0, 
                  0}, 
{0, 1, 0, 
                  0}, 
{0, 0, w, -z}/
                  t, 
{0, 0, z, w}/
                  t}
], 
{{y, -x, 0, 
                  0}/u, 
{x*z, y*z, -x^2 - y^2, 
                  0}/(u*v), 
{x*w, y*w, 
                  z*w, -x^2 - y^2 - z^2}/(v*
                    s), 
{x, y, z, w}/
                s}
]]; 
slice[$rgn_, 
          Lvar_] := Module[{rgn = $rgn, eq, nv}, 
eq = 
            First@Cases[
                rgn, _Equal]; 
nv = Coefficient[
                  Subtract @@ eq, #, 1] & /@ Lvar; 
rgn = 
            rgn /. Thread[
                Lvar -> 
                  Transpose[rotm[nv]] . Lvar]; 
rgn = 
            Reduce[rgn, Lvar, Reals] /. 
              Lvar[[-1]] == _ -> True; 
InequalityPlot3D[rgn, 
            Sequence @@ List /@ 
                Most@Lvar]
]; 
slice[
        0 ≤ x ≤ 1 && 0 ≤ y ≤ 1 && 
          0 ≤ z ≤ 1 && 0 ≤ w ≤ 1 && 
          x + y + z + w == 2, {x, y, z, 
          w}]
]

The octahedron can be parametrized as 𝕣(x, y, z) = {x, y, z, 
          2 - x - y - z}. Its volume (in the sense of a 3D solid) can be viewed as the three-dimensional area in the xyzw space. For an (n - 1)-dimensional surface in (ℝ^n) the area element dS can be expressed in terms of a vector cross product:

𝕣 = 𝕣(
          (Subscript[u, 1]), …, (Subscript[u, n - 1])),

\(TraditionalForm\`dS = \(\(|\)\(𝕣\_\(u\_1\)\%′×𝕣\_\(u\_2\)\%′×…×𝕣\_\(u\_\(n - 1\)\)\%′\)\(|\)\(du\_1  … du\_\(n - 1\)\)\)\)

\(Norm[
      Cross @@ \(\((∂\_#{x, y, z, 2 - x - y - z} &)\) /@ {x, y, 
            z}\)]\)

Then we integrate dS over the region 0 ≤ x ≤ 1, 
      0 ≤ y ≤ 1, 0 ≤ z ≤ 1, 
      0 ≤ 2 - x - y - z ≤ 1:

PiecewiseIntegrate[
      If[0 ≤ 2 - x - y - z ≤ 1, 2, 0], {x, 0, 1}, {y, 
        0, 1}, {z, 0, 1}]

The same result can be obtained using the formula for the surface integral from the above section on the delta function:

\(PiecewiseIntegrate[
      DiracDelta[\(2 - x - y - z - w\)\/Norm[\(∂\_#\((2 - x - y - z - w)\) &\) /@ {x, y, z, w}]], 
{x, 0, 1}, {y, 0, 1}, {z, 0, 
        1}, {w, 0, 1}]\)

The volume of the octahedron is 4/3 (a geometrical way to see this is to notice that the distance between its opposite vertices — for example, {1, 1, 0, 0} and {0, 0, 1, 1} — equals 2). According to a theorem by Keith Ball, the maximal possible volume is (Sqrt[2]). The maximum can be achieved by taking the hyperplane x - y == 0, which yields a cuboid with the dimensions 1×1×(Sqrt[2]) in the cross-section.

Implementation

To create a standalone package, you can either select File->Save as Special->Package Format from the main menu or select the cell below, click the right mouse button, select Copy As->Plain Text and then paste the selection into any text editor and save it in the plain text format.

Off[General::spell, 
      General::spell1]

BeginPackage["Piecewise`"];

PiecewiseIntegrate::usage = "PiecewiseIntegrate[f, {x, xmin, xmax}]
gives the definite integral of f with respect to x from xmin to xmax.
PiecewiseIntegrate[f, {x, xmin, xmax}, ...] gives the multiple definite \
integral.";

\(\(PiecewiseSum::usage = "PiecewiseSum[f, {i, imin, imax}]
evaluates the sum of f as i goes from imin to imax.
{i, imax} iterates from 1 to imax. {i, imin, imax, di} uses steps di.
PiecewiseSum[f, {i, xmin, xmax}, ...] evaluates the multiple sum.";\)
\)
NPiecewiseIntegrate::usage = "NPiecewiseIntegrate[f, {x, xmin, \
xmax}]
gives a numerical approximation to the definite integral of f with respect to \
x from xmin to xmax.
NPiecewiseIntegrate[f, {x, xmin, xmax}, ...] approximates the multiple \
definite integral.";

NPiecewiseSum::usage = "NPiecewiseSum[f, {i, imin, imax}]
numerically evaluates the sum of f as i goes from imin to imax.
{i, imax} iterates from 1 to imax. {i, imin, imax, di} uses steps di.
NPiecewiseSum[f, {i, xmin, xmax}, ...] numerically evaluates the multiple \
sum.";

ImplicitConditions::usage = "ImplicitConditions is an option for \
PiecewiseIntegrate and PiecewiseSum
which specifies whether to allow unsolved equations or inequalities in the \
generated conditions.";

Parameters::usage = "Parameters is an option for PiecewiseIntegrate \
and PiecewiseSum
which specifies the independent parameters in the integrand/summand and the \
limits.";

\(\(NIntegrateOptions::usage = "NIntegrateOptions is an option for NPiecewiseIntegrate
which specifies the list of options to be passed to NIntegrate.";\)
\)
\(\(NSumOptions::usage = "NSumOptions is an option for NPiecewiseIntegrate and NPiecewiseSum
which specifies the list of options to be passed to NSum.";\)
\)
Begin["`Private`"];
{System`Boole, System`Clip, 
        System`Piecewise};

PiecewiseIntegrate::ilim = "Invalid integration variable or \
limit(s) in `1`.";
PiecewiseIntegrate::region = "The region of integration could not \
be broken down into cylinders.";
PiecewiseIntegrate::indet = "The integrand `1` is indeterminate \
over some regions of the integration range.";
PiecewiseIntegrate::imc = "All constants in the input should be \
real.";
PiecewiseIntegrate::ddz = "The argument of DiracDelta is zero on an \
interval.";
PiecewiseIntegrate::ddlin = "The integrand `1` is not linear in \
`2`.";
PiecewiseIntegrate::ddsup = "The integrand `1` contains a product \
of DiracDelta factors
(or a multidimensional delta function) with intersecting supports.";
PiecewiseIntegrate::nonopt = "`1` is not an option or a list of \
options.";

PiecewiseSum::ilim = "Invalid iteration variable or range in \
`1`.";
PiecewiseSum::region = "The summation range could not be broken \
down into cylinders.";
PiecewiseSum::indet = "The summand `1` is indeterminate over some \
regions of the summation range.";
\(\(PiecewiseSum::irng = "The iteration range `1` is indeterminate.";\)\)
\(\(PiecewiseSum::imc = "All constants in the input should be real.";\)\)
PiecewiseSum::nonopt = "`1` is not an option or a list of \
options.";

NPiecewiseIntegrate::ilim = "Invalid integration variable or \
limit(s) in `1`.";
NPiecewiseIntegrate::region = "The region of integration could not \
be broken down into cylinders.";
NPiecewiseIntegrate::indet = "The integrand `1` is indeterminate \
over some regions of the integration range.";
NPiecewiseIntegrate::imc = "All constants in the input should be \
real.";
NPiecewiseIntegrate::nonopt = "`1` is not an option or a list of \
options.";

NPiecewiseSum::ilim = "Invalid iteration variable or range in \
`1`.";
NPiecewiseSum::region = "The summation range could not be broken \
down into cylinders.";
NPiecewiseSum::indet = "The summand `1` is indeterminate over some \
regions of the summation range.";
NPiecewiseSum::irng = "The iteration range `1` is \
indeterminate.";
NPiecewiseSum::imc = "All constants in the input should be \
real.";
NPiecewiseSum::nonopt = "`1` is not an option or a list of options.\
";

Options[
          PiecewiseIntegrate] = 
{Assumptions :> $Assumptions, ImplicitConditions -> False, 
          Parameters -> Automatic};
Options[
          PiecewiseSum] = 
{Assumptions :> $Assumptions, ImplicitConditions -> False, 
          Parameters -> Automatic};
Options[
          NPiecewiseIntegrate] = 
{NIntegrateOptions :> Options[NIntegrate], 
          NSumOptions :> 
            Options[NSum]};
Options[
          NPiecewiseSum] = 
{NSumOptions :> 
            Options[NSum]};

PiecewiseIntegrate[Sarg___] := 
        Module[{ans}, 
ans /; (ans = 
                  parse[PiecewiseIntegrate, 
                    Sarg]) =!= $Failed];
PiecewiseSum[Sarg___] := 
        Module[{ans}, 
ans /; (ans = 
                  parse[PiecewiseSum, 
                    Sarg]) =!= $Failed];
NPiecewiseIntegrate[Sarg___] := 
        Module[{ans}, 
ans /; (ans = 
                  parse[NPiecewiseIntegrate, 
                    Sarg]) =!= $Failed];
NPiecewiseSum[Sarg___] := 
        Module[{ans}, 
ans /; (ans = 
                  parse[NPiecewiseSum, 
                    Sarg]) =!= $Failed];

\(parse[op_, Sarg___] := 
      Module[
{msg, varQ, Larg = {Sarg}, Liter, Lopt, fx, 
          validopt, 
assum, impcondF, Lparam, Lnintopt, 
          Lnsumopt}, 
msg = 
          Message[MessageName[op, #], ##2] &; 
varQ = 
          Variables@# === {#} &; 
validopt[
            PiecewiseIntegrate] = 
          Assumptions | ImplicitConditions | 
            Parameters; 
validopt[PiecewiseSum] = 
          Assumptions | ImplicitConditions | 
            Parameters; 
validopt[NPiecewiseIntegrate] = 
          NIntegrateOptions | NSumOptions; 
validopt[
            NPiecewiseSum] = 
          NSumOptions; 
Catch[
If[
            Length@Larg == 0, 
msg["argm", op, 0, 
              2]; Throw[$Failed]
]; 
If[
            Length@Larg == 1, 
msg["argmu", op, 
              2]; Throw[$Failed]
]; 
If[\(! FreeQ[Larg, 
                Complex]\), 
msg["imc"]; 
            Throw[$Failed]
]; 
{fx, 
              Larg} = {First@Larg, Rest@Larg}; 
{Lopt, 
              Liter} = 
            Reverse /@ \((Reverse@
                    Larg /. 
{s___, 
                      s2___?\((\(! OptionQ[#]\) &)\)} :> {{s}, {s2}})\); 
If[
            Liter === {}, 
msg["ilim", 
              Lopt[[1]]]; 
            Throw[$Failed]
]; 
Scan[
            If[
\((op === PiecewiseIntegrate || 
                        op === 
                          NPiecewiseIntegrate)\) && 
\(! \((ListQ@# && Length@# == 3 && 
                          varQ@#[[1]])\)\) || 
\((op === PiecewiseSum || 
                        op === 
                          NPiecewiseSum)\) && 
\(! \((ListQ@# && 2 ≤ Length@# ≤ 4 && 
                          varQ@#[[1]])\)\), 
msg["ilim", #]; 
                Throw[$Failed]
] &, 
            Liter]; 
Lopt = 
            Flatten@{Lopt, Options[op]}; 
Scan[
            If[\(! OptionQ[#]\), 
msg["nonopt", #]; 
                Throw[$Failed]
] &, 
            Lopt]; 
Scan[
            If[\(! MatchQ[#[[1]], 
                    validopt[
                      op]]\), 
msg["optx", #[[1]], HoldForm@op[Sarg]]; 
                Throw[$Failed]
] &, 
            Lopt]; 
Lopt = 
            Flatten@{Lopt, 
Assumptions -> True, 
                ImplicitConditions -> False, 
                Parameters -> {}, 
NIntegrateOptions -> {}, NSumOptions -> {}}; 
{assum, impcondF, 
              Lparam, Lnintopt, 
              Lnsumopt} = 
{Assumptions, 
                ImplicitConditions, 
                Parameters, 
NIntegrateOptions, 
                NSumOptions} /. 
              Lopt; 
If[\(! MatchQ[impcondF, 
                True | False]\), 
msg["opttf", 
Sequence @@ \(Cases[
                    Lopt, 
HoldPattern[
                      ImplicitConditions -> _], 
{1}, 
                    1]\)[[1]]]; 
Throw[$Failed]
]; 
If[
            Lparam === Automatic, 
            Lparam = {}]; 
If[\(! ListQ@Lparam\), 
            Lparam = {Lparam}]; 
If[\(! VectorQ[Lparam, 
                varQ]\), 
msg["optvg", 
Sequence @@ \(Cases[
                    Lopt, 
HoldPattern[
                      Parameters -> _], 
{1}, 
                    1]\)[[1]], 
"an expression or a list of expressions that can be used as variables"]; 
Throw[$Failed]
]; 
piecewise[op, fx, 
            Liter, {Lparam, impcondF, assum}, {Lnintopt, 
              Lnsumopt}]
]\ ]
\)
\(piecewise[$op_, $fx_, $Liter_, 
{$Lparam_, $impcondF_, $assum_}, 
{$Lnintopt_, $Lnsumopt_}, $rec_:  
          0] := Module[
{init, getcondfx, convdd, chvar, 
          pwunwrap, pwstep, 
pwreduceop, rgnop, numop, 
          intdd, cad, getparam, pwrefine, msg, 
Lpwhead, 
          diracdelta, element, if, which, piece, 
          K}, 

Lpwhead = 
          Alternatives @@ 
{UnitStep, Sign, Abs, Min, 
              Max, 
Floor, Ceiling, Round, IntegerPart, 
              FractionalPart, Quotient, 
              Mod, 
DiracDelta, \(Derivative[_]\)[
                DiracDelta], DiscreteDelta, 
              KroneckerDelta, 
Element, NotElement, if, 
              which, Boole, Clip, 
              piece}; 

init[] := 
          Module[
{Lparam, Lvar, 
              LLiter}, 
If[$rec == 
                  3 || $fx === $Failed, 
              Throw[$Failed]\ ]; 
If[$Liter === {}, 
              Throw[$fx]]; 
Scan[
              If[\(! FreeQ[Rest@#, 
                      First@#]\), 
Throw[$Failed]
] &, $Liter]; 
LLiter = $Liter /. {s__, 
                      s2___} /; 
{s2} === {} || \(! FreeQ[{s}, {s2}[[1, 
                            1]]]\) :> 
{{s}, {s2}}; 
If[
              LLiter[[2]] =!= {}, 
Throw[
                piecewise[$op, 
piecewise[$op, $fx, 
                    LLiter[[2]], 
{$Lparam, $impcondF, $assum}, {$Lnintopt, $Lnsumopt}, $rec], 
LLiter[[1]], 
{$Lparam, $impcondF, $assum}, {$Lnintopt, $Lnsumopt}, $rec]\ ]
]; 
Lvar = $Liter[[All, 1]]; 
Lparam = 
              Complement[$Lparam, Lvar]; 
SetAttributes[K, 
              NHoldFirst]; 
{{}, Lparam, 
              Lvar}
]; 

getcondfx[Kpv_] := 
          Module[
{fx = $fx, Liter = $Liter, Lcondfx, 
              Lirgn, rev}, 
rev = 
              If[$op === PiecewiseIntegrate || $op === 
                    NPiecewiseIntegrate, \(-1\), 0]; 
{fx, 
                Liter} = {fx, Liter} //. 
{If -> if, 
                  Which -> which, 
                  Piecewise -> piece}; 
{fx, 
                Liter} = \(If[\(! FreeQ[#, 
                        Integrate | 
                          Sum]\), 
#, 
Refine[#, 
                      Join @@ Kpv ∈ 
                        Reals]
] &\) /@ {fx, 
                  Liter}; 
If[\(! FreeQ[Liter, 
                  DiracDelta]\), 
              Throw[$Failed]\ ]; 
If[\(! FreeQ[fx, 
                  DiracDelta]\), 
fx = 
                convdd[fx, 
                  Kpv[[3]]]
]; 
If[$op === PiecewiseSum || $op === 
                  NPiecewiseSum, 
{fx, Liter} = 
                chvar[fx, 
                  Liter]
]; 
Lirgn = 
              Distribute[
Function[{x, x1, 
                      x2}, 
Switch[
                      Refine[Sign[x2 - x1]], 
0 | 
                        1, {{x1 ≤ x ≤ x2, 
                          1}}, 
\(-1\), {{x2 ≤ 
                            x ≤ x1, 
                          rev}}, 
_, {{x1 ≤ 
                            x ≤ x2, 
                          1}, {x2 ≤ x ≤ x1 && 
                            x1 ≠ x2, rev}}
]] @@@ 
                  Liter, List]; 
Lcondfx = \({And @@ # && \(($Assumptions /. _ ∈ Integers -> 
                            True)\), 
fx*
                      Times @@ #2} &\) @@@ 
\(Transpose /@ 
                  Lirgn\); 
{Lcondfx, 
              ReplacePart[Kpv, 
                Liter[[All, 1]], 
                3]}
]; 

chvar[$$fx_, $$Liter_] := 
          Module[
{fx = $$fx, Liter = $$Liter, Lvar, 
              Ldep, 
len, x, a, b, d, i, 
              j}, 
len = 
              Length@Liter; 
Liter = 
              Replace[Liter, 
{{x_, a_} :> {x, 
                      1, a, 1}, 
{x_, a_, 
                      b_} :> {x, a, b, 
                      1}}, 
{1}]; 
Do[
{x, a, b, d} = 
                Liter[[
                  i]]; 
If[
                d == 
                    0 || 
{a, 
                      b} === {\(-∞\), \(-∞\)} || 
{a, 
                      b} === {∞, ∞} || 
a === \(-∞\) && b === ∞ && 
                    N[d] ≠ 
                      1.  || 
a === ∞ && 
                    b === \(-∞\) && 
                    N[d] ≠ \(-1. \) || 
MemberQ[{a, b}, _DirectedInfinity] && \(! NumericQ[
                        d]\), 
msg["irng", $$Liter[[i]]]; 
                Throw[$Failed]
]; 
If[\(! MatchQ[a, _DirectedInfinity]\), 
Liter[[
                    i]] = {x, 
                    0, \((b - a)\)/d}; 
{fx, 
                    Liter[[Range[i + 1, 
                        len]]]} = 
{fx, 
                      Liter[[Range[i + 1, 
                          len]]]} /. 
                    x -> a + d*x, 
If[\(! MatchQ[
                      b, _DirectedInfinity]\), 
Liter[[i]] = {x, 
                      0, \((b - a)\)/d}; 
{fx, 
                      Liter[[Range[i + 1, 
                          len]]]} = 
{fx, 
                        Liter[[Range[i + 1, 
                            len]]]} /. 
                      x -> 
                        b - d*x, 
Liter[[
                      i]] = {x, \(-∞\), ∞}
]\ ], 
{i, 
                len}
]; 
Lvar = 
              Liter[[All, 
                1]]; 
If[
              FreeQ[fx, Integrate | Sum], 
fx = 
                Refine[fx, 
                  Lvar ∈ 
                    Integers]
]; 
fx = 
              fx /. \((h : Sin | Cos | Tan | Cot | Sec | Csc)\)[
                    x_] :> 
With[{tmp = 
                        Refine[h[Expand[x]], Lvar ∈ Integers]}, 
                    tmp /; True]; 
Liter = 
              FixedPoint[
                Replace[#, 
{x_, 0, 
                          Floor[a_] + b_. } /; 
Refine[
                          b ∈ Integers] :> {x, 0, 
                        a + b}, 
{1}] &, 
Expand[Liter]\ ]; 
{fx, Liter} = 
              Together[{fx, Liter}]; 
If[$op === 
                NPiecewiseSum, 
              Return[{fx, 
                  Liter}]\ ]; 
Do[
If[
                MemberQ[
                    Drop[Liter, 
                      i], _DirectedInfinity, {2}] || 
\(! NumericQ[Liter[[i, 
                        3]]]\), 
Continue[]
]; 
Ldep = 
                Table[
getparam[
                    Rest@Liter[[j]], 
                    Take[Lvar, {i, j - 1}]], 
{j, i + 1, 
                    len}]; 
If[
                MatchQ[Ldep, {{} ... }], 
fx = 
                  Sum @@ {fx, 
                      Sequence @@ 
                        Take[Liter, {i, len}]}; 
Liter = 
                  Take[Liter, 
                    i - 1]; 
Break[]
], 
{i, len}
]; 
If[
              Liter === {}, 
Throw[\(fx /. 
Hold[f_] /; MatchQ[f, Lpwhead] :> 
                      f\) /. 
{if -> If, 
                    which -> Which, 
                    piece -> 
                      Piecewise}]
]; 
{fx, 
              Liter}
]; 

convdd[$$fx_, Lvar_] := 
          Module[
{fx = $$fx, Ldd, ddx, 
              ddh}, 
ddh = 
              DiracDelta | \(Derivative[_]\)[
                  DiracDelta]; 
If[\(! FreeQ[
                  fx, \((DiracDelta | \(Derivative[___]\)[
                          DiracDelta])\)[]]\), 
Throw[$Failed]
]; 
fx = fx /. 
{DiracDelta[s__] :> 
                    Times @@ \(DiracDelta /@ {s}\), 
\(\(Derivative[s__]\)[DiracDelta]\)[
                      s2__] :> 
\((If[
                        Length@{s} =!= Length@{s2}, 
                        Throw[$Failed]\ ]; 
Times @@ 
                        MapThread[
\(\(Derivative[#]\)[
                                DiracDelta]\)[#2] &, {{s}, {s2}}])\)\ }; 
fx = 
              Collect[fx, ddh[_], Together]; 
Ldd = 
              Cases[fx, 
                ddh[x_] /; \(! FreeQ[x, ddh[_]]\) :> 
                  x, {0, \(-1\)}]; 
If[
              Ldd =!= {}, 
msg["ddlin", $fx, 
                Ldd[[1]]]; 
              Throw[$Failed]
]; 
fx = 
              fx /. x : ddh[_] /; 
                    FreeQ[x, 
                      Alternatives @@ 
                        Lvar] :> 
\((x /. \((h : ddh)\)[y_] :> \(Hold[h]\)[y])\); 
Ldd = 
              Cases[fx, ddh[_], {0, \(-1\)}]; 
If[
              Ldd === {}, 
              Return[fx]\ ]; 
If[$op =!= 
                  PiecewiseIntegrate || 
\(! FreeQ[
                    Ldd[[All, 0]], 
                    Alternatives @@ 
                      Lvar]\), 
Throw[$Failed]
]; 
ddx = \(Select[
                  Lvar, \(! FreeQ[
                        Ldd, #]\) &]\)[[\(-1\)]]; 
fx = 
              fx /. k_*\((h : ddh)\)[y_] /; \(! FreeQ[y, 
                        ddx]\) :> 
\((k /. 
                        f_[s__] /; 
                            MatchQ[f, Lpwhead] :> \(Hold[f]\)[
                            s])\)*h[y]; 
If[\(! FreeQ[fx, 
                  x_*y_ /; 
\(! FreeQ[x, 
                          ddh[_?\((\(! FreeQ[#, 
                                      ddx]\) &)\)]]\) && 
\(! FreeQ[y, 
                          ddh[_?\((\(! FreeQ[#, 
                                      ddx]\) &)\)]]\)\ ]\), 
Throw[$Failed]
]; 
fx
]; 

pwunwrap[$Lcondfx_, $Kpv_] := 
          Module[
{pwQ, Lcondfx = $Lcondfx, 
              Kpv = $Kpv, 
Scondfx, Lcondfxcur, Lnestop, 
              Lnestvar, pw, h, i}, 
pwQ = 
              MatchQ[Head@#, 
                    Lpwhead] && 
\(! FreeQ[#, 
                      Alternatives @@ \(Join @@ 
                          Kpv\), \(-1\)]\) &; 
If[$op === 
                  PiecewiseIntegrate && 
FreeQ[
                  Lcondfx, _?pwQ], 
Throw[
                If[MemberQ[$Liter, {_, a_, b_} /; 
                      Together[b - a] == 
                        0], 
0, 
Assuming[
And @@ \(\((#[[2]] ≠ #[[3]] &)\) /@ $Liter\), 
Integrate[$fx, 
                      Sequence @@ $Liter]]
]\ ]\ ]; 
Lcondfx = 
              Lcondfx /. \((h : DiracDelta | \(Derivative[_]\)[DiracDelta])\)[
                    x_] :> 
h[
                    x /. Abs[y_] :> 
                        Sqrt[y^2]]; 
If[\(! FreeQ[
                  Lcondfx, \((DiracDelta | \(Derivative[_]\)[DiracDelta])\)[
                      x_] /; 
\(! FreeQ[x, 
                        Lpwhead]\)\ ]\), 
Throw[$Failed]
]; 
Lnestop = 
              Cases[Lcondfx, 
\((Integrate | Sum | 
                        NIntegrate | NSum | Product | 
                        NProduct | 
PiecewiseIntegrate | 
                        PiecewiseSum | 
NPiecewiseIntegrate | NPiecewiseSum)\)[
                    s__] :> {s}, \(-1\)]; 
Scan[
\((Lnestvar = \(Rest[#]\)[[All, 
                      1]]; 
If[\(! FreeQ[#, Lpwhead[s__] /; 
\(! FreeQ[{s}, 
                              Alternatives @@ 
                                Lnestvar]\)]\), 
Throw[$Failed]
])\) &, 
Lnestop]; 
Lcondfx = \(Lcondfx //. 
{if[a_, b_] :> 
                      if[a, b, 
                        Indeterminate], 
which[
                          s__] /; {s}[[\(-2\)]] =!= True :> 
                      which[s, True, 
                        Indeterminate]}\) /. 
\((h : 
                          Element | NotElement)\)[x_, 
                      y_] /; 
MatchQ[Head@x, 
                      List | Alternatives] :> 
And @@ \(\((h[#, y] &)\) /@ 
                      x\); 
Scondfx = 
              Fold[List, {}, Lcondfx]; 
Lcondfx = 
              h[]; 
While[
              Scondfx =!= {}, 
If[\(! Scondfx[[2, 1]]\) || 
                  Scondfx[[2, 
                      2]] == 
                    0, 
Scondfx = 
                  Scondfx[[1]]; 
Continue[]
]; 
pw = 
                Cases[Scondfx[[2]], 
                  x_?pwQ :> 
                    HoldPattern@x, {0, \(-1\)}]; 
If[
                pw === {}, 
Lcondfx = 
                  h[Lcondfx, 
                    Scondfx[[2]]]; 
Scondfx = 
                  Scondfx[[1]]; 
Continue[]
]; 
pw = 
                pw[[\(-1\)]]; 
{Lcondfxcur, Kpv} = 
                pwstep[pw[[1]], 
                  Kpv]; 
Lcondfxcur = 
                Table[
{Scondfx[[2, 
                          1]] && 
                        Lcondfxcur[[i, 
                          1]], 
                      Scondfx[[2, 
                        2]]} /. 
pw -> Lcondfxcur[[i, 
                        2]], 
{i, 
                    Length@Lcondfxcur}]; 
Scondfx = 
                Fold[List, Scondfx[[1]], 
                  Lcondfxcur]
]; 
Lcondfx = List @@ Flatten[Lcondfx, Infinity, h]; 
\({Lcondfx, 
                  Kpv} /. 
Hold[f_] /; 
                    MatchQ[f, Lpwhead] :> 
                  f\) /. 
{diracdelta -> DiracDelta, 
                element -> Element, 
if -> If, 
                which -> Which, 
                piece -> 
                  Piecewise}
]; 

pwstep[pw_, Kpv_] := 
          Module[
{condfx, nK, len, i, 
              j}, 
nK = 
              Length@Kpv[[1]]; 
len = Length@pw; 
condfx = 
              Switch[Head@
                  pw, 
UnitStep, 
{{And @@ \(\((# ≥ 0 &)\) /@ pw\), 
                    1}, 
{\(! And @@ \(\((# ≥ 
                                0 &)\) /@ pw\)\), 
                    0}}, 
Sign, 
{{pw[[1]] < 
                      0, \(-1\)}, 
{pw[[1]] == 0, 
                    0}, 
{pw[[1]] > 0, 
                    1}}, 
Abs, 
{{pw[[1]] < 
                      0, \(-pw[[1]]\)}, 
{pw[[1]] ≥ 0, 
                    pw[[1]]}}, 
Min, 
Table[
{And @@ 
                        Table[pw[[i]] < 
                            pw[[j]], {j, 
                            i - 1}] && 
And @@ 
                        Table[pw[[i]] ≤ pw[[j]], {j, i + 1, 
                            len}], 
pw[[i]]}, 
{i, 
                    len}], 
Max, 
Table[
{And @@ 
                        Table[pw[[i]] > 
                            pw[[j]], {j, 
                            i - 1}] && 
And @@ 
                        Table[pw[[i]] ≥ pw[[j]], {j, i + 1, 
                            len}], 
pw[[i]]}, 
{i, 
                    len}], 
Floor, 
\(nK++\); 
{{K[nK] ≤ 
                      pw[[1]] < 
                      K[nK] + 1, 
                    K[nK]}}, 
Ceiling, 
\(nK++\); 
{{K[nK] - 1 < 
                      pw[[1]] ≤ K[nK], 
                    K[nK]}}, 
Round, 
\(nK++\); 
{{K[nK] ≤ 
                      pw[[1]] < 
                      K[nK] + 1/2, 
                    K[nK]}, 
{K[nK] - 1/2 < 
                      pw[[1]] < K[nK], 
                    K[nK]}, 
{pw[[1]] == 2  K[nK] + 1/2, 
                    2  K[nK]}, 
{pw[[1]] == 2  K[nK] - 1/2, 
                    2  K[nK]}}, 
IntegerPart, 
{{pw[[1]] ≥ 0, 
                    Floor[pw[[1]]]}, 
{pw[[1]] < 0, 
                    Ceiling[
                      pw[[1]]]}}, 
FractionalPart, 
{{True, 
                    pw[[1]] - 
                      IntegerPart[
                        pw[[1]]]}}, 
Quotient, 
{{pw[[2]] ≠ 0, 
If[
                      len == 2, 
Floor[
                        pw[[1]]/
                          pw[[2]]], 
Floor[\((pw[[1]] - 
                              pw[[3]])\)/
                          pw[[2]]]\ ]\ }, 
{pw[[2]] == 
                      0, Indeterminate}}, 
Mod, 
{{True, 
                    pw[[1]] - 
                      pw[[2]]*
                        Quotient @@ 
                          pw}}, 
DiracDelta, 
{{pw[[1]] == 0, 
                    diracdelta @@ 
                      pw}, 
{pw[[1]] ≠ 0, 
                    0}}, 
\(Derivative[_]\)[
                  DiracDelta], 
{{pw[[1]] == 
                      0, \(pw[[0, 0]]\)[
                        diracdelta] @@ 
                      pw}, 
{pw[[1]] ≠ 0, 
                    0}}, 
DiscreteDelta, 
{{And @@ \(\((# == 0 &)\) /@ pw\), 
                    1}, 
{\(! And @@ \(\((# == 
                                0 &)\) /@ pw\)\), 
                    0}}, 
KroneckerDelta, 
If[
                  len == 
                    1, 
{{pw[[1]] == 0, 
                      1}, 
{pw[[1]] ≠ 0, 0}}, 
{{Equal @@ pw, 
                      1}, 
{\(! Equal @@ pw\), 
                      0}}\ ], 
Element, 
If[
                  pw[[2]] =!= 
                    Integers, 
{{True, 
                      element @@ 
                        pw}}, 
\(nK++\); 
{{pw[[1]] == 
                        K[nK], 
                      True}, 
{K[nK] < 
                        pw[[1]] < 
                        K[nK] + 1, 
                      False}}\ ], 
NotElement, 
If[
                  pw[[2]] =!= 
                    Integers, 
{{True, \(! element @@ 
                          pw\)}}, 
\(nK++\); 
{{pw[[1]] == 
                        K[nK], 
                      False}, 
{K[nK] < 
                        pw[[1]] < 
                        K[nK] + 1, 
                      True}}\ ], 
if, 
{{pw[[1]], 
                    pw[[2]]}, {\(! pw[[1]]\), 
                    pw[[3]]}}, 
which, 
Table[
{pw[[i]] && 
                      And @@ Table[\(! pw[[
                              j]]\), {j, 1, i - 2, 
                            2}], 
pw[[
                      i + 1]]}, 
{i, 1, 
                    len - 1, 
                    2}], 
Boole, 
{{pw[[1]], 
                    1}, {\(! pw[[1]]\), 
                    0}}, 
Clip, 
Switch[
                  len, 
1, {{pw[[1]] < \(-1\), \(-1\)}, {\(-1\) ≤ 
                        pw[[1]] ≤ 1, pw[[1]]}, {1 < 
                        pw[[1]], 
                      1}}, 
2, 
{{pw[[1]] < 
                        pw[[2, 1]], 
                      pw[[2, 
                        1]]}, 
{pw[[2, 1]] ≤ 
                        pw[[1]] ≤ pw[[2, 2]], 
                      pw[[1]]}, 
{pw[[2, 1]] < 
                          pw[[1]] && 
                        pw[[2, 2]] < 
                          pw[[1]], 
                      pw[[2, 
                        2]]}, 
{pw[[2, 2]] < 
                          pw[[2, 1]] && 
                        pw[[1]] == 
                          pw[[2, 1]], 
                      pw[[2, 
                        2]]}}, 
_, 
{{pw[[1]] < 
                        pw[[2, 1]], 
                      pw[[3, 
                        1]]}, 
{pw[[2, 1]] ≤ 
                        pw[[1]] ≤ pw[[2, 2]], 
                      pw[[1]]}, 
{pw[[2, 1]] < 
                          pw[[1]] && 
                        pw[[2, 2]] < 
                          pw[[1]], 
                      pw[[3, 
                        2]]}, 
{pw[[2, 2]] < 
                          pw[[2, 1]] && 
                        pw[[1]] == 
                          pw[[2, 1]], 
                      pw[[3, 
                        2]]}}\ ], 
piece, 
Switch[
                  len, 
1, 
                  piece[pw[[1]], 
                    0], 
_, {{True, 
                      which @@ Join[
                          Flatten[
                            Reverse /@ 
                              pw[[1]], 
                            1], {True, 
                            pw[[2]]}]}}
]\ ]; 
{condfx, 
              ReplacePart[Kpv, Array[K, nK], 
                1]}
]; 

pwreduceop[$cond_, $$fx_, Kpv_] := 
          Module[
{cond = $cond, fx = $$fx, LKpv, LK, 
              LKnew, Lparam, Lvar, Lcond, Lelem}, 
fx = 
              Together[fx]; 
If[
              cond === False || fx == 0, 
              Return[0]\ ]; 
{LK, Lparam, Lvar} = 
              Kpv; 
LKpv := 
              Join[LK, Lparam, 
                Lvar]; 
cond = \(cond /. 
Abs[x_] :> 
                    Sqrt[x^2]\) /. 
a_ ≠ 
                    b_ :> a < b || b < a; 
LK = 
              Union@Cases[cond, _K, \(-1\)]; 
Lparam = 
              Join[Lparam, getparam[cond, LKpv]\ ]; 
cond = 
              Reduce[cond && 
LK ∈ 
                    Integers && 
If[$op === 
                        PiecewiseSum || $op === NPiecewiseSum, 
                    Lvar ∈ Integers, True], 
LKpv, 
                Reals, 
GeneratedParameters -> K, 
                Backsubstitution -> False]; 
If[
              Head@cond === Reduce, 
msg["region"]; 
              Throw[$Failed]
]; 
If[\(! FreeQ[cond, Unequal | Re | Im | Reals]\), 
cond = 
                Reduce[cond, LKpv, 
                  Reals, 
GeneratedParameters -> K, 
                  Backsubstitution -> 
                    False]
]; 
cond = 
              Refine[cond, _ ∈ 
                    Reals] /. 
a_ ≠ 
                    b_ :> 
                  a < b || b < a; 
LKnew = 
              Union@Cases[cond, _K, \(-1\)]; 
LK = 
              Join[Complement[LKnew, LK], LK]; 
Lparam = 
              Join[Lparam, 
                getparam[cond, LKpv]\ ]; 
Lelem = 
              Union@Cases[
                  cond, _Element, \(-1\)]; 
If[\(! FreeQ[
                    cond, NotElement]\) || 
Complement[
                    Lelem[[All, 
                      2]], {Integers}] =!= {} || 
Complement[
                    Lelem[[All, 1]] /. 
                      Alternatives -> Sequence, 
LK, 
                    If[$op === PiecewiseSum || $op === NPiecewiseSum, 
                      Lvar, {}]\ ] =!= {}, 
msg["region"]; 
              Throw[$Failed]
]; 
cond = 
              cond /. x_ArcTan :> 
                  FullSimplify[x, 
                    TimeConstraint -> 10]; 
cond = 
              LogicalExpand@
                cond; 
Lcond = \(If[Head@# =!= Or, {#}, 
                    List @@ #] &\)@
                cond; 
Total[\(rgnop[#, 
                    fx, {LK, Lparam, Lvar}] &\) /@ 
                Lcond]
]; 

rgnop[cond_, $$fx_, Kpv_] := 
          Module[
{fx = $$fx, LKiter, pcond, Liter, LK, 
              Lparam, Lvar, 
floor, ceiling, n, i, 
              ans}, 
{LK, Lparam, Lvar} = 
              Kpv; 
{LKiter, pcond, Liter} = 
              cad[cond, Kpv]; 
If[\(! pcond\), 
              Return[0]]; 
Assuming[
pcond && 
                LK ∈ Integers && Lparam ∈ Reals && 
                Lvar ∈ Integers, 
{LKiter, 
                  Liter} = \({LKiter, 
                      Liter} //. 
\((h : Floor | Ceiling)\)[
                        x_] :> 
With[{tmp = 
                            h[Expand[x]]}, 
                        tmp /; True]\) //. 
\((h : 
                            Floor | Ceiling)\)[
                        a_.  + b_] /; 
Refine[
                        b ∈ 
                          Integers] :> 
With[{tmp = h[a] + b}, 
                      tmp /; True]
]; 
Do[
n = 
                Together[
                  LKiter[[i, 3]] - 
                    LKiter[[i, 
                      2]]]; 
Which[
                n < 0, fx = 0, 
n == 
                  0, {fx, 
                    Liter} = {fx, 
                      Liter} /. 
LKiter[[i, 1]] -> 
                      LKiter[[i, 
                        2]]], 
{i, 
                Length@LKiter, 
                1, \(-1\)}
]; 
If[\((fx = 
                    Together[fx])\) == 0, 
              Return[0]]; 
{fx, LKiter, 
                Liter} = {fx, LKiter, Liter} /. 
                Abs[x_] :> 
                  Sqrt[x^2]; 
If[\(($op === 
                      PiecewiseIntegrate || $op === 
                      NPiecewiseIntegrate)\) && 
FreeQ[fx, 
                  DiracDelta] && 
MemberQ[
                  Liter, {_, a_, b_} /; 
                    Together[b - a] == 0], 
Return[
                0]
]; 
If[
              fx === Indeterminate, 
              msg["indet", $fx]\ ]; 
If[$op === 
                  NPiecewiseIntegrate || $op === 
                  NPiecewiseSum, 
Return[
                numop[fx, LKiter, pcond, 
                  Liter]\ ]
]; 
ans = 
              Assuming[
If[$op === 
                    PiecewiseIntegrate, 
And @@ \(\((#[[2]] ≤ #[[1]] ≤ #[[3]] &)\) /@ 
                        LKiter\) && 
pcond && 
LK ∈ Reals && 
                    Lparam ∈ Reals, 
LK ∈ 
                      Integers && Lparam ∈ Reals && 
                    Lvar ∈ 
                      Integers
], 
If[$op === PiecewiseIntegrate, 
If[\(! FreeQ[fx, 
                        DiracDelta]\), 
intdd[fx, 
                      Liter], 
Assuming[
                      And @@ \(\((#[[2]] ≠ #[[3]] &)\) /@ 
                          Liter\), 
Integrate[fx, 
                        Sequence @@ 
                          Liter]\ ]
], 
ans = 
                    Sum @@ {fx, 
Liter[[\(-1\)]] /. {Floor -> floor, Ceiling -> 
                                ceiling}} /. 
{floor -> Floor, ceiling -> Ceiling}; 
Fold[
                    Function[{fx, iter}, 
piecewise[
                        PiecewiseSum, 
                        fx, {iter}, 
{If[
                            TrueQ[$impcondF], {}, 
                            Select[
                              Lparam, \(! FreeQ[
                                    iter, #]\) &]\ ], 
$impcondF, $Assumptions}, 
{{}, {}}]\ ], 
ans, 
                    Reverse@\(Most@
                        Liter\)]
]\ ]; 
If[ans === $Failed, 
              Throw[$Failed]\ ]; 
ans = 
              If @@ {pcond, ans, 0}; 
LKiter = 
              Select[LKiter, 
\(! FreeQ[{ans, Liter, 
                        LKiter[[All, {2, 
                            3}]]}, #[[1]]]\) &]; 
ans = 
              Fold[
                piecewise[
                    PiecewiseSum, #, {#2}, 
{{}, $impcondF, 
$Assumptions && LK ∈ Integers && 
                        Lparam ∈ 
                          Reals}, 
{{}, {}}, $rec + 
                      1] &, 
ans, 
                Reverse@LKiter]; 
If[ans === $Failed, 
              Throw[$Failed]\ ]; 
ans
]; 

intdd[$$fx_, Liter_] := 
          Module[
{fx = $$fx, ddh, helddelta, Ldd, pos, x, 
              a, b}, 
ddh = 
              DiracDelta | \(Derivative[_]\)[
                  DiracDelta]; 
If[\(! FreeQ[fx, 
                  ddh[0]]\), 
msg["ddz"]; 
              Return[Indeterminate]
]; 
pos = 
              Position[
                Liter[[All, 
                  1]], 
x_ /; \(! FreeQ[
                      fx, ddh[
                        y_ /; \(! FreeQ[y, 
                              x]\)]]\), 
{1}]; 
If[
              pos === {}, 
Return[
                Integrate[fx, 
                  Sequence @@ 
                    Liter]\ ]
]; 
pos = 
              pos[[\(-1\), 
                1]]; 
{x, a, b} = 
              Liter[[
                pos]]; 
If[
              pos < Length@Liter, 
fx = 
                piecewise[
                    PiecewiseIntegrate, 
fx /. 
                      DiracDelta -> helddelta, 
                    Drop[Liter, 
                      pos], 
{{}, $impcondF, $Assumptions}, {{}, {}}] /. 
helddelta -> 
                    DiracDelta
]; 
Ldd = 
              Cases[fx, 
                ddh[y_] /; \(! FreeQ[y, 
                      x]\), {0, \(-1\)}]; 
If[
              Ldd === {}, 
Return[
                Integrate[fx, 
                  Sequence @@ 
                    Take[Liter, 
                      pos]]\ ]
]; 
If[\(! FreeQ[fx, \((Integrate | Sum | Product)\)[
                      s__] /; 
\(! FreeQ[{s}, 
                          ddh[_]]\) || 
\(! FreeQ[\(Rest[{s}]\)[[All, 1]], 
                          x]\)\ ]\), 
Throw[$Failed]
]; 
Scan[
              If[\(! PolynomialQ[fx, #]\) || 
                    Exponent[fx, #] > 1, 
msg["ddlin", 
                    fx, #]; Throw[$Failed]
] &, 
              Ldd]; 
If[a ≠ b, 
              Return[Indeterminate]\ ]; 
fx = 
              Collect[fx, Ldd, Together]; 
fx -= 
              First@\(Flatten@
                  CoefficientList[fx, Ldd]\); 
fx = 
              fx /. k_. *h_[y_] /; 
                    MemberQ[Ldd, 
                      h[y]] :> 
Module[{L, 
                      ord}, 
L = 
                      D[#/D[y, x], x] &; 
ord = 
                      If[h === DiracDelta, 0, 
                        h[[0, 
                          1]]]; 
If[\(! IntegerQ[ord]\), Throw[$Failed]\ ]; 
\(\((\(-1\))\)^ord*
                          Nest[L, k /. DiracDelta -> helddelta, 
                              ord]/
Abs[Together[D[y, x]]] /. 
                        x -> a\) /. 
helddelta -> 
                        DiracDelta
]; 
If[\(! FreeQ[fx, 
                  ddh[0]]\), 
msg["ddsup", $fx]; 
              Return[Indeterminate]
]; 
piecewise[PiecewiseIntegrate, fx, 
              Take[Liter, 
                pos - 1], 
{{}, $impcondF, $Assumptions}, {{}, {}}]
]; 

numop[
            fx_, $LKiter_, pcond_, Liter_] := 
          Module[
{LKiter = $LKiter, f, 
              ans}, 
If[\(! TrueQ[pcond]\), 
              Throw[$Failed]\ ]; 
LKiter = 
              Select[LKiter, 
\(! FreeQ[{fx, Liter, 
                        LKiter[[
                          All, {2, 
                            3}]]}, #[[1]]]\) &]; 
ans = 
              If[$op === NPiecewiseIntegrate, 
If[
                  LKiter === {}, 
NIntegrate @@ {fx, 
                      Sequence @@ Liter, $Lnintopt}, 
f := 
                    NIntegrate @@ {fx, 
                          Sequence @@ 
                            Liter, $Lnintopt} /; 
NumericQ[
                        LKiter[[1, 
                          
                          1]]]; 
NSum @@ {f, Sequence @@ LKiter, $Lnsumopt}
], 
If[
                  LKiter === {}, 
NSum @@ {fx, 
                      Sequence @@ Liter, $Lnsumopt}, 
f := 
                    NSum @@ {fx, 
                          Sequence @@ 
                            Liter, $Lnsumopt} /; 
NumericQ[
                        LKiter[[1, 
                          1]]]; 
NSum @@ {f, Sequence @@ 
                        LKiter, $Lnsumopt}
]\ ]; 
If[\(! FreeQ[ans, NIntegrate | NSum]\), 
              Throw[$Failed]\ ]; 
ans
]; 

cad[cond_, Kpv_] := 
          Module[
{Lcond = cond, LKpv, nK, nparam, 
              nvar, 
Lcondcur, LKiter, pcond, Liter, x, 
              i}, 
Lcond = \(If[Head@# =!= And, {#}, 
                    List @@ #] &\)@Lcond; 
If[
              Lcond === {False}, 
              Return[{{}, False, {}}]\ ]; 
LKpv = 
              Join @@ Kpv; 
{nK, nparam, nvar} = 
              Length /@ Kpv; 
Liter = 
              Table[{}, {Length@LKpv}]; 
Lcond = 
              Replace[Lcond, 
{a_ < b_ | b_ > a_ /; 
                      FreeQ[{a, b}, 
                        Alternatives @@ 
                          LKpv] :> 
Together[
                        a - b] < 
                      0, 
a_ ≤ b_ | 
                        b_ ≥ a_ /; 
                      FreeQ[{a, b}, 
                        Alternatives @@ 
                          LKpv] :> 
Together[
                        a - b] ≤ 
                      0}, 
{1}]; 
Lcond = 
              Map[Together, Lcond, {2}]; 
Lcond = 
              Replace[Lcond, {b_ > a_ :> a < b, 
                  b_ ≥ a_ :> 
                    a ≤ b}, {1}]; 
pcond = 
              And @@ Select[Lcond, 
                  FreeQ[#, 
                      Alternatives @@ 
                        LKpv] &]; 
Do[
x = LKpv[[
                  i]]; 
Lcondcur = 
                Cases[Lcond, 
x < a_ | a_ < x | 
                      x ≤ a_ | a_ ≤ x | x == a_ | 
                      a_ == x /; 
FreeQ[a, 
                      Alternatives @@ 
                        Drop[LKpv, i]]\ ]; 
If[
                nK + 1 ≤ i ≤ 
                  nK + nparam, 
pcond = 
                  pcond && And @@ 
                      Lcondcur; 
Continue[]
]; 
If[
                i ≤ 
                    nK || 
i > 
                      nK + nparam && \(($op === PiecewiseSum || $op === 
                          NPiecewiseSum)\), 
Lcondcur = 
                  Lcondcur /. 
{a_ < x :> 
                        Floor[a] + 1 ≤ x, 
x < 
                          a_ :> 
                        x ≤ 
                          Ceiling[a] - 1, 
a_ ≤ 
                          x :> 
                        Ceiling[a] ≤ 
                          x, 
x ≤ 
                          a_ :> 
                        x ≤ 
                          Floor[a], 
a_ == x | 
                          x == a_ :> 
                        Sequence[Ceiling[a] ≤ x, 
                          x ≤ 
                            Floor[a]]}
]; 
Lcondcur = 
                Lcondcur /. 
{a_ < x :> 
                      a ≤ x, 
x < 
                        a_ :> 
                      x ≤ a, 
a_ == x | 
                        x == a_ :> 
                      Sequence[a ≤ x, 
                        x ≤ 
                          a]}; 
Liter[[
                  i]] = {x, 
Max@
                    Cases[Lcondcur, 
                      a_ ≤ x :> 
                        a], 
Min@
                    Cases[Lcondcur, 
                      x ≤ a_ :> 
                        a]}; 
If[
                MatchQ[Liter[[
                    i]], 
{_, 
                        a_Max, _} | {_, _, 
                        a_Min} /; 
\(! FreeQ[a, 
                        Alternatives @@ 
                          LKpv, \(-1\)]\)\ ], 
Throw[$Failed]
], 
{i, 
                Length@LKpv}
]; 
If[\(! pcond\), Return[{{}, False, {}}]\ ]; 
{LKiter, 
                Liter} = \(Take[
                    
                    Liter, #] &\) /@ {nK, \(-nvar\)}; 
{LKiter, pcond, 
              Liter}
]; 

getparam[$expr_, Lindep_] := 
          Module[
{expr = $expr, 
              Lparam}, 
expr = \(\(Flatten@{expr} /. 
x_Inequality :> 
                      LogicalExpand@x\) /. 
Element | 
                      NotElement -> \((# &)\)\) /. 
Alternatives -> 
                  Sequence; 
Lparam = 
              Variables@
                FixedPoint[
Replace[#, 
s_ /; \(! FreeQ[s, 
                                Alternatives @@ 
                                  Lindep, \(-1\)]\) || 
\(! AtomQ@s\) && Variables@s === {} :> 
Sequence @@ 
                          s, 
{1}] &, 
expr]; 
DeleteCases[
              Lparam, 
Alternatives @@ 
                  Lindep | 
True | False | 
                Indeterminate | _Slot | _Function]
]; 

pwrefine[$expr_] := 
          Module[
{expr = $expr, 
              simp}, 
simp = 
              Simplify[#, Assumptions -> True, 
                    TimeConstraint -> 
                      5] /. 
{b_ > a_ :> a < b, 
                    b_ ≥ a_ :> 
                      a ≤ b} &; 
expr = 
              expr /. {If -> if, Which -> which, 
                  Piecewise -> piece}; 
expr = 
              expr //. k_. *if[cond_, e_, 0] :> 
                  if[simp[cond], k*e, 0]; 
expr = 
              expr //. if[cond_, k_. *\((if[cond2_, e_, 0] + a_. )\) + b_. , 
                      0] /; 
Count[{a, 
                          b}, _if, \(-1\)] ≤ 
                      4 :> 
if[
                      simp[cond && cond2], k*\((e + a)\) + b, 
                      0] + 
if[simp[cond && \(! cond2\)], 
                      k*a + b, 0]; 
expr = 
              expr //. if[a_ < b_, e_, 0] + 
                    if[a_ ≤ b_, e2_, 
                      0] :> 
if[a < b, e + e2, 
                      0] + if[a == b, e2, 0]; 
expr = 
              expr //. if[cond_, e_, 0] + 
                    if[cond_, e2_, 0] :> 
if[
                    cond, e + e2, 0]; 
expr = 
              expr //. if[False, _, 0] -> 0; 
expr = 
              expr /. if[cond_, e_, 0] :> 
if[
                    cond, Together[Refine[e, cond]], 
                    0]; 
expr = 
              expr /. if[_, 0, 0] | if[False, _, 0] -> 
                  0; 
Together[
                Refine[expr]] /. 
{if -> If, 
                which -> Which, 
                piece -> 
                  Piecewise}
]; 

msg = 
          Message[MessageName[$op, #], ##2] &; 

Module[
{Kpv, Lcondfx, assum = $assum, 
            Lmsg, msgonF, ans}, 
Lmsg = 
            Hold /@ Unevaluated@
{Reduce::nsmet, 
                  Reduce::ratnz, Reduce::naqs, 
                  FullSimplify::infd, 
General::ivar, 
                  General::cas, General::infy, 
                  General::indet}; 
\(\((msgonF[#] = 
                  StringQ @@ #)\) &\) /@ 
            Lmsg; 
CheckAbort[
Off @@@ 
              Lmsg; 
If[ListQ@assum, 
              assum = And @@ 
                  assum]; 
Block[
{$Assumptions = 
                  assum}, 
ans = 
                Catch[
Kpv = 
                    init[]; 
{Lcondfx, Kpv} = 
                    getcondfx[Kpv]; 
{Lcondfx, Kpv} = 
                    pwunwrap[Lcondfx, Kpv]; 
ans = 
                    Total[\(pwreduceop[##, Kpv] &\) @@@ 
                        Lcondfx]; 
If[\(! FreeQ[ans, 
                        Integrate | Sum]\), ans, 
                    pwrefine[
                      ans]]
]\ ]; 
\(If[
                  msgonF[#], On @@ #] &\) /@ 
              Lmsg; 
ans, 
\(If[
                  msgonF[#], On @@ #] &\) /@ 
              Lmsg; 
Abort[]
]\ ]\ ]
\)
End[];
EndPackage[]

On[General::spell, General::spell1]

Credits

Copyright © 2003-2004, Maxim Rytin