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    *)