programming - 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?

Answer

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]