Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

This question already has an answer here:

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?

share|improve this question
    
Surely you know that aggressive reordering of expressions can change (often decrease) rounding accuracy of IEEE754 end result? If your goal is to reduce amount of operations in code generated by C compiler, you can use standards-noncompliant compiler option such as --ffast-math which probably knows very well how to reorder this kind of expressions. Otherwise - well, I'm intrigued to hear of a solution in Mathematica. –  kirma Sep 11 '13 at 14:42
    
@kirma Yes, I know. But I don't believe compilers dare to do this kind of temporary substitution; IMHO the existing operations are just reordered in time in order to optimize usage of various FP units and reduce congestion of the pipelines. Nevertheless, you raise a valid point, but the way how this expression was obtained is by deriving it first on paper and then some additional work was done with Mathematica and Simplify. You can see that during all these stages, the expressions get reordered randomly, either by human standards of beauty or length. –  Darko Veberic Sep 11 '13 at 15:02
    
@kirma This could be an interesting question, what is the form of this expression so that the rounding error is minimal? –  Darko Veberic Sep 11 '13 at 15:04
    
Compilers such as gcc have --ffast-math specifically for this purpose. It's enabled only on "unsafe" optimization levels by default if ever. The funniest part about IEEE754 arithmetic is that even summation is an ordered operation; a+b+c is likely to produce different result from a+c+b. Compilers wear the consistency straightjacket by default, but if you use --ffast-math, are relieved (at least to some useful extent) of it. Of course, this is relatively irrelevant point if your question is considered purely on the context of Mathematica. –  kirma Sep 11 '13 at 15:10
add comment

marked as duplicate by rm -rf Sep 11 '13 at 16:39

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

2 Answers

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.

share|improve this answer
    
Excellent. Remarkably, your output is really identical to the temporaries I found above. Thanks! –  Darko Veberic Sep 11 '13 at 14:45
    
I guess this output could be with some pattern matching simply transformed into a nested With[] statement... –  Darko Veberic Sep 11 '13 at 15:15
add comment

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
share|improve this answer
    
Thanks! This is a bit more obfuscated, but usefully more verbose. –  Darko Veberic Sep 11 '13 at 14:45
add comment

Not the answer you're looking for? Browse other questions tagged or ask your own question.