Adapting linearExpand from my answer to How to do algebra on unevaluated integrals?, we can come up with some transformations to factor separable multiple integrals. The function someFunction internally deals with and returns Inactive integrals, which can be evaluated with Activate, if appropriate or desired.
Examples
someFunction[Integrate[p[x] p[y], {x, -1, 1}, {y, -1, 1}]]
someFunction[Integrate[p[x] p[y], {x, -1, 1}, {y, -1, 100}]]
someFunction[Integrate[p[x] p[y], {x, -1, 1}, {y, -1, x}]]
someFunction[Integrate[p[x] (p[y] + q[y]), {x, -1, 1}, {y, -1, 1}]^2]
(* Vitaliy Kaurov's example *)
test = Integrate[p[x]*p[y]*q[z]*r[s]*r[u],
{x, -1, 1}, {y, -1, 1}, {z, 0, 2}, {s, 3, 10}, {u, 3, 10}];
someFunction[test]
% /. changeVar[x]
someFunction[Integrate[q[x], {x, 0, 2}] + test]
Outline
We need some auxilliary functions.
iterated converts a multiple integral into an iterated integral.linearExpand applies linearity properties to an integral, distributing the integral over a sum and factoring out constants.factorConstants recursively factors constants out of nested integrals using linearExpand.changeVar changes the variable of integration.
These can be applied by Simplify via the TransformationFunctions option. The trick, and it can be tricky, is devising a ComplexityFunction that will prefer a result in a factored out form. Note that we use changeVar only to combine two integrals that are equivalent. At the end we can change the variable in the independent integrals to the same one.
Code dump
ClearAll[linearExpand, iterated, factorConstants, changeVar, someFunction];
linearExpand[e_, x_, head_] :=
e //. {op : head[arg_Plus, __] :> Distribute[op],
head[arg1_Times, var_, rest___] /; ! FreeQ[var, x] :>
With[{dependencies = Internal`DependsOnQ[#, x] & /@ List @@ arg1},
Pick[arg1, dependencies, False] *
head[Pick[arg1, dependencies, True], var, rest]]};
iterated[Integrate[f_, dom : {_Symbol, _, _} ..]] :=
Fold[Inactive[Integrate], f, Reverse@{dom}];
iterated[Inactive[Integrate][f_, dom : {_Symbol, _, _} ..]] :=
Fold[Inactive[Integrate], f, Reverse@{dom}];
factorConstants[Inactive[Integrate][f_, {v_Symbol, a_, b_}]] :=
Inactive[Integrate][
f /. j : Inactive[Integrate][_, _] :> factorConstants[j], {v, a, b}] /.
i : Inactive[Integrate][_, {v, _, _}] :>
linearExpand[i, v, Inactive[Integrate]];
factorConstants[x_] := x;
changeVar[Inactive[Integrate][f_, {v_, a_, b_}]] :=
i : Inactive[Integrate][g_, {w_, a, b}] /;
Simplify[f == g /. v -> w] :> (i /. w -> v);
changeVar[v_] :=
i : Inactive[Integrate][g_, {w_Symbol, _, _}] :> (i /. w /; FreeQ[g, Integrate] -> v);
someFunction[integral_] :=
With[{v =
FirstCase[Hold[integral],
HoldPattern[(Integrate | Inactive[Integrate])[_, {x_, _, _}, ___]] :> x,
Missing[],
Infinity]},
(Simplify[
integral /.
{i : (Integrate | Inactive[Integrate])[_, {_Symbol, _, _}, {_Symbol, _, _} ..] :>
iterated@Inactivate[i, Integrate],
i : (Integrate | Inactive[Integrate])[_, {_Symbol, _, _}] :>
Inactivate[i, Integrate]},
TransformationFunctions -> {Automatic,
# /. i : Inactive[Integrate][f_, {_, _, _}] :> factorConstants[i] &,
# /. changeVar /@
DeleteDuplicates@
Cases[#, Inactive[Integrate][f_, {v_, a_, b_}], Infinity] &},
ComplexityFunction -> (LeafCount[#] +
10 Count[{#},
Inactive[Integrate][Inactive[Integrate][__], __], Infinity] +
5 Count[{#},
Inactive[Integrate][f_, __] /; ! FreeQ[f, Integrate], Infinity] -
Count[{#}, Power[Inactive[Integrate][f_, __], _], Infinity] +
Length@DeleteDuplicates@
Cases[{#}, Inactive[Integrate][f_, {x_, _, _}] :> x, Infinity] &)
] /. changeVar[v]) /; FreeQ[v, Missing]
];
someFunction[expr_] := expr;
The first two terms after LeafCount of the ComplexityFunction penalize nested integrals, the first one penalizing nested integrals in which the constants have not been factored out. The next term, being subtracted, rewards gathering powers of integrals. The last term penalizes extra variables of integration, which drives Simplify to prefer changing variables to the same thing.
