How to reduce the number of operations needed to evaluate an expression? [duplicate]

Question

For some numerical calculations in C++ I have, for example, this complicated expression to evaluate. Given as Mathematica input, in its original form it is

w - 4(w - y)((w - y)^2 y + 6(1 + y)((w - y)y + (1 + y)^2)) /
((w - y)^2 ((w - y)y + 6y^2 + 8y(1 + y)) + (1 + y)^2 (36(w - y)y + 24(1 + y)^2))

The goal is to reduce the number of required floating point operations to evaluate it in C++. Using useful code from Counting multiplications (complexity function) we can see it requires the following operations and their respective counts: {{Times, 22}, {Plus, 18}, {Power, 7}}. The FullSimplify'ed expression is a bit shorter

w - 4(w - y)(6 + y(18 + 6w + w^2 + 4(3 + w)y + y^2)) /
(24 + y(w^3 + w^2 (8 + 11y) + (2 + y)(48 + y(30 + y)) + w(36 + y(56 + 11y))))

and requires the following operation counts: {{Times, 16}, {Plus, 18}, {Power, 5}}.

One way to reduce the number of operations is to identify common components and evaluate them beforehand as temporaries. After staring at the expression for an hour I could find the following substitutions:

f0 = w - y;
f1 = 1 + y;
f00 = f0 f0;
f11 = f1 f1;
f0y = f0 y;

and the expression now reads

w - 4f0(6f1(f11 + f0y) + f00 y) / (f11(24f11 + 36f0y) + f00(6y y + 8f1 y + f0y))

The operation count is now: {{Times, 19}, {Plus, 9}, {Power, 0}}.

Is there a way to express this process of finding sub-expressions more formally in Mathematica? And automatize it at least for such simple expressions involving only multiplication and addition?

Answer 1

You can use the expression optimizer that is used within compile:

in = w - 4 (w - 
     y) ((w - y)^2 y + 
      6 (1 + y) ((w - y) y + (1 + y)^2))/((w - y)^2 ((w - y) y + 
         6 y^2 + 8 y (1 + y)) + (1 + y)^2 (36 (w - y) y + 
         24 (1 + y)^2));

Experimental`OptimizeExpression[in]
(*
Experimental`OptimizedExpression[
 Block[{Compile`$1, Compile`$2, Compile`$5, Compile`$3, Compile`$6, 
       Compile`$7}, Compile`$1 = -y; Compile`$2 = w + Compile`$1; 
      Compile`$5 = 1 + y; Compile`$3 = Compile`$2^2; 
      Compile`$6 = Compile`$2 y; Compile`$7 = Compile`$5^2; 
      w - (4 Compile`$2 (Compile`$3 y + 
        6 Compile`$5 (Compile`$6 + 
               Compile`$7)))/(Compile`$3 (Compile`$6 + 6 y^2 + 
             8 y Compile`$5) + 
      Compile`$7 (36 Compile`$2 y + 24 Compile`$7))]]
*)

hope this helps.

Answer 2

Maybe we can learn something from what Compile produces.

cf = Compile[
  {{w, _Real}, {y, _Real}}
  ,
  w - 4 (w - 
      y) ((w - y)^2 y + 
       6 (1 + y) ((w - y) y + (1 + y)^2))/((w - y)^2 ((w - y) y + 
          6 y^2 + 8 y (1 + y)) + (1 + y)^2 (36 (w - y) y + 
          24 (1 + y)^2))
  ]

In the compiled code we see that at least the number of squares was reduced. The option "ExpressionOptimization" is responsible for this. The default setting is Automatic, which if there are no external calls to the kernel, is basically equivalent to True. The compiled code can be viewed as follows

Needs["CompiledFunctionTools`"]
cf // CompilePrint

giving

        2 arguments
        6 Integer registers
        12 Real registers
        Underflow checking off
        Overflow checking off
        Integer overflow checking on
        RuntimeAttributes -> {}

        R0 = A1
        R1 = A2
        I0 = 1
        I2 = 6
        I3 = 8
        I5 = 24
        I1 = 4
        I4 = 36
        Result = R10


1   R2 = - R1
2   R3 = R0 + R2
3   R4 = I0
4   R4 = R4 + R1
5   R5 = Square[ R3]
6   R6 = R3 * R1
7   R7 = Square[ R4]
8   R8 = R5 * R1
9   R9 = R6 + R7
10  R10 = I2
11  R10 = R10 * R4 * R9
12  R8 = R8 + R10
13  R10 = Square[ R1]
14  R9 = I2
15  R9 = R9 * R10
16  R10 = I3
17  R10 = R10 * R1 * R4
18  R11 = R6 + R9 + R10
19  R9 = R5 * R11
20  R11 = I4
21  R11 = R11 * R3 * R1
22  R10 = I5
23  R10 = R10 * R7
24  R11 = R11 + R10
25  R10 = R7 * R11
26  R9 = R9 + R10
27  R10 = Reciprocal[ R9]
28  R8 = R8 * R10
29  R10 = I1
30  R10 = R10 * R3 * R8
31  R8 = - R10
32  R10 = R0 + R8
33  Return