a / (b * c) is more efficient than a / b / c

[syclik]

Summary:

See @bob-carpenter's comment

This appears in at least:

• stan/math/prim/scal/prob/neg_binomial_lccdf.hpp

Description:

This will make the math library a little faster.

Current Version:

v2.16.0

[seantalts]

On my machine with primitives, these takes exactly the same amount of time and we can do a billion of them in a second. Is this more important with autodiff variables for some reason?

[bob-carpenter]

Two things. (1) Yes, more important with autodiff. Follow the chain rule---the expression for the derivatives involves more operations.

(2) I used the following evaluation code to make sure the compiler doesn't optimize away anything. Multiplication is way faster. What surprised me is that addition took about the same time as multiplication. The loop arithmetic's all basically free since the CPU has a separate pipeline for integers, though there is some overhead from the ongoing additions.

#include <iostream>

int main() {
  double a;
  for (double x = 1; x < 100000; ++x)
    for (double y = 1; y < 100000; ++y)
      a += x / y;
  std::cout << "a = " << a << std::endl;
}

Division:

~/temp2$ clang++ -O3 race.cpp
~/temp2$ time ./a.out
a = 6.04501e+10

real	0m45.319s
user	0m45.218s
sys	0m0.069s

Multiplication:

~/temp2$ time ./a.out
a = 2.49995e+19

real	0m9.608s
user	0m9.593s
sys	0m0.010s

[bob-carpenter]

I didn't test a / b / c---that might get unfolded into a / (b * c). They're not quite equivalent w.r.t. overflow behavior at the size limits.

[seantalts]

I tested this as well with the following code:

#include <stdlib.h>

int main (int argc, char* argv[]) {
  double a = atof(argv[0]);
  double b = atof(argv[1]);
  double c = atof(argv[2]);
  double d = 0;
  for (int i = 0; i < 1000000000; ++i) {
    //d += a / b / c;
    d += a / (b * c);
  }
  return d;
}

and it was the same time to two decimal places with -O3. I also checked with a disassembler to make sure that there were two division instructions vs one division and one multiplication. My theory here is that you can actually pipeline two divisions at once for the same cost as one division and one multiplication? Pipelining throws a lot of microbenchmark results out the window in non-obvious ways.

I'll work out the chain rule stuff on paper when I get home tonight.

[bbbales2]

An issue with this ``` for (int i = 0; i < 1000000000; ++i) { //d += a / b / c; d += a / (b * c); } ``` is that the compiler could recognize a / (b * c) is a constant expression across the loop and pull it out. I think the clang++ on my computer was doing this.

…

On 8/10/17, seantalts ***@***.***> wrote: I tested this as well with the following code: ```cpp #include <stdlib.h> int main (int argc, char* argv[]) { double a = atof(argv[0]); double b = atof(argv[1]); double c = atof(argv[2]); double d = 0; for (int i = 0; i < 1000000000; ++i) { //d += a / b / c; d += a / (b * c); } return d; } ``` and it was the same time to two decimal places with -O3. I also checked with a disassembler to make sure that there were two division instructions vs one division and one multiplication. My theory here is that you can actually pipeline two divisions at once for the same cost as one division and one multiplication? Pipelining throws a lot of microbenchmark results out the window in non-obvious ways. I'll work out the chain rule stuff on paper when I get home tonight. -- You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub: #596 (comment)

[bob-carpenter]

@bbbales2 hit the nail on the head here. You'll see that I used the loop variables and saved the result for printing at the end so the whole thing couldn't be compiled away.

[seantalts]

How embarrassing, I must have been looking at the -O0 disassembly - the O3 does in fact do loop-invariant code motion. I feel like so often we have to do little tricks that it seems like the compiler should be able to optimize (like the one for this ticket) that I'd lost some of the respect due to the compiler's optimizer.

I converted Bob's code to a test of the original question (two divisions vs one division and one multiplication) here:

#include <iostream>

int main() {
  double a = 0.0;
  for (double x = 1; x < 1000; ++x)
    for (double y = 1; y < 1000; ++y)
      for (double z = 1000; z > 1; --z)
        //a += x / y / z;
        a += x / (y * z);
  std::cout << "a = " << a << std::endl;
}

Do you guys see anything wrong with that?

It's also the same speed at -O3 for me. In the assembly, I actually see the proposed enhancement as having 3 additional instructions (all mulsds) and retaining the divisions, for some reason.

Here's the loop for double division:

0000000100000980         xorps      xmm4, xmm4                                  ; CODE XREF=_main+133
0000000100000983         cvtsi2sd   xmm4, edx
0000000100000987         movapd     xmm5, xmm3
000000010000098b         divsd      xmm5, xmm4
000000010000098f         addsd      xmm5, xmm7
0000000100000993         lea        esi, dword [rdx-1]
0000000100000996         xorps      xmm4, xmm4
0000000100000999         cvtsi2sd   xmm4, esi
000000010000099d         movapd     xmm6, xmm3
00000001000009a1         divsd      xmm6, xmm4
00000001000009a5         addsd      xmm6, xmm5
00000001000009a9         lea        esi, dword [rdx-2]
00000001000009ac         xorps      xmm4, xmm4
00000001000009af         cvtsi2sd   xmm4, esi
00000001000009b3         movapd     xmm7, xmm3
00000001000009b7         divsd      xmm7, xmm4
00000001000009bb         addsd      xmm7, xmm6
00000001000009bf         add        edx, 0xfffffffd
00000001000009c2         cmp        edx, 0x1
00000001000009c5         jg         loc_100000980

And here's the loop with the proposed enhancement (b*c):

