Distributing PlusMinus consistently

I would like to give the built-in function PlusMinus, which is undefined by default, the intuitive meaning: considering all possible combinations of $+$ and $-$.

For simple expressions, this can be implemented as follows:

PlusMinus[a_,b_] := Flatten[{a+b,a-b}]
a \[PlusMinus] b \[PlusMinus] c
{a+b+c, a-b+c, a+b-c, a-b-c}

However, there are of course more complicated examples, such as when $\pm$ appears in a function's argument, for example I would like to do the same thing for arbitrary (non-listable) functions $f(a\pm b\pm\cdots)$, which should by the same principle yield

{f[a+b+...], f[a-b+...], ...}

To complicate things even further, consider

• $f(a\pm b, c)$
• $a+f(b+c\pm d,e)$
• $\frac{-b\pm\sqrt{b^2-4ac}}{2a}$
• ...

Can this be implemented purely by giving PlusMinus the appropriate definitions?

• However be aware that this is not the meaning most commonly assigned to it. Rather, in most cases $a\pm b\pm c$ is considered to be either $a+b+c$ or $a-b-c$. That's why there's a corresponding $\mp$ (MinusPlus in Mathematica): $a\pm b\mp c$ is either $a+b-c$ or $a-b+c$). – celtschk Jun 18 '12 at 16:36
• @celtschk I think your definition (the standard one) is more useful, but programming it seems more complicated, since you have to keep track of all the appearances of $\pm$ in an expression and change them consistently. I have no idea how to do that. – becko Feb 25 '16 at 10:06

Here's my shot, haven't tested it much yet, but it's pretty weird. Flattens out in NumericFunctions, I thought that made sense.

ClearAll[PlusMinus];
Module[{PlusMinusList},
SetAttributes[PlusMinus, {Flat, OneIdentity, NumericFunction}];
SetAttributes[PlusMinusList, {Flat, OneIdentity}];
PlusMinus[a_] := PlusMinusList[a, -a];
PlusMinus[a_, b_] := PlusMinusList[a + b, a - b];
PlusMinusList /:
h_Symbol?(MemberQ[Attributes[#], NumericFunction] &)[b___,
pm_PlusMinusList, a___] :=
Block[{PlusMinusList}, h[b, #, a] & /@ pm];
PlusMinusList[exp___] := {exp} /; Length@Stack[] === 4;
PlusMinusList /:
h_?(Head[#] =!= Symbol || !
MemberQ[Attributes[#], NumericFunction] &)[bef___,
PlusMinusList[pm___], aft___] :=
h[bef, {pm}, aft];
]

The idea is that it splits the results not in a regular List but in a PlusMinusList, which flattens itself out. With UpValues, it distributes over NumericFunctions. Then I put a couple of weird definitions to turn the PlusMinusList into a List in two cases: when it's already in the highest level of the stack, and when it is wrapped up by a non numeric function. But I'm already seing that it doesn't Flatten properly nested PlusMinus because it doens't have the NumericFunction Attribute

EDIT

I added the attribute NumericFunction to PlusMinus, and removed the condition "/; h =!= PlusMinusList;" which I think serves no purpose. Also added a line to consider the single argument case PlusMinus[x]

• $a+f(b\pm c)$ will result in a + f[{b + c, b - c}] :-( – David Jan 26 '12 at 17:07
• @David, That's because it is set to distribute on NumericFunctions. I thought it was a good idea, otherwise it would end up doing big messes since you coundn't ever treat it as a list. Just before doign that do SetAttributes[f, NumericFunction] so MMA knows that it's a function that's supposed to return a numeric value if given numeric inputs – Rojo Jan 26 '12 at 17:20
• Ah, now it works, thank you very much! – David Jan 26 '12 at 17:26

I managed to throw together sort of a solution. Note however that this is not what I am looking for (since I'm defining a helper function, whereas I would like the symbol itself to behave like this), the question remains open. However, I thought this might be useful for others nevertheless.

The algorithm works by searching all occurrences of PlusMinus, and replaces them with Slot. The resulting expression is then converted to a Function, and this function is applied to all the tuples of Plus and Subtract up to the length of how many PlusMinus there were in the first place.

distributePlusMinus[expression_] := Module[{count, tuples, i = 0, functionBody},
count = Count[expression, PlusMinus, ∞, Heads -> True];
tuples = Tuples[{Plus, Subtract}, count];
functionBody = Replace[expression, PlusMinus :> Slot[++i], ∞, Heads -> True];
Function[Evaluate@functionBody] @@@ tuples
]

Test results:

$a\pm b:$

test1 = a ± b;
distributePlusMinus[test1]
{a + b, a - b}

$f(a\pm b):$

test2 = f[a ± b];
distributePlusMinus[test2]
{f[a + b], f[a - b]}

$a\pm f(b\pm c):$

test3 = a ± f[b ± c];
distributePlusMinus[test3]
{a + f[b + c], a + f[b - c], a - f[b + c], a - f[b - c]}

$a\pm b\pm f(c\pm d,e):$

test4 = a ± b ± f[c ± d, e];
distributePlusMinus[test4]
{a+b+f[c+d,e],a+b+f[c-d,e],a-b+f[c+d,e],a-b+f[c-d,e],a+b-f[c+d,e],a+b-f[c-d,e],a-b-f[c+d,e],a-b-f[c-d,e]}

Note: These results depend on the fact that Mathematica splits up $a\pm b\pm c$ to $(a\pm b)\pm c$. Therefore, the algorithm needs adjustments for functions that do not exhibit this behavior. Especially, entering PlusMinus[a,b,c] is left as it is, while expressions involving $\pm$ are paired up as mentioned. The whole thing therefore works for the pretty-printed notation only.