calculus and analysis - Solution to a specific problem caused by generic simplification

I'm trying to get MMA to help me evaluate certain integrals of trig functions. Here is an example: (The actual expressions I want to evaluate are more complicated than this one, but this illustrates the problem.)

Assuming[
    n ∈ Integers && m ∈ Integers,
    Integrate[Cos[n π x] Cos[m π x], {x, 0, 1}]
]
(* 0 *)

This answer is of course wrong, strictly speaking. The correct answer is $\frac{1}{2}\left(\delta_{n,m}+\delta_{n,-m}\right)$. As discussed here, for instance, MMA aims to produce generically correct results, and the special case $m=\pm n$ ends up being overlooked. I understand all that.

My question is whether anyone can suggest a straightforward workaround to get MMA to produce a more generally correct result for slightly more complex cases such as Integrate[Sin[k π x] Cos[n π x] Sin[m π x], {x,0,1}].

To clarify, I will add that there is no difficulty evaluating the special cases, if you know what they are. For instance:

Assuming[

    n ∈ Integers,
    Integrate[Cos[n π x] Cos[n π x], {x, 0, 1}]
]
(* 1/2 *)

So this is one of those vexing questions where it's easy to find the answer, once you know what it is. And of course you can use trig identities to get the answer, but getting the signs right is a tedious, fiddly business. I'd like to let MMA do it for me.

In case anyone cares, these integrals arise from PDEs when the solutions are represented as a cosine series.

Thanks.

Answer

Here is a fully automated solution, using ideas from my previous answer and Dr. Hintze's. It is in the form of two functions: trigIntegralToPiecewise and piecewiseToDelta. trigIntegralToPiecewise takes the thing to be integrated (as a pure function of the integration variable) and a list of the integer variables. It returns a Function that evaluates to a Piecewise function of the variables. If trigIntegralToPiecewise returns successfully, this piecewise function should always be equivalent to the integral.

The function can then be used as input to piecewiseToDelta. This is a bit of a kludge but seems to work fairly well. It is called using something like:

{df, check1, check2} = piecewiseToDelta[pwf, {k, m, n}];

df1 is now a pure function of k, m, n that evaluates to a sum of Kronecker deltas. (This is useful not only because the delta form is typically more compact and easier to understand than the Piecewise form, but also because MMA is good at simplifying sums over Kronecker deltas.) It is not guaranteed to be equivalent to pwf. The check1 and check2 returns allow you to reassure yourself. check1 is an Inactive expression using Reduce that, activated, should evaluate to True if the two functions are equivalent. I have not always had good results with Reduce, so check2 is something else. check2[10] evaluates pwf[k,m,n]==df[k,m,n] for $k$, $m$, and $n$ from -10 to 10 and returns True if all equalities hold.

trigIntegralToPiecewise::nonzero = 
  "works only for generically zero integrals";

trigIntegralToPiecewise[
  func_,
  vars_

  ] := Module[{ri, num, den, dvs, cases, res, pw, pwfunc},
  ri = Integrate[func[x], {x, 0, 1}];
  num = Numerator[Together[ri]];
  den = Denominator[Together[ri]];
  If[0 =!= Simplify[num, Assumptions -> vars \[Element] Integers],
   Message[trigIntegralToPiecewise::nonzero];
   Return[$Failed]
   ];
  cases = BooleanConvert[Reduce[0 == den, vars]];
  cases = If[Head[cases] === Or,

List @@ cases,
{cases}
];
  cases = Reverse[Subsets[cases]];
  cases = And @@@ cases;
  cases = DeleteDuplicates[Reduce[#, vars] & /@ cases];
  pw = Table[
res = Simplify[
  Integrate[func[x], {x, 0, 1},
   Assumptions -> case],

  vars \[Element] Integers && case
  ];
{res, case},
{case, cases}
];
  pw = Append[pw, {0, True}]; (* To be safe *)
  pwfunc = Function @@ {
 vars,
 Piecewise[
  pw

  ]
 };
  pwfunc
  ]

piecewiseToDelta[
  pwf_,
  vars_,
  max_: 3,
  extrabasis_: {}

  ] := Module[{nv, basis, kds, iterators, lb, ub, cs, eqs, cvals, 
   dfunc, check1, check2},
  nv = Length[vars];
  iterators = 
   Transpose[{vars, ConstantArray[lb, nv], ConstantArray[ub, nv]}];
  basis = Join[{1}, vars, extrabasis];
  kds = Plus @@@ Tuples[Outer[Times, Drop[vars, 1], {1, -1}]];
  kds = Sum[
(Sum[C[i, j] basis[[j]], {j, Length[basis]}]) KroneckerDelta[
  vars[[1]], kds[[i]]],

{i, Length[kds]}
];
  cs = Flatten@Table[
 C[i, j],
 {j, Length[basis]}, {i, Length[kds]}
 ];
  eqs = Table @@ 
Prepend[iterators /. {lb -> -max, ub -> max}, 
 kds == pwf @@ vars];
  eqs = DeleteCases[Flatten[eqs], True];

  cvals = Solve[eqs, cs][[1]];
  dfunc = Function @@ {
 vars,
 kds /. cvals /. C[_, _] -> 0
 };
  check1 = Inactivate[
Assuming[
 vars \[Element] Integers,
 Simplify[
  Reduce[pwf @@ vars == dfunc @@ vars]

  ]
 ]
];
  check2[cmax_, it_: iterators] := Module[{iters},
iters = it /. {lb -> -cmax, ub -> cmax};
And @@ Flatten[
  Table @@ Prepend[
    iters,
    pwf @@ vars == dfunc @@ vars
    ]

  ]
];
  {dfunc, check1, check2}
  ]

Example:

pwf1 = trigIntegralToPiecewise[x ↦ Sin[k 𝝿 x] Cos[n 𝝿 x] Sin[m 𝝿 x], {k, m, n}];

{df1, check11, check21} = piecewiseToDelta[pwf1, {k, m, n}];
df1 // TraditionalForm

$$\{k,m,n\}↦-\frac{1}{4} \delta _{k,n-m}-\frac{1}{4} \delta _{k,-m-n}+\frac{\delta _{k,m-n}}{4}+\frac{\delta _{k,m+n}}{4}$$

check11[20]
(* True *)

Activate@check21
(* True *)