How much Algebra does C2 Know? Part 1: Associativity

Making loops execute faster is firmly rooted in algebra, but how much does C2 know or care about? When building a highly optimised query engine, a critical concern is the quality of assembly code generated for loops. There is a lot more to JIT compilation than loop optimisation; inlining, class hierarchy analysis, escape analysis to name but a few. Moreover, everything it does has to be fast since it shares resources with the application itself; it can’t spend time unless it brings a net benefit. Being such a generalist, does C2, the JIT compiler used in server applications, know high school algebra?

Specific knowledge of maths is not always worthwhile to program execution, even when it leads to high performance gains. As a motivating example, there is no way to refer directly to the natural numbers, for which there is a body of knowledge, in any programming language I have ever used. For instance, the sum of the first n natural numbers is \frac{(n+1)n}{2}, and most high school students know it. This expression is intuitively much faster to evaluate than the equivalent algorithm:

int sum(int n) {
    int total = 0;
    for (int i = 0; i <= n; ++i) {
        total += i;
    }
    return total;
}

But would this loop rewrite be a worthwhile optimisation? The expression takes about 3.5ns to compute the sum of the first million natural numbers, whereas the loop takes 350µs, so we can conclude that C2 does not know this formula and prefers brute force. I would be aghast if time had been spent on optimisations like this: unless your application spends a lot of time adding up contiguous ranges of natural numbers, the marginal benefit is negligible. If this is what your application does most, you should do it yourself. The possibility of an optimisation doesn’t imply its viability: there needs to be a benefit when considering engineering effort, speed improvement, reliability and ubiquity. While this optimisation fails miserably on the grounds of ubiquity, there’s useful schoolboy maths that C2 does seem to know.

Associativity and Dependencies

Each x86 instruction has a throughput – the number of cycles it takes to complete – and a latency – the number of cycles it takes before the result is available to the next instruction in a chain. These numbers are produced by processor vendors, but there are independent numbers like these from Agner Fog, the introduction of which goes into more detail defining terms like latency. At first, the latency number feels a bit like a scam: what use is an advertised throughput if we can’t use the result immediately? This is where pipelining comes in: independent instructions can be interleaved. If a loop operation is associative and there are no dependencies between iterations, then it can be unrolled, which enables pipelining. If a loop operation is also commutative, then out of order execution is permitted. Evidence of an unrolled loop suggests that the compiler has realised that an operation is at least associative.

To see this in action it’s necessary to find an associative loop reduction that the compiler can’t vectorise. I took an example from the RoaringBitmap library – computing the cardinality of a bitmap container – which is a perfect example to capture this behaviour, because bit counts cannot be vectorised in Java.

/**
   * Recomputes the cardinality of the bitmap.
   */
  protected void computeCardinality() {
    this.cardinality = 0;
    for (int k = 0; k < this.bitmap.length; k++) {
      this.cardinality += Long.bitCount(this.bitmap[k]);
    }
  }

we can see evidence of loop unrolling and out of order execution when looking at the assembly code emitted. The popcnt instructions are executed on the array out of order, and do not wait for the addition to the accumulator.

  0x0000018287aeef40: popcnt  r9,qword ptr [rbx+r13*8+10h]
                                                ;...f3
                                                ;...4e
                                                ;...0f
                                                ;...b8
                                                ;...4c
                                                ;...eb
                                                ;...10

  0x0000018287aeef47: movsxd  r8,r13d           ;...4d
                                                ;...63
                                                ;...c5

  0x0000018287aeef4a: popcnt  r10,qword ptr [rbx+r8*8+28h]
                                                ;...f3
                                                ;...4e
                                                ;...0f
                                                ;...b8
                                                ;...54
                                                ;...c3
                                                ;...28

  0x0000018287aeef51: popcnt  r11,qword ptr [rbx+r8*8+18h]
                                                ;...f3
                                                ;...4e
                                                ;...0f
                                                ;...b8
                                                ;...5c
                                                ;...c3
                                                ;...18

  0x0000018287aeef58: popcnt  rdx,qword ptr [rbx+r8*8+20h]
                                                ;...f3
                                                ;...4a
                                                ;...0f
                                                ;...b8
                                                ;...54
                                                ;...c3
                                                ;...20

  0x0000018287aeef5f: movsxd  r8,r9d            ;...4d
                                                ;...63
                                                ;...c1

  0x0000018287aeef62: add     r8,rbp            ;...4c
                                                ;...03
                                                ;...c5

  0x0000018287aeef65: movsxd  r9,edx            ;...4c
                                                ;...63
                                                ;...ca

To generate this assembly code you can run the project at github with the arguments

--include .*popcnt.* 
--print-assembly

The compiler does a very good job in this case: you can try unrolling the loop yourself, but you can only match performance if you guess the loop stride correctly. It’s impossible to prove a negative proposition, but it’s likely you’ll only make it worse if you try. C2 graduates with flying colours here: it definitely understands associativity and dependence.

The catch with pipelining is that an instruction must always wait for its operands. While the operation is associative, there is no way to reorder the code below.

    private int[] prefixSum(int[] data) {
        int[] result = new int[data.length];
        for (int i = 1; i < result.length; ++i) {
            result[i] = result[i - 1] + data[i];
        }
        return result;
    }

What happens with a prefix sum? There’s no unrolling: you can see the loop control statements have not been removed (look for commands like cmp ebx, inc ebx). The loop is also executed in order because it is sequentially dependent.