0000000100000970         xorps      xmm3, xmm3                                  ; CODE XREF=_main+145
0000000100000973         cvtsi2sd   xmm3, edx
0000000100000977         mulsd      xmm3, xmm2
000000010000097b         movapd     xmm4, xmm1
000000010000097f         divsd      xmm4, xmm3
0000000100000983         addsd      xmm4, xmm6
0000000100000987         lea        esi, dword [rdx-1]
000000010000098a         xorps      xmm3, xmm3
000000010000098d         cvtsi2sd   xmm3, esi
0000000100000991         mulsd      xmm3, xmm2
0000000100000995         movapd     xmm5, xmm1
0000000100000999         divsd      xmm5, xmm3
000000010000099d         addsd      xmm5, xmm4
00000001000009a1         lea        esi, dword [rdx-2]
00000001000009a4         xorps      xmm3, xmm3
00000001000009a7         cvtsi2sd   xmm3, esi
00000001000009ab         mulsd      xmm3, xmm2
00000001000009af         movapd     xmm6, xmm1
00000001000009b3         divsd      xmm6, xmm3
00000001000009b7         addsd      xmm6, xmm5
00000001000009bb         add        edx, 0xfffffffd
00000001000009be         cmp        edx, 0x1
00000001000009c1         jg         loc_100000970

[seantalts]

So I think there's a little weirdness going on in there due to the compiler deciding loop variables should be integers and then converting them to doubles every iteration. So I just took my original test and added volatile keywords, which tell the compiler that the variables could change at any time and to basically avoid optimizing those variables:

#include <stdlib.h>

int main (int argc, char* argv[]) {
  volatile double a = atof(argv[0]);
  volatile double b = atof(argv[1]);
  volatile double c = atof(argv[2]);
  volatile double d = 0;
  for (int i = 0; i < 1000000000; ++i) {
    //d += a / b / c;
    d += a / (b * c);
  }
  return d;
}

I still see the same times observed for each version but I have simpler assembly for us to look at. Something weird is still going on in the (b*c) version as it again has 2 multiplications and 2 divisions for some reason.
Here's the a / b / c loop

0000000100000f60         movsd      xmm0, qword [rbp+var_28]                    ; CODE XREF=_main+107
0000000100000f65         divsd      xmm0, qword [rbp+var_20]
0000000100000f6a         divsd      xmm0, qword [rbp+var_18]
0000000100000f6f         addsd      xmm0, qword [rbp+var_10]
0000000100000f74         movsd      qword [rbp+var_10], xmm0
0000000100000f79         dec        eax
0000000100000f7b         jne        loc_100000f60

Here's the a / (b * c) loop:

0000000100000f40         movsd      xmm0, qword [rbp+var_28]                    ; CODE XREF=_main+141
0000000100000f45         movsd      xmm1, qword [rbp+var_20]
0000000100000f4a         mulsd      xmm1, qword [rbp+var_18]
0000000100000f4f         divsd      xmm0, xmm1
0000000100000f53         addsd      xmm0, qword [rbp+var_10]
0000000100000f58         movsd      qword [rbp+var_10], xmm0
0000000100000f5d         movsd      xmm0, qword [rbp+var_28]
0000000100000f62         movsd      xmm1, qword [rbp+var_20]
0000000100000f67         mulsd      xmm1, qword [rbp+var_18]
0000000100000f6c         divsd      xmm0, xmm1
0000000100000f70         addsd      xmm0, qword [rbp+var_10]
0000000100000f75         add        eax, 0xfffffffe
0000000100000f78         movsd      qword [rbp+var_10], xmm0
0000000100000f7d         jne        loc_100000f40

[bbbales2]

This code might be tricky to read cause I think it's unrolled a little. ``` 00000001000009bb add edx, 0xfffffffd ``` That is incrementing the loop by -3 maybe? At least that's what I saw when I looked at the assembly produced by Clang. Maybe just look at the assembly with -O1 to avoid the super duper optimizations. Also: ``` for (double z = 1000; z > 1; --z) a += x / y / z; ``` For this loop ```x / y``` is getting computed and cached outside of the z loop for me haha, which would make the results confusing still. Might be worth just taking the assembly and directly modifying it to test the ops you're curious about. What Bob is saying though should be true. Division is relatively hard (maybe 4-5x slower? -- this seems right: https://stackoverflow.com/questions/39095993/does-each-floating-point-operation-take-the-same-time/39135689#39135689) than multiplication and addition which usually cost about the same.

[bob-carpenter]

Neat. I'm surprised the code for the multiplication in the denominator is so much worse. It winds up doing just as many divisions, which I really don't understand.

[seantalts]

@bbbales2 I know that division is slower than multiplication, I would just be surprised that the compiler writers wouldn't know that, you know? It's a very easy peephole optimization, easier even than the loop-invariant code motion we saw above. That said, I need to reason through the autodiff situation... seems like it's probably still a useful optimization in that case. I just want to make sure that when we undertake work like this that we understand how valuable it is and when it's applicable.

[seantalts]

Also interesting - the (b *c) would seem to be less numerically stable (more likely to overflow)? I don't have a good grasp of that. And this points out that the way intermediate calculations are represented as doubles could affect the calculation differently, though it's not clear that matters to us.

[bbbales2]

I would just be surprised that the compiler writers wouldn't know that, you know? It's a very easy peephole optimization, easier even than the loop-invariant code motion we saw above.

Haha, that's a good point. Hadn't thought of it like that. Guess it all comes down to the autodiff? The compiler isn't gonna be able to do the optimizations with stan::var s, and the adjoints are going to be slightly more annoying with the division.