# Simplify for precision

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.

• I think the details are actually relevant here. The expression you give can be written as a matrix-vector contraction {i, j} . {{2, 5}, {3, 7}} . {x, y}. As usual in such expressions, the order of contraction matters. Contracting from the left allows staying in the integers longer (your ansC); contracting from the right means switching to floats right away (your ansB). So if your expression in general can be written in this form (or similar), then arguing through contraction-order will help. Commented Jul 18, 2021 at 6:36

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.

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

• This could be excellent. I will test and try to make it recognise addition and subtraction. Maybe adding a LeafCount÷1000000. Commented Jul 18, 2021 at 16:01
• @jdaw1 Please note that subtraction is represented as addition internally: for example, FullForm[a - b] shows an internal representation of Plus[a, Times[-1, b]], which is $a + (-1)\cdot b$. So there is no need to treat subtraction differently from what I've presented. Commented Jul 18, 2021 at 18:54
• Please, why Orderless? Consider 7+i+x, which is one floating operation, but x+i+7 is two. They aren’t the same. Commented Jul 18, 2021 at 22:03
• Orderless because of commutativity. As for multi-argument operations like $7+i+x$, what I've proposed is a bare-minimum simplistic implementation that you are free to extend in the directions you suggest, if that's what gets you closer to a solid heuristic for precision. Commented Jul 19, 2021 at 7:16
• @yarchik I think I need to convert to a form like yours, but strictly binary (which I’ll name, PostScript-style, add and mul), and then tell Mathematica that the binary add and mul are each commutative and associative, and jointly distributive (please, clues about how to this would be welcomed). That would allow Mathematica to hold differently 7 i add f add (better) and 7 i f add add` (worse) (written in PostScript’s reverse Polish notation). Then your complexity function can guide M’s optimisation in the wanted direction. This is a plan. Commented Jul 19, 2021 at 19:54