com/openkappa/simd/prefix/PrefixSum.impl([I)[I  [0x000001c21215bb20, 0x000001c21215bd78]  600 bytes
Argument 0 is unknown.RIP: 0x1c21215bb20 Code size: 0x00000258
[Entry Point]
[Verified Entry Point]
[Constants]
  # {method} {0x000001c2244379d0} 'impl' '([I)[I' in 'com/openkappa/simd/prefix/PrefixSum'
  0x000001c21215bb20: int3                      ;...cc

  0x000001c21215bb21: nop     word ptr [rax+rax+0h]  ;...66
                                                ;...66
                                                ;...66
                                                ;...0f
                                                ;...1f
                                                ;...84
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00

  0x000001c21215bb2c: nop                       ;...66
                                                ;...66
                                                ;...66
                                                ;...90

  0x000001c21215bb30: mov     dword ptr [rsp+0ffffffffffff9000h],eax
                                                ;...89
                                                ;...84
                                                ;...24
                                                ;...00
                                                ;...90
                                                ;...ff
                                                ;...ff

  0x000001c21215bb37: push    rbp               ;...55

  0x000001c21215bb38: sub     rsp,60h           ;...48
                                                ;...83
                                                ;...ec
                                                ;...60

  0x000001c21215bb3c: mov     rbp,qword ptr [rdx+10h]  ;...48
                                                ;...8b
                                                ;...6a
                                                ;...10

  0x000001c21215bb40: mov     r13,qword ptr [rdx+8h]  ;...4c
                                                ;...8b
                                                ;...6a
                                                ;...08

  0x000001c21215bb44: mov     ebx,dword ptr [rdx]  ;...8b
                                                ;...1a

  0x000001c21215bb46: mov     rcx,rdx           ;...48
                                                ;...8b
                                                ;...ca

  0x000001c21215bb49: mov     r10,51da8d20h     ;...49
                                                ;...ba
                                                ;...20
                                                ;...8d
                                                ;...da
                                                ;...51
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00

  0x000001c21215bb53: call indirect r10         ;...41
                                                ;...ff
                                                ;...d2

  0x000001c21215bb56: test    rbp,rbp           ;...48
                                                ;...85
                                                ;...ed

  0x000001c21215bb59: je      1c21215bcb9h      ;...0f
                                                ;...84
                                                ;...5a
                                                ;...01
                                                ;...00
                                                ;...00

  0x000001c21215bb5f: mov     r10d,dword ptr [rbp+8h]  ;...44
                                                ;...8b
                                                ;...55
                                                ;...08

  0x000001c21215bb63: cmp     r10d,0f800016dh   ;...41
                                                ;...81
                                                ;...fa
                                                ;...6d
                                                ;...01
                                                ;...00
                                                ;...f8
                                                ;   {metadata({type array int})}
  0x000001c21215bb6a: jne     1c21215bd1dh      ;...0f
                                                ;...85
                                                ;...ad
                                                ;...01
                                                ;...00
                                                ;...00

  0x000001c21215bb70: mov     r8,rbp            ;...4c
                                                ;...8b
                                                ;...c5

  0x000001c21215bb73: mov     r10d,dword ptr [r13+8h]  ;...45
                                                ;...8b
                                                ;...55
                                                ;...08
                                                ; implicit exception: dispatches to 0x000001c21215bd41
  0x000001c21215bb77: cmp     r10d,0f800016dh   ;...41
                                                ;...81
                                                ;...fa
                                                ;...6d
                                                ;...01
                                                ;...00
                                                ;...f8
                                                ;   {metadata({type array int})}
  0x000001c21215bb7e: jne     1c21215bd1dh      ;...0f
                                                ;...85
                                                ;...99
                                                ;...01
                                                ;...00
                                                ;...00
                                                ;*iload_3 {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@7 (line 21)

  0x000001c21215bb84: mov     edi,dword ptr [r13+0ch]  ;...41
                                                ;...8b
                                                ;...7d
                                                ;...0c
                                                ;*arraylength {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@9 (line 21)

  0x000001c21215bb88: cmp     ebx,edi           ;...3b
                                                ;...df

  0x000001c21215bb8a: jnl     1c21215bcaah      ;...0f
                                                ;...8d
                                                ;...1a
                                                ;...01
                                                ;...00
                                                ;...00
                                                ;*if_icmpge {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@10 (line 21)

  0x000001c21215bb90: mov     r11d,ebx          ;...44
                                                ;...8b
                                                ;...db

  0x000001c21215bb93: inc     r11d              ;...41
                                                ;...ff
                                                ;...c3

  0x000001c21215bb96: xor     r9d,r9d           ;...45
                                                ;...33
                                                ;...c9

  0x000001c21215bb99: mov     ecx,1h            ;...b9
                                                ;...01
                                                ;...00
                                                ;...00
                                                ;...00

  0x000001c21215bb9e: cmp     r11d,ecx          ;...44
                                                ;...3b
                                                ;...d9

  0x000001c21215bba1: cmovl   r11d,ecx          ;...44
                                                ;...0f
                                                ;...4c
                                                ;...d9

  0x000001c21215bba5: cmp     r11d,edi          ;...44
                                                ;...3b
                                                ;...df

  0x000001c21215bba8: cmovnle r11d,edi          ;...44
                                                ;...0f
                                                ;...4f
                                                ;...df

  0x000001c21215bbac: cmp     r11d,r9d          ;...45
                                                ;...3b
                                                ;...d9

  0x000001c21215bbaf: cmovl   r11d,r9d          ;...45
                                                ;...0f
                                                ;...4c
                                                ;...d9

  0x000001c21215bbb3: cmp     r11d,edi          ;...44
                                                ;...3b
                                                ;...df

  0x000001c21215bbb6: cmovnle r11d,edi          ;...44
                                                ;...0f
                                                ;...4f
                                                ;...df

  0x000001c21215bbba: nop                       ;...66
                                                ;...90
                                                ;*iastore {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@24 (line 22)

  0x000001c21215bbbc: mov     ebp,ebx           ;...8b
                                                ;...eb

  0x000001c21215bbbe: dec     ebp               ;...ff
                                                ;...cd
                                                ;*isub {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@18 (line 22)

  0x000001c21215bbc0: cmp     ebp,edi           ;...3b
                                                ;...ef

  0x000001c21215bbc2: jnb     1c21215bcc1h      ;...0f
                                                ;...83
                                                ;...f9
                                                ;...00
                                                ;...00
                                                ;...00

  0x000001c21215bbc8: mov     r9d,dword ptr [r8+0ch]  ;...45
                                                ;...8b
                                                ;...48
                                                ;...0c
                                                ;*iaload {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@22 (line 22)
                                                ; implicit exception: dispatches to 0x000001c21215bd31
  0x000001c21215bbcc: mov     ebp,dword ptr [r13+rbx*4+0ch]
                                                ;...41
                                                ;...8b
                                                ;...6c
                                                ;...9d
                                                ;...0c
                                                ;*iaload {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@19 (line 22)

  0x000001c21215bbd1: cmp     ebx,r9d           ;...41
                                                ;...3b
                                                ;...d9

  0x000001c21215bbd4: jnb     1c21215bce1h      ;...0f
                                                ;...83
                                                ;...07
                                                ;...01
                                                ;...00
                                                ;...00

  0x000001c21215bbda: add     ebp,dword ptr [r8+rbx*4+10h]
                                                ;...41
                                                ;...03
                                                ;...6c
                                                ;...98
                                                ;...10
                                                ;*iadd {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@23 (line 22)

  0x000001c21215bbdf: cmp     ebx,edi           ;...3b
                                                ;...df

  0x000001c21215bbe1: jnb     1c21215bd01h      ;...0f
                                                ;...83
                                                ;...1a
                                                ;...01
                                                ;...00
                                                ;...00

  0x000001c21215bbe7: mov     dword ptr [r13+rbx*4+10h],ebp
                                                ;...41
                                                ;...89
                                                ;...6c
                                                ;...9d
                                                ;...10
                                                ;*iastore {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@24 (line 22)

  0x000001c21215bbec: inc     ebx               ;...ff
                                                ;...c3
                                                ;*iinc {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@25 (line 21)

  0x000001c21215bbee: cmp     ebx,r11d          ;...41
                                                ;...3b
                                                ;...db

  0x000001c21215bbf1: jl      1c21215bbbch      ;...7c
                                                ;...c9
                                                ;*if_icmpge {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@10 (line 21)

  0x000001c21215bbf3: cmp     edi,r9d           ;...41
                                                ;...3b
                                                ;...f9

  0x000001c21215bbf6: mov     r11d,edi          ;...44
                                                ;...8b
                                                ;...df

  0x000001c21215bbf9: cmovnle r11d,r9d          ;...45
                                                ;...0f
                                                ;...4f
                                                ;...d9

  0x000001c21215bbfd: mov     r10d,r11d         ;...45
                                                ;...8b
                                                ;...d3

  0x000001c21215bc00: add     r10d,0fffffff9h   ;...41
                                                ;...83
                                                ;...c2
                                                ;...f9

  0x000001c21215bc04: mov     ecx,80000000h     ;...b9
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...80

  0x000001c21215bc09: cmp     r11d,r10d         ;...45
                                                ;...3b
                                                ;...da

  0x000001c21215bc0c: cmovl   r10d,ecx          ;...44
                                                ;...0f
                                                ;...4c
                                                ;...d1

  0x000001c21215bc10: cmp     ebx,r10d          ;...41
                                                ;...3b
                                                ;...da

  0x000001c21215bc13: jnl     1c21215bc80h      ;...7d
                                                ;...6b

  0x000001c21215bc15: nop     word ptr [rax+rax+0h]  ;...66
                                                ;...66
                                                ;...66
                                                ;...0f
                                                ;...1f
                                                ;...84
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;*isub {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@18 (line 22)

  0x000001c21215bc20: mov     ecx,dword ptr [r8+rbx*4+10h]
                                                ;...41
                                                ;...8b
                                                ;...4c
                                                ;...98
                                                ;...10

  0x000001c21215bc25: add     ecx,dword ptr [r13+rbx*4+0ch]
                                                ;...41
                                                ;...03
                                                ;...4c
                                                ;...9d
                                                ;...0c
                                                ;*iadd {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@23 (line 22)

  0x000001c21215bc2a: mov     dword ptr [r13+rbx*4+10h],ecx
                                                ;...41
                                                ;...89
                                                ;...4c
                                                ;...9d
                                                ;...10
                                                ;*iastore {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@24 (line 22)

  0x000001c21215bc2f: movsxd  r11,ebx           ;...4c
                                                ;...63
                                                ;...db

  0x000001c21215bc32: add     ecx,dword ptr [r8+r11*4+14h]
                                                ;...43
                                                ;...03
                                                ;...4c
                                                ;...98
                                                ;...14
                                                ;*iadd {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@23 (line 22)

  0x000001c21215bc37: mov     dword ptr [r13+r11*4+14h],ecx
                                                ;...43
                                                ;...89
                                                ;...4c
                                                ;...9d
                                                ;...14
                                                ;*iastore {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@24 (line 22)

  0x000001c21215bc3c: add     ecx,dword ptr [r8+r11*4+18h]
                                                ;...43
                                                ;...03
                                                ;...4c
                                                ;...98
                                                ;...18
                                                ;*iadd {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@23 (line 22)

  0x000001c21215bc41: mov     dword ptr [r13+r11*4+18h],ecx
                                                ;...43
                                                ;...89
                                                ;...4c
                                                ;...9d
                                                ;...18
                                                ;*iastore {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@24 (line 22)

  0x000001c21215bc46: add     ecx,dword ptr [r8+r11*4+1ch]
                                                ;...43
                                                ;...03
                                                ;...4c
                                                ;...98
                                                ;...1c
                                                ;*iadd {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@23 (line 22)

  0x000001c21215bc4b: mov     dword ptr [r13+r11*4+1ch],ecx
                                                ;...43
                                                ;...89
                                                ;...4c
                                                ;...9d
                                                ;...1c
                                                ;*iastore {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@24 (line 22)

  0x000001c21215bc50: add     ecx,dword ptr [r8+r11*4+20h]
                                                ;...43
                                                ;...03
                                                ;...4c
                                                ;...98
                                                ;...20
                                                ;*iadd {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@23 (line 22)

  0x000001c21215bc55: mov     dword ptr [r13+r11*4+20h],ecx
                                                ;...43
                                                ;...89
                                                ;...4c
                                                ;...9d
                                                ;...20
                                                ;*iastore {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@24 (line 22)

  0x000001c21215bc5a: add     ecx,dword ptr [r8+r11*4+24h]
                                                ;...43
                                                ;...03
                                                ;...4c
                                                ;...98
                                                ;...24
                                                ;*iadd {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@23 (line 22)

  0x000001c21215bc5f: mov     dword ptr [r13+r11*4+24h],ecx
                                                ;...43
                                                ;...89
                                                ;...4c
                                                ;...9d
                                                ;...24
                                                ;*iastore {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@24 (line 22)

  0x000001c21215bc64: add     ecx,dword ptr [r8+r11*4+28h]
                                                ;...43
                                                ;...03
                                                ;...4c
                                                ;...98
                                                ;...28
                                                ;*iadd {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@23 (line 22)

  0x000001c21215bc69: mov     dword ptr [r13+r11*4+28h],ecx
                                                ;...43
                                                ;...89
                                                ;...4c
                                                ;...9d
                                                ;...28

  0x000001c21215bc6e: add     ecx,dword ptr [r8+r11*4+2ch]
                                                ;...43
                                                ;...03
                                                ;...4c
                                                ;...98
                                                ;...2c

  0x000001c21215bc73: mov     dword ptr [r13+r11*4+2ch],ecx
                                                ;...43
                                                ;...89
                                                ;...4c
                                                ;...9d
                                                ;...2c
                                                ;*iastore {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@24 (line 22)

  0x000001c21215bc78: add     ebx,8h            ;...83
                                                ;...c3
                                                ;...08
                                                ;*iinc {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@25 (line 21)

  0x000001c21215bc7b: cmp     ebx,r10d          ;...41
                                                ;...3b
                                                ;...da

  0x000001c21215bc7e: jl      1c21215bc20h      ;...7c
                                                ;...a0
                                                ;*if_icmpge {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@10 (line 21)

  0x000001c21215bc80: cmp     ebx,edi           ;...3b
                                                ;...df

  0x000001c21215bc82: jnl     1c21215bcaah      ;...7d
                                                ;...26
                                                ;*iastore {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@24 (line 22)

  0x000001c21215bc84: mov     ebp,ebx           ;...8b
                                                ;...eb

  0x000001c21215bc86: dec     ebp               ;...ff
                                                ;...cd
                                                ;*isub {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@18 (line 22)

  0x000001c21215bc88: cmp     ebp,edi           ;...3b
                                                ;...ef

  0x000001c21215bc8a: jnb     1c21215bcc1h      ;...73
                                                ;...35

  0x000001c21215bc8c: mov     ebp,dword ptr [r13+rbx*4+0ch]
                                                ;...41
                                                ;...8b
                                                ;...6c
                                                ;...9d
                                                ;...0c
                                                ;*iaload {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@19 (line 22)

  0x000001c21215bc91: cmp     ebx,r9d           ;...41
                                                ;...3b
                                                ;...d9

  0x000001c21215bc94: jnb     1c21215bce1h      ;...73
                                                ;...4b

  0x000001c21215bc96: add     ebp,dword ptr [r8+rbx*4+10h]
                                                ;...41
                                                ;...03
                                                ;...6c
                                                ;...98
                                                ;...10
                                                ;*iadd {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@23 (line 22)

  0x000001c21215bc9b: cmp     ebx,edi           ;...3b
                                                ;...df

  0x000001c21215bc9d: jnb     1c21215bd01h      ;...73
                                                ;...62

  0x000001c21215bc9f: mov     dword ptr [r13+rbx*4+10h],ebp
                                                ;...41
                                                ;...89
                                                ;...6c
                                                ;...9d
                                                ;...10
                                                ;*iastore {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@24 (line 22)

  0x000001c21215bca4: inc     ebx               ;...ff
                                                ;...c3
                                                ;*iinc {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@25 (line 21)

  0x000001c21215bca6: cmp     ebx,edi           ;...3b
                                                ;...df

  0x000001c21215bca8: jl      1c21215bc84h      ;...7c
                                                ;...da
                                                ;*if_icmpge {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.prefix.PrefixSum::impl@10 (line 21)

  0x000001c21215bcaa: mov     rax,r13           ;...49
                                                ;...8b
                                                ;...c5

  0x000001c21215bcad: add     rsp,60h           ;...48
                                                ;...83
                                                ;...c4
                                                ;...60

  0x000001c21215bcb1: pop     rbp               ;...5d

  0x000001c21215bcb2: test    dword ptr [1c27b720000h],eax
                                                ;...85
                                                ;...05
                                                ;...48
                                                ;...43
                                                ;...5c
                                                ;...69
                                                ;   {poll_return}
  0x000001c21215bcb8: ret                       ;...c3

Does this harm performance? add takes 0.33 cycles, whereas popcnt takes 1 cycle per instruction. Shouldn’t a prefix sum be faster to calculate than a population count, on the same length of array and same width of integer? They can be compared head to head (implementing prefix sum for long[] to keep word width constant)

--include .*prefix.PrefixSum.PrefixSumLong|.*popcnt.PopCount.PopCount$

Far from having 3x throughput, the prefix sum is much worse. This is entirely because there is no loop unrolling and no pipelining. When possible, C2 applies aggressive unrolling optimisations unavailable to the programmer. For vectorisable operations (requiring linear independence and countability), loop unrolling further marks the loop as a candidate for auto-vectorisation.

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size
PopCount thrpt 1 10 9.174499 0.394487 ops/ms 100000
PopCount thrpt 1 10 1.217521 0.513734 ops/ms 1000000
PrefixSumLong thrpt 1 10 6.807279 0.925282 ops/ms 100000
PrefixSumLong thrpt 1 10 0.443974 0.053544 ops/ms 1000000

If the dependencies need to fetch data from RAM the latency can be much higher than loading from registers or from prefetched cache. Even when fetching from RAM, the worst case scenario, during this delay independent instructions can complete, unless they have a false dependency.

Zeroing Negative Values in Arrays Efficiently

Replacing negatives with zeroes in large arrays of values is a primitive function of several complex financial risk measures, including potential future exposure (PFE) and the liquidity coverage ratio (LCR). While this is not an interesting operation by any stretch of the imagination, it is useful and there is significant benefit in its performance. This is an operation that can be computed very efficiently using the instruction VMAXPD. For Intel Xeon processors, this instruction requires half a cycle to calculate and has a latency (how long before another instruction can use its result) of four cycles. There is currently no way to trick Java into using this instruction for this simple operation, though there is a placeholder implementation on the current DoubleVector prototype in Project Panama which may do so.

C++ Intel Intrinsics

It’s possible to target instructions from different processor vendors, in my case Intel, by using intrinsic functions which expose instructions as high level functions. The code looks incredibly ugly but it works. Here is a C++ function for 256 bit ymm registers:

void zero_negatives(const double* source, double* target, const size_t length) {
	for (size_t i = 0; i + 3 < length; i += 4) {
		__m256d vector = _mm256_load_pd(source + i);
		__m256d zeroed = _mm256_max_pd(vector, _mm256_setzero_pd());
		_mm256_storeu_pd(target + i, zeroed);
	}
}

The function loads doubles into 256 bit vectors, within each vector replaces the negative values with zero, and writes them back into an array. It generates the following assembly code (which, incidentally, is less of a shit show to access than in Java):

void zero_negatives(const double* source, double* target, const size_t length) {
00007FF746EE5110  mov         qword ptr [rsp+18h],r8  
00007FF746EE5115  mov         qword ptr [rsp+10h],rdx  
00007FF746EE511A  mov         qword ptr [rsp+8],rcx  
00007FF746EE511F  push        r13  
00007FF746EE5121  push        rbp  
00007FF746EE5122  push        rdi  
00007FF746EE5123  sub         rsp,250h  
00007FF746EE512A  mov         r13,rsp  
00007FF746EE512D  lea         rbp,[rsp+20h]  
00007FF746EE5132  and         rbp,0FFFFFFFFFFFFFFE0h  
00007FF746EE5136  mov         rdi,rsp  
00007FF746EE5139  mov         ecx,94h  
00007FF746EE513E  mov         eax,0CCCCCCCCh  
00007FF746EE5143  rep stos    dword ptr [rdi]  
00007FF746EE5145  mov         rcx,qword ptr [rsp+278h]  
	for (size_t i = 0; i + 3 < length; i += 4) {
00007FF746EE514D  mov         qword ptr [rbp+8],0  
00007FF746EE5155  jmp         zero_negatives+53h (07FF746EE5163h)  
00007FF746EE5157  mov         rax,qword ptr [rbp+8]  
00007FF746EE515B  add         rax,4  
00007FF746EE515F  mov         qword ptr [rbp+8],rax  
00007FF746EE5163  mov         rax,qword ptr [rbp+8]  
00007FF746EE5167  add         rax,3  
00007FF746EE516B  cmp         rax,qword ptr [length]  
00007FF746EE5172  jae         zero_negatives+0DDh (07FF746EE51EDh)  
		__m256d vector = _mm256_load_pd(source + i);
00007FF746EE5174  mov         rax,qword ptr   
00007FF746EE517B  mov         rcx,qword ptr [rbp+8]  
00007FF746EE517F  lea         rax,[rax+rcx*8]  
00007FF746EE5183  vmovupd     ymm0,ymmword ptr [rax]  
00007FF746EE5187  vmovupd     ymmword ptr [rbp+180h],ymm0  
00007FF746EE518F  vmovupd     ymm0,ymmword ptr [rbp+180h]  
00007FF746EE5197  vmovupd     ymmword ptr [rbp+40h],ymm0  
		__m256d zeroed = _mm256_max_pd(vector, _mm256_setzero_pd());
00007FF746EE519C  vxorpd      xmm0,xmm0,xmm0  
00007FF746EE51A0  vmovupd     ymmword ptr [rbp+200h],ymm0  
00007FF746EE51A8  vmovupd     ymm0,ymmword ptr [rbp+40h]  
00007FF746EE51AD  vmaxpd      ymm0,ymm0,ymmword ptr [rbp+200h]  
00007FF746EE51B5  vmovupd     ymmword ptr [rbp+1C0h],ymm0  
00007FF746EE51BD  vmovupd     ymm0,ymmword ptr [rbp+1C0h]  
00007FF746EE51C5  vmovupd     ymmword ptr [rbp+80h],ymm0  
		_mm256_storeu_pd(target + i, zeroed);
00007FF746EE51CD  mov         rax,qword ptr [target]  
00007FF746EE51D4  mov         rcx,qword ptr [rbp+8]  
00007FF746EE51D8  lea         rax,[rax+rcx*8]  
00007FF746EE51DC  vmovupd     ymm0,ymmword ptr [rbp+80h]  
00007FF746EE51E4  vmovupd     ymmword ptr [rax],ymm0  
	}
00007FF746EE51E8  jmp         zero_negatives+47h (07FF746EE5157h)  
}
00007FF746EE51ED  lea         rsp,[r13+250h]  
00007FF746EE51F4  pop         rdi  
00007FF746EE51F5  pop         rbp  
00007FF746EE51F6  pop         r13  
00007FF746EE51F8  ret    

This code is noticeably fast. I measured the throughput averaged over 1000 iterations, with an array of 10 million doubles (800MB) uniformly distributed between +/- 1E7, to quantify the throughput in GB/s and iterations/s. This code does between 4.5 and 5 iterations per second, which translates to processing approximately 4GB/s. This seems high, and since I am unaware of best practices in C++, if the measurement is flawed, I would gratefully be educated in the comments.

void benchmark() {
	const size_t length = 1E8;
	double* values = new double[length];
	fill_array(values, length);
	double* zeroed = new double[length];
	auto start = std::chrono::high_resolution_clock::now();
	int iterations = 1000;
	for (int i = 0; i < iterations; ++i) {
		zero_negatives(values, zeroed, length);
	}
	auto end = std::chrono::high_resolution_clock::now();
	auto nanos = std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count();
	double thrpt_s = (iterations * 1E9) / nanos;
	double thrpt_gbps = (thrpt_s * sizeof(double) * length) / 1E9;
	std::cout << thrpt_s << "/s" << std::endl;
	std::cout << thrpt_gbps << "GB/s" << std::endl;
	delete[] values;
	delete[] zeroed;
}

While I am sure there are various ways an expert could tweak this for performance, this code can’t get much faster unless there are 512 bit zmm registers available, in which case it would be wasteful. While the code looks virtually the same for AVX512 (just replace “256” with “512”), portability and efficiency are at odds. Handling the mess of detecting the best instruction set for the deployed architecture is the main reason for using Java in performance sensitive (but not critical) applications. But this is not the code the JVM generates.

Java Auto-Vectorisation (Play Your Cards Right)

There is currently no abstraction modelling vectorisation in Java. The only access available is if the compiler engineers implement an intrinsic, or auto-vectorisation, which will try, and sometimes succeed admirably, to translate your code to a good vector implementation. There is currently a prototype project for explicit vectorisation in Project Panama. There are a few ways to skin this cat, and it’s worth looking at the code they generate and the throughput available from each approach.

There is a choice between copying the array and zeroing out the negatives, and allocating a new array and only writing the non-negative values. There is another choice between an if statement and branchless code using Math.max. This results in the following four implementations which I measure on comparable data to the C++ benchmark (10 million doubles, normally distributed with mean zero). To be fair to the Java code, as in the C++ benchmarks, the cost of allocation is isolated by writing into an array pre-allocated once per benchmark. This penalises the approaches where the array is copied first and then zeroed wherever the value is negative. The code is online at github.

    @Benchmark
    @CompilerControl(CompilerControl.Mode.DONT_INLINE)
    public double[] BranchyCopyAndMask(ArrayWithNegatives state) {
        double[] data = state.data;
        double[] result = state.target;
        System.arraycopy(data, 0, result, 0, data.length);
        for (int i = 0; i < result.length; ++i) {
            if (result[i] < 0D) {
                result[i] = 0D;
            }
        }
        return result;
    }

    @Benchmark
    @CompilerControl(CompilerControl.Mode.DONT_INLINE)
    public double[] BranchyNewArray(ArrayWithNegatives state) {
        double[] data = state.data;
        double[] result = state.target;
        for (int i = 0; i < result.length; ++i) {
            result[i] = data[i] < 0D ? 0D : data[i];
        }
        return result;
    }

    @Benchmark
    @CompilerControl(CompilerControl.Mode.DONT_INLINE)
    public double[] NewArray(ArrayWithNegatives state) {
        double[] data = state.data;
        double[] result = state.target;
        for (int i = 0; i < result.length; ++i) {
            result[i] = Math.max(data[i], 0D);
        }
        return result;
    }

    @Benchmark
    @CompilerControl(CompilerControl.Mode.DONT_INLINE)
    public double[] CopyAndMask(ArrayWithNegatives state) {
        double[] data = state.data;
        double[] result = state.target;
        System.arraycopy(data, 0, result, 0, data.length);
        for (int i = 0; i < result.length; ++i) {
            result[i] = Math.max(result[i], 0D);
        }
        return result;
    }

None of these implementations comes close to the native code above. The best implementation performs 1.8 iterations per second which equates to processing approximately 1.4GB/s, vastly inferior to the 4GB/s achieved with Intel intrinsics. The results are below:

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit
BranchyCopyAndMask thrpt 1 10 1.314845 0.061662 ops/s
BranchyNewArray thrpt 1 10 1.802673 0.061835 ops/s
CopyAndMask thrpt 1 10 1.146630 0.018903 ops/s
NewArray thrpt 1 10 1.357020 0.116481 ops/s

As an aside, there is a very interesting observation to make, worthy of its own post: if the array consists only of positive values, the “branchy” implementations run very well, at speeds comparable to the zero_negatives (when it ran with 50% negatives). The ratio of branch hits to misses is an orthogonal explanatory variable, and the input data, while I often don’t think about it enough, is very important.

I only looked at the assembly emitted for the fastest version (BranchyNewArray) and it doesn’t look anything like zero_negatives, though it does use some vectorisation – as pointed out by Daniel Lemire in the comments, this code has probably not been vectorised and is probably using SSE2 (indeed only quad words are loaded into 128 bit registers):

com/openkappa/simd/positive/PositiveValues.BranchyNewArray(Lcom/openkappa/simd/positive/ArrayWithNegatives;)[D  [0x000002ae309c3ce0, 0x000002ae309c3ff8]  792 bytes
Argument 0 is unknown.RIP: 0x2ae309c3ce0 Code size: 0x00000318
[Entry Point]
[Constants]
  # {method} {0x000002ae4d13ec70} 'BranchyNewArray' '(Lcom/openkappa/simd/positive/ArrayWithNegatives;)[D' in 'com/openkappa/simd/positive/PositiveValues'
  0x000002ae309c3ce0: mov     r10d,dword ptr [rdx+8h]  ;...44
                                                ;...8b
                                                ;...52
                                                ;...08

  0x000002ae309c3ce4: shl     r10,3h            ;...49
                                                ;...c1
                                                ;...e2
                                                ;...03

  0x000002ae309c3ce8: cmp     r10,rax           ;...4c
                                                ;...3b
                                                ;...d0

  0x000002ae309c3ceb: jne     2ae3042c200h      ;...0f
                                                ;...85
                                                ;...0f
                                                ;...85
                                                ;...a6
                                                ;...ff
                                                ;   {runtime_call ic_miss_stub}
  0x000002ae309c3cf1: nop     word ptr [rax+rax+0h]  ;...66
                                                ;...66
                                                ;...66
                                                ;...0f
                                                ;...1f
                                                ;...84
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002ae309c3cfc: nop                       ;...66
                                                ;...66
                                                ;...66
                                                ;...90

[Verified Entry Point]
  0x000002ae309c3d00: mov     dword ptr [rsp+0ffffffffffff9000h],eax
                                                ;...89
                                                ;...84
                                                ;...24
                                                ;...00
                                                ;...90
                                                ;...ff
                                                ;...ff

  0x000002ae309c3d07: push    rbp               ;...55

  0x000002ae309c3d08: sub     rsp,60h           ;...48
                                                ;...83
                                                ;...ec
                                                ;...60

  0x000002ae309c3d0c: mov     rcx,2ae4d163880h  ;...48
                                                ;...b9
                                                ;...80
                                                ;...38
                                                ;...16
                                                ;...4d
                                                ;...ae
                                                ;...02
                                                ;...00
                                                ;...00
                                                ;   {metadata(method data for {method} {0x000002ae4d13ec70} 'BranchyNewArray' '(Lcom/openkappa/simd/positive/ArrayWithNegatives;)[D' in 'com/openkappa/simd/positive/PositiveValues')}
  0x000002ae309c3d16: mov     esi,dword ptr [rcx+0fch]
                                                ;...8b
                                                ;...b1
                                                ;...fc
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002ae309c3d1c: add     esi,8h            ;...83
                                                ;...c6
                                                ;...08

  0x000002ae309c3d1f: mov     dword ptr [rcx+0fch],esi
                                                ;...89
                                                ;...b1
                                                ;...fc
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002ae309c3d25: and     esi,1ff8h         ;...81
                                                ;...e6
                                                ;...f8
                                                ;...1f
                                                ;...00
                                                ;...00

  0x000002ae309c3d2b: cmp     esi,0h            ;...83
                                                ;...fe
                                                ;...00

  0x000002ae309c3d2e: je      2ae309c3ec1h      ;...0f
                                                ;...84
                                                ;...8d
                                                ;...01
                                                ;...00
                                                ;...00
                                                ;*aload_1 {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.positive.PositiveValues::BranchyNewArray@0 (line 29)

  0x000002ae309c3d34: mov     edx,dword ptr [r8+0ch]  ;...41
                                                ;...8b
                                                ;...50
                                                ;...0c
                                                ; implicit exception: dispatches to 0x000002ae309c3ee2
  0x000002ae309c3d38: shl     rdx,3h            ;...48
                                                ;...c1
                                                ;...e2
                                                ;...03
                                                ;*getfield data {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.positive.PositiveValues::BranchyNewArray@1 (line 29)

  0x000002ae309c3d3c: mov     ecx,dword ptr [r8+10h]  ;...41
                                                ;...8b
                                                ;...48
                                                ;...10

  0x000002ae309c3d40: shl     rcx,3h            ;...48
                                                ;...c1
                                                ;...e1
                                                ;...03
                                                ;*getfield target {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.positive.PositiveValues::BranchyNewArray@6 (line 30)

  0x000002ae309c3d44: mov     esi,0h            ;...be
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002ae309c3d49: jmp     2ae309c3e27h      ;...e9
                                                ;...d9
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;*iload {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.positive.PositiveValues::BranchyNewArray@13 (line 31)

  0x000002ae309c3d4e: nop                       ;...66
                                                ;...90

  0x000002ae309c3d50: movsxd  rax,esi           ;...48
                                                ;...63
                                                ;...c6

  0x000002ae309c3d53: cmp     esi,dword ptr [rdx+0ch]  ;...3b
                                                ;...72
                                                ;...0c
                                                ; implicit exception: dispatches to 0x000002ae309c3ee7
  0x000002ae309c3d56: jnb     2ae309c3ef1h      ;...0f
                                                ;...83
                                                ;...95
                                                ;...01
                                                ;...00
                                                ;...00

  0x000002ae309c3d5c: vmovsd  xmm0,qword ptr [rdx+rax*8+10h]
                                                ;...c5
                                                ;...fb
                                                ;...10
                                                ;...44
                                                ;...c2
                                                ;...10
                                                ;*daload {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.positive.PositiveValues::BranchyNewArray@26 (line 32)

  0x000002ae309c3d62: vxorpd  xmm1,xmm1,xmm1    ;...c5
                                                ;...f1
                                                ;...57
                                                ;...c9

  0x000002ae309c3d66: vucomisd xmm0,xmm1        ;...c5
                                                ;...f9
                                                ;...2e
                                                ;...c1

  0x000002ae309c3d6a: mov     eax,1h            ;...b8
                                                ;...01
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002ae309c3d6f: jp      2ae309c3d88h      ;...0f
                                                ;...8a
                                                ;...13
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002ae309c3d75: jnbe    2ae309c3d88h      ;...0f
                                                ;...87
                                                ;...0d
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002ae309c3d7b: mov     eax,0h            ;...b8
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002ae309c3d80: je      2ae309c3d88h      ;...0f
                                                ;...84
                                                ;...02
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002ae309c3d86: dec     eax               ;...ff
                                                ;...c8
                                                ;*dcmpg {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.positive.PositiveValues::BranchyNewArray@28 (line 32)

  0x000002ae309c3d88: cmp     eax,0h            ;...83
                                                ;...f8
                                                ;...00

  0x000002ae309c3d8b: mov     rax,2ae4d163880h  ;...48
                                                ;...b8
                                                ;...80
                                                ;...38
                                                ;...16
                                                ;...4d
                                                ;...ae
                                                ;...02
                                                ;...00
                                                ;...00
                                                ;   {metadata(method data for {method} {0x000002ae4d13ec70} 'BranchyNewArray' '(Lcom/openkappa/simd/positive/ArrayWithNegatives;)[D' in 'com/openkappa/simd/positive/PositiveValues')}
  0x000002ae309c3d95: mov     rdi,158h          ;...48
                                                ;...bf
                                                ;...58
                                                ;...01
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002ae309c3d9f: jnl     2ae309c3dafh      ;...0f
                                                ;...8d
                                                ;...0a
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002ae309c3da5: mov     rdi,168h          ;...48
                                                ;...bf
                                                ;...68
                                                ;...01
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002ae309c3daf: mov     rbx,qword ptr [rax+rdi]  ;...48
                                                ;...8b
                                                ;...1c
                                                ;...38

  0x000002ae309c3db3: lea     rbx,[rbx+1h]      ;...48
                                                ;...8d
                                                ;...5b
                                                ;...01

  0x000002ae309c3db7: mov     qword ptr [rax+rdi],rbx  ;...48
                                                ;...89
                                                ;...1c
                                                ;...38

  0x000002ae309c3dbb: jnl     2ae309c3dd5h      ;...0f
                                                ;...8d
                                                ;...14
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;*ifge {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.positive.PositiveValues::BranchyNewArray@29 (line 32)

  0x000002ae309c3dc1: mov     rax,2ae4d163880h  ;...48
                                                ;...b8
                                                ;...80
                                                ;...38
                                                ;...16
                                                ;...4d
                                                ;...ae
                                                ;...02
                                                ;...00
                                                ;...00
                                                ;   {metadata(method data for {method} {0x000002ae4d13ec70} 'BranchyNewArray' '(Lcom/openkappa/simd/positive/ArrayWithNegatives;)[D' in 'com/openkappa/simd/positive/PositiveValues')}
  0x000002ae309c3dcb: inc     dword ptr [rax+178h]  ;...ff
                                                ;...80
                                                ;...78
                                                ;...01
                                                ;...00
                                                ;...00

  0x000002ae309c3dd1: vxorpd  xmm0,xmm0,xmm0    ;...c5
                                                ;...f9
                                                ;...57
                                                ;...c0
                                                ;*goto {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.positive.PositiveValues::BranchyNewArray@33 (line 32)

  0x000002ae309c3dd5: movsxd  rax,esi           ;...48
                                                ;...63
                                                ;...c6

  0x000002ae309c3dd8: cmp     esi,dword ptr [rcx+0ch]  ;...3b
                                                ;...71
                                                ;...0c

  0x000002ae309c3ddb: jnb     2ae309c3efah      ;...0f
                                                ;...83
                                                ;...19
                                                ;...01
                                                ;...00
                                                ;...00

  0x000002ae309c3de1: vmovsd  qword ptr [rcx+rax*8+10h],xmm0
                                                ;...c5
                                                ;...fb
                                                ;...11
                                                ;...44
                                                ;...c1
                                                ;...10
                                                ;*dastore {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.positive.PositiveValues::BranchyNewArray@40 (line 32)

  0x000002ae309c3de7: inc     esi               ;...ff
                                                ;...c6

  0x000002ae309c3de9: mov     rax,2ae4d163880h  ;...48
                                                ;...b8
                                                ;...80
                                                ;...38
                                                ;...16
                                                ;...4d
                                                ;...ae
                                                ;...02
                                                ;...00
                                                ;...00
                                                ;   {metadata(method data for {method} {0x000002ae4d13ec70} 'BranchyNewArray' '(Lcom/openkappa/simd/positive/ArrayWithNegatives;)[D' in 'com/openkappa/simd/positive/PositiveValues')}
  0x000002ae309c3df3: mov     edi,dword ptr [rax+100h]
                                                ;...8b
                                                ;...b8
                                                ;...00
                                                ;...01
                                                ;...00
                                                ;...00

  0x000002ae309c3df9: add     edi,8h            ;...83
                                                ;...c7
                                                ;...08

  0x000002ae309c3dfc: mov     dword ptr [rax+100h],edi
                                                ;...89
                                                ;...b8
                                                ;...00
                                                ;...01
                                                ;...00
                                                ;...00

  0x000002ae309c3e02: and     edi,0fff8h        ;...81
                                                ;...e7
                                                ;...f8
                                                ;...ff
                                                ;...00
                                                ;...00

  0x000002ae309c3e08: cmp     edi,0h            ;...83
                                                ;...ff
                                                ;...00

  0x000002ae309c3e0b: je      2ae309c3f03h      ;...0f
                                                ;...84
                                                ;...f2
                                                ;...00
                                                ;...00
                                                ;...00
                                                ; ImmutableOopMap{rcx=Oop rdx=Oop }
                                                ;*goto {reexecute=1 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.positive.PositiveValues::BranchyNewArray@44 (line 31)

  0x000002ae309c3e11: test    dword ptr [2ae24520000h],eax
                                                ;...85
                                                ;...05
                                                ;...e9
                                                ;...c1
                                                ;...b5
                                                ;...f3
                                                ;   {poll}
  0x000002ae309c3e17: mov     rax,2ae4d163880h  ;...48
                                                ;...b8
                                                ;...80
                                                ;...38
                                                ;...16
                                                ;...4d
                                                ;...ae
                                                ;...02
                                                ;...00
                                                ;...00
                                                ;   {metadata(method data for {method} {0x000002ae4d13ec70} 'BranchyNewArray' '(Lcom/openkappa/simd/positive/ArrayWithNegatives;)[D' in 'com/openkappa/simd/positive/PositiveValues')}
  0x000002ae309c3e21: inc     dword ptr [rax+190h]  ;...ff
                                                ;...80
                                                ;...90
                                                ;...01
                                                ;...00
                                                ;...00
                                                ;*goto {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.positive.PositiveValues::BranchyNewArray@44 (line 31)

  0x000002ae309c3e27: mov     eax,dword ptr [rcx+0ch]  ;...8b
                                                ;...41
                                                ;...0c
                                                ;*arraylength {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.positive.PositiveValues::BranchyNewArray@16 (line 31)
                                                ; implicit exception: dispatches to 0x000002ae309c3f24
  0x000002ae309c3e2a: cmp     esi,eax           ;...3b
                                                ;...f0

  0x000002ae309c3e2c: mov     rax,2ae4d163880h  ;...48
                                                ;...b8
                                                ;...80
                                                ;...38
                                                ;...16
                                                ;...4d
                                                ;...ae
                                                ;...02
                                                ;...00
                                                ;...00
                                                ;   {metadata(method data for {method} {0x000002ae4d13ec70} 'BranchyNewArray' '(Lcom/openkappa/simd/positive/ArrayWithNegatives;)[D' in 'com/openkappa/simd/positive/PositiveValues')}
  0x000002ae309c3e36: mov     rdi,148h          ;...48
                                                ;...bf
                                                ;...48
                                                ;...01
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002ae309c3e40: jl      2ae309c3e50h      ;...0f
                                                ;...8c
                                                ;...0a
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002ae309c3e46: mov     rdi,138h          ;...48
                                                ;...bf
                                                ;...38
                                                ;...01
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002ae309c3e50: mov     rbx,qword ptr [rax+rdi]  ;...48
                                                ;...8b
                                                ;...1c
                                                ;...38

  0x000002ae309c3e54: lea     rbx,[rbx+1h]      ;...48
                                                ;...8d
                                                ;...5b
                                                ;...01

  0x000002ae309c3e58: mov     qword ptr [rax+rdi],rbx  ;...48
                                                ;...89
                                                ;...1c
                                                ;...38

  0x000002ae309c3e5c: jl      2ae309c3d50h      ;...0f
                                                ;...8c
                                                ;...ee
                                                ;...fe
                                                ;...ff
                                                ;...ff
                                                ;*if_icmpge {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.positive.PositiveValues::BranchyNewArray@17 (line 31)

  0x000002ae309c3e62: mov     rax,rcx           ;...48
                                                ;...8b
                                                ;...c1

  0x000002ae309c3e65: add     rsp,60h           ;...48
                                                ;...83
                                                ;...c4
                                                ;...60

  0x000002ae309c3e69: pop     rbp               ;...5d

  0x000002ae309c3e6a: test    dword ptr [2ae24520000h],eax
                                                ;...85
                                                ;...05
                                                ;...90
                                                ;...c1
                                                ;...b5
                                                ;...f3
                                                ;   {poll_return}
  0x000002ae309c3e70: ret                       ;...c3
                                                ;*areturn {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.simd.positive.PositiveValues::BranchyNewArray@48 (line 34)

I don’t really understand, and haven’t thought about, the intent of the emitted code, but it makes extensive use of the instruction VUCOMISD for comparisons with zero, which has a lower latency but lower throughput than VMAXPD. It would certainly be interesting to see how Project Panama does this. Perhaps this should just be made available as a fail-safe intrinsic like Arrays.mismatch?

Still True in Java 9: Handwritten Hash Codes are Faster

One of the things to keep in mind with Java is that the best performance advice for one version may not apply to the next. The JVM has an optimiser which works by detecting intent in byte-code; it does a better job when the programmer is skilled enough to make that intent clear. There have been times when the optimiser has done a bad job, either through bugs or feature gaps, and compensatory idioms emerge. The value of these idioms degrades over time as the optimiser improves; a stroke of ingenuity in one version can become ritualistic nonsense in the next. It’s important not to be that guy who fears adding strings together because it was costly a decade ago, but does it always get better, even with low hanging fruit?

I reproduced the results of an extremely astute optimisation presented in 2015 by Daniel Lemire. I was hoping to see that an improved optimiser in Java 9, having observed several cases of automatic vectorisation, would render this optimisation null.

Hash Codes for Arrays

I encourage you to read the original blog post because it is informative, but for the sake of coherence I will summarise the key point here. Arrays.hashCode implements the following hash code:

    public static int hashCode(int[] a) {
        if (a == null)
            return 0;

        int result = 1;
        for (int element : a)
            result = 31 * result + element;

        return result;
    }

This results in a good hash code, but a scalar internal representation of this code has a problem: a data dependency on the multiplication, which is slower than moving data into registers. A significant and reproducible speed up can be observed when unrolling the loop, which enforces a prefetch of 128 bits of the array into registers before doing the slower multiplications:

    public static int hashCode(int[] a) {
        if (a == null)
            return 0;

        int result = 1;
        int i = 0;
        for (; i + 3 < a.length; i += 4) {
            result = 31 * 31 * 31 * 31 * result
                   + 31 * 31 * 31 * a[i]
                   + 31 * 31 * a[i + 1]
                   + 31 * a[i + 2]
                   + a[i + 3];
        }
        for (; i < a.length; i++) {
            result = 31 * result + a[i];
        }
        return result;
    }

The improvement in performance from this simple change can be confirmed with Java 8 by running a simple benchmark.

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size
BuiltIn thrpt 1 10 9.537736 0.382617 ops/us 100
BuiltIn thrpt 1 10 0.804620 0.103037 ops/us 1000
BuiltIn thrpt 1 10 0.092297 0.003947 ops/us 10000
Unrolled thrpt 1 10 14.974398 0.628522 ops/us 100
Unrolled thrpt 1 10 1.514986 0.035759 ops/us 1000
Unrolled thrpt 1 10 0.137408 0.010200 ops/us 10000

The performance improvement is so obvious, and the change so easy to make, that one wonders why JVM vendors didn’t make the change themselves.

Java 9: Universally Better Automatic Optimisations?

As the comments on the blog post suggest, this is a prime candidate for vectorisation. Auto-vectorisation is a thing in Java 9. Using intrinsics or code clean enough to express intent clearly, you can really expect to see good usage of SIMD in Java 9. I was really expecting the situation to reverse in Java 9; but it doesn’t.

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size
BuiltIn thrpt 1 10 9.822568 0.381087 ops/us 100
BuiltIn thrpt 1 10 0.949273 0.024021 ops/us 1000
BuiltIn thrpt 1 10 0.093171 0.004502 ops/us 10000
Unrolled thrpt 1 10 13.762617 0.440135 ops/us 100
Unrolled thrpt 1 10 1.501106 0.094897 ops/us 1000
Unrolled thrpt 1 10 0.139963 0.011487 ops/us 10000

This is a still a smart optimisation two years later, but it offends my sensibilities in the same way hints do in SQL – a layer of abstraction has been punctured. I will have to try again in Java 10.

New Methods in Java 9: Math.fma and Arrays.mismatch

There are two noteworthy new methods in Java 9: Arrays.mismatch and Math.fma.

Arrays.mismatch

This method takes two primitive arrays, and returns the index of the first differing values. This effectively computes the longest common prefix of the two arrays. This is really quite useful, mostly for text processing but also for Bioinformatics (protein sequencing and so on, much more interesting than the sort of thing I work on). Having worked extensively with Apache HBase (where a vast majority of the API involves manipulating byte arrays) I can think of lots of less interesting use cases for this method.

Looking carefully, you can see that the method calls into the internal ArraysSupport utility class, which will try to perform a vectorised mismatch (an intrinsic candidate). Since this will use AVX instructions, this is very fast; much faster than a handwritten loop.

Let’s measure the boost versus a handwritten loop, testing across a range of common prefices and array lengths of byte[].

    @Benchmark
    @CompilerControl(CompilerControl.Mode.DONT_INLINE)
    public void Mismatch_Intrinsic(BytePrefixData data, Blackhole bh) {
        bh.consume(Arrays.mismatch(data.data1, data.data2));
    }


    @Benchmark
    @CompilerControl(CompilerControl.Mode.DONT_INLINE)
    public void Mismatch_Handwritten(BytePrefixData data, Blackhole bh) {
        byte[] data1 = data.data1;
        byte[] data2 = data.data2;
        int length = Math.min(data1.length, data2.length);
        int mismatch = -1;
        for (int i = 0; i < length; ++i) {
            if (data1[i] != data2[i]) {
                mismatch = i;
                break;
            }
        }
        bh.consume(mismatch);
    }

The results speak for themselves. Irritatingly, there is some duplication in output because I haven’t figured out how to make JMH use a subset of the Cartesian product of its parameters.

Benchmark (prefix) (size) Mode Cnt Score Error Units
Mismatch_Handwritten 10 100 thrpt 10 22.360 ± 0.938 ops/us
Mismatch_Handwritten 10 1000 thrpt 10 2.459 ± 0.256 ops/us
Mismatch_Handwritten 10 10000 thrpt 10 0.255 ± 0.009 ops/us
Mismatch_Handwritten 100 100 thrpt 10 22.763 ± 0.869 ops/us
Mismatch_Handwritten 100 1000 thrpt 10 2.690 ± 0.044 ops/us
Mismatch_Handwritten 100 10000 thrpt 10 0.273 ± 0.008 ops/us
Mismatch_Handwritten 1000 100 thrpt 10 24.970 ± 0.713 ops/us
Mismatch_Handwritten 1000 1000 thrpt 10 2.791 ± 0.066 ops/us
Mismatch_Handwritten 1000 10000 thrpt 10 0.281 ± 0.007 ops/us
Mismatch_Intrinsic 10 100 thrpt 10 89.169 ± 2.759 ops/us
Mismatch_Intrinsic 10 1000 thrpt 10 26.995 ± 0.501 ops/us
Mismatch_Intrinsic 10 10000 thrpt 10 3.553 ± 0.065 ops/us
Mismatch_Intrinsic 100 100 thrpt 10 83.037 ± 5.590 ops/us
Mismatch_Intrinsic 100 1000 thrpt 10 26.249 ± 0.714 ops/us
Mismatch_Intrinsic 100 10000 thrpt 10 3.523 ± 0.122 ops/us
Mismatch_Intrinsic 1000 100 thrpt 10 87.921 ± 6.566 ops/us
Mismatch_Intrinsic 1000 1000 thrpt 10 25.812 ± 0.442 ops/us
Mismatch_Intrinsic 1000 10000 thrpt 10 4.177 ± 0.059 ops/us

Why is there such a big difference? Look at how the score decreases as a function of array length, even when the common prefix, and therefore the number of comparisons required, is small: clearly the performance of this algorithm depends on the efficiency of memory access. Arrays.mismatch optimises this, reading qwords of the array into SIMD registers. Working one long at a time, it is possible to compute the XOR in a single instruction to determine if it’s even necessary to look at each byte.

java/util/ArraysSupport.vectorizedMismatch(Ljava/lang/Object;JLjava/lang/Object;JII)I  [0x000002bd9215a820, 0x000002bd9215aa78]  600 bytes
Argument 0 is unknown.RIP: 0x2bd9215a820 Code size: 0x00000258
[Entry Point]
[Verified Entry Point]
[Constants]
  # {method} {0x000002bda79cbf68} 'vectorizedMismatch' '(Ljava/lang/Object;JLjava/lang/Object;JII)I' in 'java/util/ArraysSupport'
  # parm0:    rdx:rdx   = 'java/lang/Object'
  # parm1:    r8:r8     = long
  # parm2:    r9:r9     = 'java/lang/Object'
  # parm3:    rdi:rdi   = long
  # parm4:    rsi       = int
  # parm5:    rcx       = int
  #           [sp+0x60]  (sp of caller)
  0x000002bd9215a820: mov     dword ptr [rsp+0ffffffffffff9000h],eax
                                                ;...89
                                                ;...84
                                                ;...24
                                                ;...00
                                                ;...90
                                                ;...ff
                                                ;...ff

  0x000002bd9215a827: push    rbp               ;...55

  0x000002bd9215a828: sub     rsp,50h           ;...48
                                                ;...83
                                                ;...ec
                                                ;...50
                                                ;*synchronization entry
                                                ; - java.util.ArraysSupport::vectorizedMismatch@-1 (line 120)

  0x000002bd9215a82c: mov     r10,rdi           ;...4c
                                                ;...8b
                                                ;...d7

  0x000002bd9215a82f: vmovq   xmm2,r9           ;...c4
                                                ;...c1
                                                ;...f9
                                                ;...6e
                                                ;...d1

  0x000002bd9215a834: vmovq   xmm1,rdx          ;...c4
                                                ;...e1
                                                ;...f9
                                                ;...6e
                                                ;...ca

  0x000002bd9215a839: mov     r14d,ecx          ;...44
                                                ;...8b
                                                ;...f1

  0x000002bd9215a83c: vmovd   xmm0,esi          ;...c5
                                                ;...f9
                                                ;...6e
                                                ;...c6

  0x000002bd9215a840: mov     r9d,3h            ;...41
                                                ;...b9
                                                ;...03
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002bd9215a846: sub     r9d,ecx           ;...44
                                                ;...2b
                                                ;...c9
                                                ;*isub {reexecute=0 rethrow=0 return_oop=0}
                                                ; - java.util.ArraysSupport::vectorizedMismatch@5 (line 120)

  0x000002bd9215a849: mov     edx,esi           ;...8b
                                                ;...d6

  0x000002bd9215a84b: mov     ecx,r9d           ;...41
                                                ;...8b
                                                ;...c9

  0x000002bd9215a84e: sar     edx,cl            ;...d3
                                                ;...fa
                                                ;*ishr {reexecute=0 rethrow=0 return_oop=0}
                                                ; - java.util.ArraysSupport::vectorizedMismatch@17 (line 122)

  0x000002bd9215a850: mov     eax,1h            ;...b8
                                                ;...01
                                                ;...00
                                                ;...00
                                                ;...00

  0x000002bd9215a855: xor     edi,edi           ;...33
                                                ;...ff

  0x000002bd9215a857: test    edx,edx           ;...85
                                                ;...d2

  0x000002bd9215a859: jle     2bd9215a97ah      ;...0f
                                                ;...8e
                                                ;...1b
                                                ;...01
                                                ;...00
                                                ;...00

The code for this benchmark is at github.

Math.fma

In comparison to users of some languages, Java programmers are lackadaisical about floating point errors. It’s a good job that historically Java hasn’t been considered suitable for the implementation of numerical algorithms. But all of a sudden there is a revolution of data science on the JVM, albeit mostly driven by the Scala community, with JVM implementations of structures like recurrent neural networks abounding. It matters less for machine learning than root finding, but how accurate can these implementations be without JVM level support for minimising the propagation floating point errors? With Math.fma this is improving, by allowing two common operations to be performed before rounding.

Math.fma fuses a multiplication and an addition into a single floating point operation to compute expressions like ab + c. This has two key benefits:

  1. There’s only one operation, and only one rounding error
  2. This is explicitly supported in AVX2 by the VFMADD* instructions

Newton’s Method

To investigate any superior suppression of floating point errors, I use a toy implementation of Newton’s method to compute the root of a quadratic equation, which any teenager could calculate analytically (the error is easy to quantify).

I compare these two implementations for 4x^2 - 12x + 9 (there is a repeated root at 1.5) to get an idea for the error (defined by |1.5 - x_n|) after a large number of iterations.

I implemented this using FMA:

public class NewtonsMethodFMA {

    private final double[] coefficients;

    public NewtonsMethodFMA(double[] coefficients) {
        this.coefficients = coefficients;
    }


    public double evaluateF(double x) {
        double f = 0D;
        int power = coefficients.length - 1;
        for (int i = 0; i < coefficients.length; ++i) {
            f = Math.fma(coefficients[i], Math.pow(x, power--), f);
        }
        return f;
    }

    public double evaluateDF(double x) {
        double df = 0D;
        int power = coefficients.length - 2;
        for (int i = 0; i < coefficients.length - 1; ++i) {
            df = Math.fma((power + 1) * coefficients[i],  Math.pow(x, power--), df);
        }
        return df;
    }

    public double solve(double initialEstimate, int maxIterations) {
        double result = initialEstimate;
        for (int i = 0; i < maxIterations; ++i) {
            result -= evaluateF(result)/evaluateDF(result);
        }
        return result;
    }
}

And an implementation with normal operations:


public class NewtonsMethod {

    private final double[] coefficients;

    public NewtonsMethod(double[] coefficients) {
        this.coefficients = coefficients;
    }


    public double evaluateF(double x) {
        double f = 0D;
        int power = coefficients.length - 1;
        for (int i = 0; i < coefficients.length; ++i) {
            f += coefficients[i] * Math.pow(x, power--);
        }
        return f;
    }

    public double evaluateDF(double x) {
        double df = 0D;
        int power = coefficients.length - 2;
        for (int i = 0; i < coefficients.length - 1; ++i) {
            df += (power + 1) * coefficients[i] * Math.pow(x, power--);
        }
        return df;
    }

    public double solve(double initialEstimate, int maxIterations) {
        double result = initialEstimate;
        for (int i = 0; i < maxIterations; ++i) {
            result -= evaluateF(result)/evaluateDF(result);
        }
        return result;
    }
}

When I run this code for 1000 iterations, the FMA version results in 1.5000000083575202, whereas the vanilla version results in 1.500000017233207. It’s completely unscientific, but seems plausible and confirms my prejudice so… In fact, it’s not that simple, and over a range of initial values, there is only a very small difference in FMA’s favour. There’s not even a performance improvement – clearly this method wasn’t added so you can start implementing numerical root finding algorithms – the key takeaway is that the results are slightly different because a different rounding strategy has been used.

Benchmark (maxIterations) Mode Cnt Score Error Units
NM_FMA 100 thrpt 10 93.805 ± 5.174 ops/ms
NM_FMA 1000 thrpt 10 9.420 ± 1.169 ops/ms
NM_FMA 10000 thrpt 10 0.962 ± 0.044 ops/ms
NM_HandWritten 100 thrpt 10 93.457 ± 5.048 ops/ms
NM_HandWritten 1000 thrpt 10 9.274 ± 0.483 ops/ms
NM_HandWritten 10000 thrpt 10 0.928 ± 0.041 ops/ms

Choosing the Right Radix: Measurement or Mathematics?

I recently wrote a post about radix sorting, and found that for large arrays of unsigned integers a handwritten implementation beats Arrays.sort. However, I paid no attention to the choice of radix and used a default of eight bits. It turns out this was a lucky choice: modifying my benchmark to parametrise the radix, I observed a maximum at one byte, regardless of the size of array.

Is this an algorithmic or technical phenomenon? Is this something that could have been predicted on the back of an envelope without running a benchmark?

Extended Benchmark Results

Size Radix Score Score Error (99.9%) Unit
100 2 135.559923 7.72397 ops/ms
100 4 262.854579 37.372678 ops/ms
100 8 345.038234 30.954927 ops/ms
100 16 7.717496 1.144967 ops/ms
1000 2 13.892142 1.522749 ops/ms
1000 4 27.712719 4.057162 ops/ms
1000 8 52.253497 4.761172 ops/ms
1000 16 7.656033 0.499627 ops/ms
10000 2 1.627354 0.186948 ops/ms
10000 4 3.620869 0.029128 ops/ms
10000 8 6.435789 0.610848 ops/ms
10000 16 3.703248 0.45177 ops/ms
100000 2 0.144575 0.014348 ops/ms
100000 4 0.281837 0.025707 ops/ms
100000 8 0.543274 0.031553 ops/ms
100000 16 0.533998 0.126949 ops/ms
1000000 2 0.011293 0.001429 ops/ms
1000000 4 0.021128 0.003137 ops/ms
1000000 8 0.037376 0.005783 ops/ms
1000000 16 0.031053 0.007987 ops/ms

Modeling

To model the execution time of the algorithm, we can write t = f(r, n), where n \in \mathbb{N} is the length of the input array, and r \in [1, 32) is the size in bits of the radix. We can inspect if the model predicts non-monotonic execution time with a minimum (corresponding to maximal throughput), or if t increases indefinitely as a function of r. If we find a plausible model predicting a minimum, temporarily treating r as continuous, we can solve \frac{\partial f}{\partial r}|_{n=N, r \in [1,32)} = 0 to find the theoretically optimal radix. This pre-supposes we derive a non-monotonic model.

Constructing a Model

We need to write down an equation before we can do any calculus, which requires two dangerous assumptions.

  1. Each operation has the same cost, making the execution time proportional to the number of operation.
  2. The costs of operations do not vary as a function of n or r.

This means all we need to do is find a formula for the number of operations, and then vary n and r. The usual pitfall in this approach relates to the first assumption, in that memory accesses are modelled as uniform cost; memory access can vary widely in cost ranging from registers to RAM on another socket. We are about to fall foul of both assumptions constructing an intuitive model of the algorithm’s loop.

        while (shift < Integer.SIZE) {
            Arrays.fill(histogram, 0);
            for (int i = 0; i < data.length; ++i) {
                ++histogram[((data[i] & mask) >> shift) + 1];
            }
            for (int i = 0; i < 1 << radix; ++i) {
                histogram[i + 1] += histogram[i];
            }
            for (int i = 0; i < data.length; ++i) {
                copy[histogram[(data[i] & mask) >> shift]++] = data[i];
            }
            for (int i = 0; i < data.length; ++i) {
                data[i] = copy[i];
            }
            shift += radix;
            mask <<= radix;
        }

The outer loop depends on the choice of radix while the inner loops depend on the size of the array and the choice of radix. There are five obvious aspects to capture:

  • The first inner loop takes time proportional to n
  • The third and fourth inner loops take time proportional to n
  • We can factor the per-element costs of loops 1, 3 and 4 into a constant a
  • The second inner loop takes time proportional to 2^r, modeled with by the term b2^r
  • The body of the loop executes 32/r times

This can be summarised as the formula:

f(r, n) = 32\frac{(3an + b2^r)}{r}

It was claimed the algorithm had linear complexity in n and it only has a linear term in n. Good. However, the exponential r term in the numerator dominates the linear term in the denominator, making the function monotonic in r. The model fails to predict the observed throughput maximising radix. There are clearly much more complicated mechanisms at play than can be captured counting operations.

Sorting Unsigned Integers Faster in Java

I discovered a curious resource for audio-visualising sort algorithms, which is exciting for two reasons. The first is that I finally feel like I understand Alexander Scriabin: he was not a composer. He discovered Tim Sort 80 years before Tim Peters and called it Black Mass. (If you aren’t familiar with the piece, fast-forward to 1:40 to hear the congruence.)

The second reason was that I noticed Radix Sort (LSD). While it was an affront to my senses, it used a mere 800 array accesses and no comparisons! I was unaware of this algorithm so delved deeper and implemented it for integers, and benchmarked my code against Arrays.sort.

Radix Sort

It is taken as given by many (myself included, or am I just projecting my thoughts on to others?) that O(n \log n) is the best you can do in a sort algorithm. But this is actually only true for sort algorithms which depend on comparison. If you can afford to restrict the data types your sort algorithm supports to types with a positional interpretation (java.util can’t because it needs to be ubiquitous and maintainable), you can get away with a linear time algorithm.

Radix sort, along with the closely related counting sort, does not use comparisons. Instead, the data is interpreted as a fixed length string of symbols. For each position, the cumulative histogram of symbols is computed to calculate sort indices. While the data needs to be scanned several times, the algorithm scales linearly and the overhead of the multiple scans is amortised for large arrays.

As you can see on Wikipedia, there are two kinds of radix sort: Least Significant Digit and Most Significant Digit. This dichotomy relates to the order the (representational) string of symbols is traversed in. I implemented and benchmarked the LSD version for integers.

Implementation

The implementation interprets an integer as the concatenation of n bit string symbols of fixed size size 32/n. It performs n passes over the array, starting with the least significant bits, which it modifies in place. For each pass the data is scanned three times, in order to:

  1. Compute the cumulative histogram over the symbols in their natural sort order
  2. Copy the value with symbol k to the mth position in a buffer, where m is defined by the cumulative density of k.
  3. Copy the buffer back into the original array

The implementation, which won’t work unless the chunks are proper divisors of 32, is below. The bonus (or caveat) is that it automatically supports unsigned integers. The code could be modified slightly to work with signed integers at a performance cost.

import java.util.Arrays;

public class RadixSort {

    private final int radix;

    public RadixSort() {
        this(Byte.SIZE);
    }

    public RadixSort(int radix) {
        assert 32 % radix== 0;
        this.radix= radix;
    }

    public void sort(int[] data) {
        int[] histogram = new int[(1 << radix) + 1];
        int shift = 0;
        int mask = (1 << radix) - 1;
        int[] copy = new int[data.length];
        while (shift < Integer.SIZE) {
            Arrays.fill(histogram, 0);
            for (int i = 0; i < data.length; ++i) {
                ++histogram[((data[i] & mask) >> shift) + 1];
            }
            for (int i = 0; i < 1 << radix; ++i) {
                histogram[i + 1] += histogram[i];
            }
            for (int i = 0; i < data.length; ++i) {
                copy[histogram[(data[i] & mask) >> shift]++] = data[i];
            }
            for (int i = 0; i < data.length; ++i) {
                data[i] = copy[i];
            }
            shift += radix;
            mask <<= radix;
        }
    }
}

The time complexity is obviously linear, a temporary buffer is allocated, but in comparison to Arrays.sort it looks fairly spartan. Instinctively, cache locality looks fairly poor because the second inner loop of the three jumps all over the place. Will this implementation beat Arrays.sort (for integers)?

Benchmark

The algorithm is measured using arrays of random positive integers, for which both algorithms are equivalent, from a range of sizes. This isn’t always the best idea (the Tim Sort algorithm comes into its own on nearly sorted data), so take the result below with a pinch of salt. Care must be taken to copy the array in the benchmark since both algorithms are in-place.

public void launchBenchmark(String... jvmArgs) throws Exception {
        Options opt = new OptionsBuilder()
                .include(this.getClass().getName() + ".*")
                .mode(Mode.SampleTime)
                .mode(Mode.Throughput)
                .timeUnit(TimeUnit.MILLISECONDS)
                .measurementTime(TimeValue.seconds(10))
                .warmupIterations(10)
                .measurementIterations(10)
                .forks(1)
                .shouldFailOnError(true)
                .shouldDoGC(true)
                .jvmArgs(jvmArgs)
                .resultFormat(ResultFormatType.CSV) 
                .build();

        new Runner(opt).run();
    }


    @Benchmark
    public void Arrays_Sort(Data data, Blackhole bh) {
        int[] array = Arrays.copyOf(data.data, data.size);
        Arrays.sort(array);
        bh.consume(array);
    }

    @Benchmark
    public void Radix_Sort(Data data, Blackhole bh) {
        int[] array = Arrays.copyOf(data.data, data.size);
        data.radixSort.sort(array);
        bh.consume(array);
    }


    @State(Scope.Thread)
    public static class Data {

        @Param({"100", "1000", "10000", "100000", "1000000"})
        int size;

        int[] data;
        RadixSort radixSort = new RadixSort();

        @Setup(Level.Trial)
        public void init() {
            data = createArray(size);
        }
    }

    public static int[] createArray(int size) {
        int[] array = new int[size];
        Random random = new Random(0);
        for (int i = 0; i < size; ++i) {
            array[i] = Math.abs(random.nextInt());
        }
        return array;
    }
Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size
Arrays_Sort thrpt 1 10 1304.687189 147.380334 ops/ms 100
Arrays_Sort thrpt 1 10 78.518664 9.339994 ops/ms 1000
Arrays_Sort thrpt 1 10 1.700208 0.091836 ops/ms 10000
Arrays_Sort thrpt 1 10 0.133835 0.007146 ops/ms 100000
Arrays_Sort thrpt 1 10 0.010560 0.000409 ops/ms 1000000
Radix_Sort thrpt 1 10 404.807772 24.930898 ops/ms 100
Radix_Sort thrpt 1 10 51.787409 4.881181 ops/ms 1000
Radix_Sort thrpt 1 10 6.065590 0.576709 ops/ms 10000
Radix_Sort thrpt 1 10 0.620338 0.068776 ops/ms 100000
Radix_Sort thrpt 1 10 0.043098 0.004481 ops/ms 1000000
Arrays_Sort sample 1 3088586 0.000902 0.000018 ms/op 100
Arrays_Sort·p0.00 sample 1 1 0.000394 ms/op 100
Arrays_Sort·p0.50 sample 1 1 0.000790 ms/op 100
Arrays_Sort·p0.90 sample 1 1 0.000791 ms/op 100
Arrays_Sort·p0.95 sample 1 1 0.001186 ms/op 100
Arrays_Sort·p0.99 sample 1 1 0.001974 ms/op 100
Arrays_Sort·p0.999 sample 1 1 0.020128 ms/op 100
Arrays_Sort·p0.9999 sample 1 1 0.084096 ms/op 100
Arrays_Sort·p1.00 sample 1 1 4.096000 ms/op 100
Arrays_Sort sample 1 2127794 0.011876 0.000037 ms/op 1000
Arrays_Sort·p0.00 sample 1 1 0.007896 ms/op 1000
Arrays_Sort·p0.50 sample 1 1 0.009872 ms/op 1000
Arrays_Sort·p0.90 sample 1 1 0.015408 ms/op 1000
Arrays_Sort·p0.95 sample 1 1 0.024096 ms/op 1000
Arrays_Sort·p0.99 sample 1 1 0.033920 ms/op 1000
Arrays_Sort·p0.999 sample 1 1 0.061568 ms/op 1000
Arrays_Sort·p0.9999 sample 1 1 0.894976 ms/op 1000
Arrays_Sort·p1.00 sample 1 1 4.448256 ms/op 1000
Arrays_Sort sample 1 168991 0.591169 0.001671 ms/op 10000
Arrays_Sort·p0.00 sample 1 1 0.483840 ms/op 10000
Arrays_Sort·p0.50 sample 1 1 0.563200 ms/op 10000
Arrays_Sort·p0.90 sample 1 1 0.707584 ms/op 10000
Arrays_Sort·p0.95 sample 1 1 0.766976 ms/op 10000
Arrays_Sort·p0.99 sample 1 1 0.942080 ms/op 10000
Arrays_Sort·p0.999 sample 1 1 2.058273 ms/op 10000
Arrays_Sort·p0.9999 sample 1 1 7.526102 ms/op 10000
Arrays_Sort·p1.00 sample 1 1 46.333952 ms/op 10000
Arrays_Sort sample 1 13027 7.670135 0.021512 ms/op 100000
Arrays_Sort·p0.00 sample 1 1 6.356992 ms/op 100000
Arrays_Sort·p0.50 sample 1 1 7.634944 ms/op 100000
Arrays_Sort·p0.90 sample 1 1 8.454144 ms/op 100000
Arrays_Sort·p0.95 sample 1 1 8.742502 ms/op 100000
Arrays_Sort·p0.99 sample 1 1 9.666560 ms/op 100000
Arrays_Sort·p0.999 sample 1 1 12.916883 ms/op 100000
Arrays_Sort·p0.9999 sample 1 1 28.037900 ms/op 100000
Arrays_Sort·p1.00 sample 1 1 28.573696 ms/op 100000
Arrays_Sort sample 1 1042 96.278673 0.603645 ms/op 1000000
Arrays_Sort·p0.00 sample 1 1 86.114304 ms/op 1000000
Arrays_Sort·p0.50 sample 1 1 94.896128 ms/op 1000000
Arrays_Sort·p0.90 sample 1 1 104.293990 ms/op 1000000
Arrays_Sort·p0.95 sample 1 1 106.430464 ms/op 1000000
Arrays_Sort·p0.99 sample 1 1 111.223767 ms/op 1000000
Arrays_Sort·p0.999 sample 1 1 134.172770 ms/op 1000000
Arrays_Sort·p0.9999 sample 1 1 134.742016 ms/op 1000000
Arrays_Sort·p1.00 sample 1 1 134.742016 ms/op 1000000
Radix_Sort sample 1 2240042 0.002941 0.000033 ms/op 100
Radix_Sort·p0.00 sample 1 1 0.001578 ms/op 100
Radix_Sort·p0.50 sample 1 1 0.002368 ms/op 100
Radix_Sort·p0.90 sample 1 1 0.003556 ms/op 100
Radix_Sort·p0.95 sample 1 1 0.004344 ms/op 100
Radix_Sort·p0.99 sample 1 1 0.011056 ms/op 100
Radix_Sort·p0.999 sample 1 1 0.027232 ms/op 100
Radix_Sort·p0.9999 sample 1 1 0.731127 ms/op 100
Radix_Sort·p1.00 sample 1 1 5.660672 ms/op 100
Radix_Sort sample 1 2695825 0.018553 0.000038 ms/op 1000
Radix_Sort·p0.00 sample 1 1 0.013424 ms/op 1000
Radix_Sort·p0.50 sample 1 1 0.016576 ms/op 1000
Radix_Sort·p0.90 sample 1 1 0.025280 ms/op 1000
Radix_Sort·p0.95 sample 1 1 0.031200 ms/op 1000
Radix_Sort·p0.99 sample 1 1 0.050944 ms/op 1000
Radix_Sort·p0.999 sample 1 1 0.082944 ms/op 1000
Radix_Sort·p0.9999 sample 1 1 0.830295 ms/op 1000
Radix_Sort·p1.00 sample 1 1 6.660096 ms/op 1000
Radix_Sort sample 1 685589 0.145695 0.000234 ms/op 10000
Radix_Sort·p0.00 sample 1 1 0.112512 ms/op 10000
Radix_Sort·p0.50 sample 1 1 0.128000 ms/op 10000
Radix_Sort·p0.90 sample 1 1 0.196608 ms/op 10000
Radix_Sort·p0.95 sample 1 1 0.225792 ms/op 10000
Radix_Sort·p0.99 sample 1 1 0.309248 ms/op 10000
Radix_Sort·p0.999 sample 1 1 0.805888 ms/op 10000
Radix_Sort·p0.9999 sample 1 1 1.818141 ms/op 10000
Radix_Sort·p1.00 sample 1 1 14.401536 ms/op 10000
Radix_Sort sample 1 60843 1.641961 0.005783 ms/op 100000
Radix_Sort·p0.00 sample 1 1 1.251328 ms/op 100000
Radix_Sort·p0.50 sample 1 1 1.542144 ms/op 100000
Radix_Sort·p0.90 sample 1 1 2.002944 ms/op 100000
Radix_Sort·p0.95 sample 1 1 2.375680 ms/op 100000
Radix_Sort·p0.99 sample 1 1 3.447030 ms/op 100000
Radix_Sort·p0.999 sample 1 1 5.719294 ms/op 100000
Radix_Sort·p0.9999 sample 1 1 8.724165 ms/op 100000
Radix_Sort·p1.00 sample 1 1 13.074432 ms/op 100000
Radix_Sort sample 1 4846 20.640787 0.260926 ms/op 1000000
Radix_Sort·p0.00 sample 1 1 14.893056 ms/op 1000000
Radix_Sort·p0.50 sample 1 1 18.743296 ms/op 1000000
Radix_Sort·p0.90 sample 1 1 26.673152 ms/op 1000000
Radix_Sort·p0.95 sample 1 1 30.724915 ms/op 1000000
Radix_Sort·p0.99 sample 1 1 40.470446 ms/op 1000000
Radix_Sort·p0.999 sample 1 1 63.016600 ms/op 1000000
Radix_Sort·p0.9999 sample 1 1 136.052736 ms/op 1000000
Radix_Sort·p1.00 sample 1 1 136.052736 ms/op 1000000

The table tells an interesting story. Arrays.sort is vastly superior for small arrays (the arrays most people have), but for large arrays the custom implementation comes into its own. Interestingly, this is consistent with the computer science. If you need to sort large arrays of (unsigned) integers and care about performance, think about implementing radix sort.

Comparing Streams with Arrays

Java 8 added some great features which, combined with the vigilance of a responsible adult, make it easier to keep a Java code base civilised. Streams and lambdas, especially the limited support offered for primitive types, are a fantastic addition. Is there an observable performance cost to using these features? For a simple calculation amenable to instruction level parallelism, I compare modern and traditional implementations and observe the differences in instructions generated.

Sum of Squares

The sum of squares is the building block of a linear regression analysis so is ubiquitous in statistical computing. It is associative and therefore data-parallel. I compare four implementations: a sequential stream wrapping an array, a parallel stream wrapping an array, a generative sequential stream and a traditional for loop. The benchmark code is on github.

    @Benchmark
    public void SS_SequentialStream(Data state, Blackhole bh) {
        bh.consume(DoubleStream.of(state.data)
                               .map(x -> x * x)
                               .reduce((x, y) -> x + y));
    }

    @Benchmark
    public void SS_ParallelStream(Data state, Blackhole bh) {
        bh.consume(DoubleStream.of(state.data)
                               .parallel()
                               .map(x -> x * x)
                               .reduce((x, y) -> x + y));
    }

    @Benchmark
    public void SS_ForLoop(Data state, Blackhole bh) {
        double result = 0D;
        final double[] data = state.data;
        for (int i = 0; i < data.length; ++i) {
            result += data[i] * data[i];
        }
        bh.consume(result);
    }

    @Benchmark
    public void SS_GenerativeSequentialStream(Data state, Blackhole bh) {
        double[] data = state.data;
        bh.consume(IntStream.iterate(0, i -> i + 1)
                            .limit(1_000_000)
                            .mapToDouble(i -> data[i])
                            .map(x -> x * x)
                            .reduce((x, y) -> x + y));
    }

I must admit I prefer the readability of the stream versions, but let’s see if there is a comedown after the syntactic sugar rush.

Running a Benchmark

I compare the three implementations on an array of one million doubles. I am using JDK 9-ea, VM 9-ea+166 on a fairly powerful laptop with 8 processors:

$ cat /proc/cpuinfo
processor       : 0
vendor_id       : GenuineIntel
cpu family      : 6
model           : 94
model name      : Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz
stepping        : 3
cpu MHz         : 2592.000
cache size      : 256 KB
physical id     : 0
siblings        : 8
core id         : 0
cpu cores       : 4
apicid          : 0
initial apicid  : 0
fpu             : yes
fpu_exception   : yes
cpuid level     : 22
wp              : yes
flags           : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe pni dtes64 monitor ds_cpl vmx est tm2 ssse3 fma cx16 xtpr pdcm sse4_1 sse4_2 x2apic movbe popcnt aes xsave osxsave avx f16c rdrand lahf_lm ida arat epb xsaveopt pln pts dtherm fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm mpx rdseed adx smap clflushopt
clflush size    : 64
cache_alignment : 64
address sizes   : 39 bits physical, 48 bits virtual
power management:

processor       : 1
vendor_id       : GenuineIntel
cpu family      : 6
model           : 94
model name      : Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz
stepping        : 3
cpu MHz         : 2592.000
cache size      : 256 KB
physical id     : 0
siblings        : 8
core id         : 0
cpu cores       : 4
apicid          : 1
initial apicid  : 1
fpu             : yes
fpu_exception   : yes
cpuid level     : 22
wp              : yes
flags           : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe pni dtes64 monitor ds_cpl vmx est tm2 ssse3 fma cx16 xtpr pdcm sse4_1 sse4_2 x2apic movbe popcnt aes xsave osxsave avx f16c rdrand lahf_lm ida arat epb xsaveopt pln pts dtherm fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm mpx rdseed adx smap clflushopt
clflush size    : 64
cache_alignment : 64
address sizes   : 39 bits physical, 48 bits virtual
power management:

processor       : 2
vendor_id       : GenuineIntel
cpu family      : 6
model           : 94
model name      : Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz
stepping        : 3
cpu MHz         : 2592.000
cache size      : 256 KB
physical id     : 0
siblings        : 8
core id         : 1
cpu cores       : 4
apicid          : 2
initial apicid  : 2
fpu             : yes
fpu_exception   : yes
cpuid level     : 22
wp              : yes
flags           : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe pni dtes64 monitor ds_cpl vmx est tm2 ssse3 fma cx16 xtpr pdcm sse4_1 sse4_2 x2apic movbe popcnt aes xsave osxsave avx f16c rdrand lahf_lm ida arat epb xsaveopt pln pts dtherm fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm mpx rdseed adx smap clflushopt
clflush size    : 64
cache_alignment : 64
address sizes   : 39 bits physical, 48 bits virtual
power management:

processor       : 3
vendor_id       : GenuineIntel
cpu family      : 6
model           : 94
model name      : Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz
stepping        : 3
cpu MHz         : 2592.000
cache size      : 256 KB
physical id     : 0
siblings        : 8
core id         : 1
cpu cores       : 4
apicid          : 3
initial apicid  : 3
fpu             : yes
fpu_exception   : yes
cpuid level     : 22
wp              : yes
flags           : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe pni dtes64 monitor ds_cpl vmx est tm2 ssse3 fma cx16 xtpr pdcm sse4_1 sse4_2 x2apic movbe popcnt aes xsave osxsave avx f16c rdrand lahf_lm ida arat epb xsaveopt pln pts dtherm fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm mpx rdseed adx smap clflushopt
clflush size    : 64
cache_alignment : 64
address sizes   : 39 bits physical, 48 bits virtual
power management:

processor       : 4
vendor_id       : GenuineIntel
cpu family      : 6
model           : 94
model name      : Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz
stepping        : 3
cpu MHz         : 2592.000
cache size      : 256 KB
physical id     : 0
siblings        : 8
core id         : 2
cpu cores       : 4
apicid          : 4
initial apicid  : 4
fpu             : yes
fpu_exception   : yes
cpuid level     : 22
wp              : yes
flags           : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe pni dtes64 monitor ds_cpl vmx est tm2 ssse3 fma cx16 xtpr pdcm sse4_1 sse4_2 x2apic movbe popcnt aes xsave osxsave avx f16c rdrand lahf_lm ida arat epb xsaveopt pln pts dtherm fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm mpx rdseed adx smap clflushopt
clflush size    : 64
cache_alignment : 64
address sizes   : 39 bits physical, 48 bits virtual
power management:

processor       : 5
vendor_id       : GenuineIntel
cpu family      : 6
model           : 94
model name      : Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz
stepping        : 3
cpu MHz         : 2592.000
cache size      : 256 KB
physical id     : 0
siblings        : 8
core id         : 2
cpu cores       : 4
apicid          : 5
initial apicid  : 5
fpu             : yes
fpu_exception   : yes
cpuid level     : 22
wp              : yes
flags           : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe pni dtes64 monitor ds_cpl vmx est tm2 ssse3 fma cx16 xtpr pdcm sse4_1 sse4_2 x2apic movbe popcnt aes xsave osxsave avx f16c rdrand lahf_lm ida arat epb xsaveopt pln pts dtherm fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm mpx rdseed adx smap clflushopt
clflush size    : 64
cache_alignment : 64
address sizes   : 39 bits physical, 48 bits virtual
power management:

processor       : 6
vendor_id       : GenuineIntel
cpu family      : 6
model           : 94
model name      : Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz
stepping        : 3
cpu MHz         : 2592.000
cache size      : 256 KB
physical id     : 0
siblings        : 8
core id         : 3
cpu cores       : 4
apicid          : 6
initial apicid  : 6
fpu             : yes
fpu_exception   : yes
cpuid level     : 22
wp              : yes
flags           : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe pni dtes64 monitor ds_cpl vmx est tm2 ssse3 fma cx16 xtpr pdcm sse4_1 sse4_2 x2apic movbe popcnt aes xsave osxsave avx f16c rdrand lahf_lm ida arat epb xsaveopt pln pts dtherm fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm mpx rdseed adx smap clflushopt
clflush size    : 64
cache_alignment : 64
address sizes   : 39 bits physical, 48 bits virtual
power management:

processor       : 7
vendor_id       : GenuineIntel
cpu family      : 6
model           : 94
model name      : Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz
stepping        : 3
cpu MHz         : 2592.000
cache size      : 256 KB
physical id     : 0
siblings        : 8
core id         : 3
cpu cores       : 4
apicid          : 7
initial apicid  : 7
fpu             : yes
fpu_exception   : yes
cpuid level     : 22
wp              : yes
flags           : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe pni dtes64 monitor ds_cpl vmx est tm2 ssse3 fma cx16 xtpr pdcm sse4_1 sse4_2 x2apic movbe popcnt aes xsave osxsave avx f16c rdrand lahf_lm ida arat epb xsaveopt pln pts dtherm fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm mpx rdseed adx smap clflushopt
clflush size    : 64
cache_alignment : 64
address sizes   : 39 bits physical, 48 bits virtual
power management:

Before running the benchmark we might expect the for loop and stream to have similar performance, and the parallel version to be about eight times faster. The generative version is very similar to the for loop so a slow down might not be expected. To isolate vectorised optimisations, I first run with SuperWord disabled.

Benchmark Mode Count Score Units
SS_ForLoop thrpt 10 0.540 ± 0.072 ops/ms
SS_ParallelStream thrpt 10 2.336 ± 0.412 ops/ms
SS_SequentialStream thrpt 10 0.565 ± 0.052 ops/ms
SS_GenerativeSequentialStream thrpt 10 0.247 ± 0.036 ops/ms

The for loop and stream are similar but the stream is better. The parallel version is faster but is either not utilising all of the cores or incurring a cost; the improvement in throughput is only a factor of 4.3. It has never felt like a good idea to execute code blindly on a default fork join pool without safe guards. There is a huge difference between streams which wrap arrays and a stream which iterates over an array, most likely related to cache locality.

What happens when SuperWord optimisations are enabled? Throughput improves except for the generative stream (suggesting no SIMD execution at all), but the for loop improves the most.

Benchmark Mode Count Score Units
SS_ForLoop thrpt 10 0.614 ± 0.035 ops/ms
SS_ParallelStream thrpt 10 2.551 ± 0.095 ops/ms
SS_SequentialStream thrpt 10 0.596 ± 0.043 ops/ms
SS_GenerativeSequentialStream thrpt 10 0.240 ± 0.049 ops/ms

The three nines sample is actually worse than the for loop in the parallel stream version, whereas the generative stream is off the chart.

Benchmark Mode Score Units
SS_ForLoop·p0.999 sample 5.424 ms/op
SS_ParallelStream·p0.999 sample 6.077 ms/op
SS_SequentialStream·p0.999 sample 6.340 ms/op
SS_GenerativeSequentialStream·p0.999 sample 14.516 ms/op

Inspecting the assembly code generated, it is clear that the for loop body is actually only compiling down to SSE here: it’s only loading one double at a time into registers.

0x000001aada783fca: vmovsd  xmm0,qword ptr [rbp+r13*8+10h] ...
0x000001aada783fd1: vmulsd  xmm0,xmm0,xmm0 ...
0x000001aada783fd5: vmovq   xmm1,rbx ...
0x000001aada783fda: vaddsd  xmm0,xmm1,xmm0 ...

But so is the sequential stream – here is the lambda implementing the square:

0x00000170c11228be: vmulsd  xmm1,xmm0,xmm0    ;...c5
                                                ;...fb
                                                ;...59
                                                ;...c8
                                                ;*dmul {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.mythbusters.streams.StreamReduction::lambda$SS_SequentialStream$0@2 (line 66)
                                                ; - com.openkappa.mythbusters.streams.StreamReduction$$Lambda$40/1123395594::applyAsDouble@1
                                                ; - java.util.stream.DoublePipeline$2$1::accept@12 (line 213)
                                                ; - java.util.Spliterators$DoubleArraySpliterator::forEachRemaining@53 (line 1198)
                                                ; - java.util.Spliterator$OfDouble::forEachRemaining@12 (line 828)
                                                ; - java.util.stream.AbstractPipeline::copyInto@32 (line 484)
                                                ; - java.util.stream.AbstractPipeline::wrapAndCopyInto@13 (line 474)
                                                ; - java.util.stream.ReduceOps$ReduceOp::evaluateSequential@6 (line 913)
                                                ; - java.util.stream.AbstractPipeline::evaluate@88 (line 234)

And here is the reduce:

  0x00000170c1118610: vaddsd  xmm1,xmm1,xmm2    ;...c5
                                                ;...f3
                                                ;...58
                                                ;...ca
                                                ;*dadd {reexecute=0 rethrow=0 return_oop=0}
                                                ; - com.openkappa.mythbusters.streams.StreamReduction::lambda$SS_SequentialStream$1@2 (line 66)
                                                ; - com.openkappa.mythbusters.streams.StreamReduction$$Lambda$41/1355169748::applyAsDouble@2
                                                ; - java.util.stream.ReduceOps$12ReducingSink::accept@30 (line 693)
                                                ; - java.util.stream.DoublePipeline$2$1::accept@17 (line 213)
                                                ; - java.util.Spliterators$DoubleArraySpliterator::forEachRemaining@53 (line 1198)

Is it mechanically sympathetic to separate these instructions like this? Probably not, which explains the degradation in performance relative to the for loop, but it’s not so drastic that a completely different instruction set is being used. Use of the same SSE instructions can be seen with Stream.parallelStream.

Conclusion

The same optimisations are applied to streams when they are countable. Arrays seems to do slightly better, but not well enough to eschew streams. Uncountable or opaquely countable streams perform significantly worse than countable streams or traditional code. Stream.parallelStream can improve performance but does not make the best use of the cores available and using a global thread pool leads to unpredictability. Balancing performance and readability leads to an approach where data is stored in paginated arrays and accessed and transformed via the streams API.