Simplify for precision

Question

I have an a set of equations to solve — details irrelevant — the result of which will be converted into PostScript, a 1980s printer-control language the only floating precision of which is single. Therefore I care about precision.

Consider a simple case, in which there are floating-point variables (or parameters) x and y, and small integer parameters i and j (small in the sense that any relevant combination of i and j is below integer overflow). Consider the following:

(* Seven floating-point operations, four multiplies, three additions *)
ansA = 2 x i + 3 x j + 5 y i + 7 y j;

(* Seven floating-point operations, four multiplies, three additions *)
ansB = (2 x + 5 y) i + (3 x + 7 y) j;

(* Three floating-point operations, two multiplies, one addition *)
ansC = (2 i + 3 j) x + (5 i + 7 j) y;

And[ansA == ansB, ansA == ansC, ansB == ansC] // Simplify

Clearly, precision arguments favour the last of these. Please, how can Mathematica’s Simplify (or FullSimplify) best be commanded to find the precision-optimal expression? If relevant, I have four integers parameters, and three floating-point parameters.

Answer 1

I don't think MA can do this kind of optimization for such small expressions. But you can play with Compile in order to see how it is evaluated:

cf = Compile[{{x, _Real}, {y, _Real}, {i, _Integer}, {j, _Integer}}, 
             Evaluate@ansA, 
             "RuntimeOptions" -> "Speed", 
             CompilationOptions -> {"ExpressionOptimization" -> True}]

You can printed generated code:

Needs["CompiledFunctionTools`"]
CompilePrint[cf]

    4 arguments
    6 Integer registers
    7 Real registers
    Underflow checking off
    Overflow checking off
    Integer overflow checking off
    RuntimeAttributes -> {}

    R0 = A1
    R1 = A2
    I0 = A3
    I1 = A4
    I4 = 5
    I2 = 2
    I5 = 7
    I3 = 3
    Result = R2

1   R2 = I2
2   R3 = I0
3   R2 = R2 * R3 * R0
4   R3 = I3
5   R4 = I1
6   R3 = R3 * R4 * R0
7   R4 = I4
8   R5 = I0
9   R4 = R4 * R5 * R1
10  R5 = I5
11  R6 = I1
12  R5 = R5 * R6 * R1
13  R2 = R2 + R3 + R4 + R5
14  Return

Notice that in more complicated cases this approach works.

Answer 2

Here's an approach that counts the involved operations (addition, multiplication) and remembers which ones of them involve floats and which ones involve only integers. This is of course a simplistic model of precision or complexity; but maybe it can be of use for you.

SetAttributes[{myplus, mytimes}, Orderless];
operations[A_] := Module[{B, R},
  (* substitute integers by "i" and floats by "f" *)
  B = A /. {Plus -> myplus, Times -> mytimes} /.
           {i | j -> "i", x | y -> "f", _Integer -> "i", _?NumericQ -> "f"};
  (* execute the formula and remember operations *)
  R = Reap[B //. {myplus[a : ("i"|"f")..] :>
                    (Sow[p[a]];If[Union[{a}]==={"i"},"i","f"]),
                  mytimes[a : ("i"|"f")..] :>
                    (Sow[t[a]];If[Union[{a}]==={"i"},"i","f"])}][[2, 1]]];

Notice how in the definition of R the floating-point numbers are "contagious" in that a sum or product involving at least one float results in a float.

As an example, we define the complexity of a formula as the number of multiplications involving at least one floating-point number:

complexity[A_] := Count[operations[A], t["f", __]]

Let's try it out: the three given formulas contain different numbers of floating-point multiplications:

ansA = 2 x i + 3 x j + 5 y i + 7 y j;
ansB = (2 x + 5 y) i + (3 x + 7 y) j;
ansC = (2 i + 3 j) x + (5 i + 7 j) y;

complexity[ansA]
(*    4    *)

complexity[ansB]
(*    6    *)

complexity[ansC]
(*    2    *)

We can use this complexity function in FullSimplify to discover ansC automatically:

FullSimplify[ansA, ComplexityFunction -> complexity]
(*    (2 i + 3 j) x + (5 i + 7 j) y    *)

Different choices for complexity will give different results here, and some experimentation may be needed.

Update: more detailed complexity function

A more detailed complexity function would count the total number of binary operations that involve at least one floating-point number. For example, $x+y$ would be one such operation; $x+i$ would be one; $x+i+2=x+(i+2)$ would be one (because $i+2$ can be done in the integers); $x+y+2$ would be two floating-point operations.

complexity[A_] := 
  Total[operations[A] /. (t | p)[a : ("i" | "f") ..] :>
          Count[{a}, "f"] - Boole[FreeQ[{a}, "i"]]]

Try out this complexity function:

complexity[ansA]
(*    7    *)

complexity[ansB]
(*    9    *)

complexity[ansC]
(*    3    *)

FullSimplify[ansA, ComplexityFunction -> complexity]
(*    (2 i + 3 j) x + (5 i + 7 j) y    *)