Why R Math Functions on Windows are Slow, and How to Fix It

R on windows has much slower versions of the log, sine and cosine functions than are available on other platforms, and this can be a serious performance bottleneck for programs which frequently call these math functions.  The reason for this is that the library R uses to obtain the log function on windows (libmingwex.a) contains a version of the log implementation which is out of date relative to more modern code and much slower than other available versions.  That the glibc implementations of common math functions are slow is a known issue that others have discussed on the internet.  This post uses the log function as an example to show specifically why it is slow, and then suggest some quick work arounds for these math functions.

The log Function on MinGW / Windows

Comparing the assembly code of the log function as generated with MinGW on windows with the assembly generated for the log function on Mac OSX shows why the code is slow. Below is the assembly code for log function as decompiled from MinGW (their implementation is basically a cut/paste from GNU libc). The  crucial feature here is that most instructions start with an “f” indicating that they are using the floating point registers, and there is one instruction the fyl2xp1 that is one of the most expensive operations out there. This instruction takes a log using hardware, and is known to be slower than most other possible ways to calculate log for modern machines.

0000000000000010 <__logl_internal>:
  10:	d9 ed                	fldln2 
  12:	db 2a                	fldt   (%rdx)
  14:	d9 c0                	fld    %st(0)
  16:	dc 25 e4 ff ff ff    	fsubl  -0x1c(%rip)        # 0 <one>
  1c:	d9 c0                	fld    %st(0)
  1e:	d9 e1                	fabs   
  20:	dc 1d e2 ff ff ff    	fcompl -0x1e(%rip)        # 8 <limit>
  26:	df e0                	fnstsw %ax
  28:	80 e4 45             	and    $0x45,%ah
  2b:	74 12                	je     3f <__logl_internal+0x2f>
  2d:	dd d9                	fstp   %st(1)
  2f:	d9 f9                	fyl2xp1 
  31:	48 89 c8             	mov    %rcx,%rax
  34:	48 c7 41 08 00 00 00 	movq   $0x0,0x8(%rcx)
  3b:	00 
  3c:	db 39                	fstpt  (%rcx)
  3e:	c3                   	retq   
  3f:	dd d8                	fstp   %st(0)
  41:	d9 f1                	fyl2x  
  43:	48 89 c8             	mov    %rcx,%rax
  46:	48 c7 41 08 00 00 00 	movq   $0x0,0x8(%rcx)
  4d:	00 
  4e:	db 39                	fstpt  (%rcx)

These floating point register instructions were state of the art way back when, but nowadays most hardware is using special XMM/YMM registers and instructions that can make for much faster code to calculate logarithms (See Section 17 for more information on this).

The log Function on Mac OSX

Demonstrating faster assembly code is the Mac OSX implementation of log shown below.  The key feature here is that it uses XMM registers and has a more modern and performant implementation.  Note that many other implementations (e.g. the Microsoft compiler) also use faster versions like this. That is, Windows is not slow, R is slow on Windows. It’s all the same Intel chips under the hood.

+0x00	vmovq               %xmm0, %rax
+0x05	shrq                $32, %rax
+0x09	movl                %eax, %edx
+0x0b	subl                $1048576, %eax
+0x10	cmpl                $2145386496, %eax
+0x15	jae                 "0x7fff9793e100+0xe4"
+0x1b	    andl                $1048575, %eax
+0x20	    addl                $4096, %eax
+0x25	    andl                $2088960, %eax
+0x2a	    shrl                $9, %eax
+0x2d	    addl                $3222802432, %edx
+0x33	    sarl                $20, %edx
+0x36	    vcvtsi2sdl          %edx, %xmm1, %xmm1
+0x3a	    vmovsd              112942(%rip), %xmm6
+0x42	    vandpd              %xmm6, %xmm0, %xmm0
+0x46	    vmovsd              74890(%rip), %xmm7
+0x4e	    vorpd               %xmm7, %xmm0, %xmm0
+0x52	    leaq                112999(%rip), %rdx
+0x59	    vmovss              (%rdx,%rax), %xmm3
+0x5e	    vpsllq              $32, %xmm3, %xmm3
+0x63	    vfmsub213sd         %xmm7, %xmm0, %xmm3
+0x68	    vmovddup            %xmm3, %xmm3
+0x6c	    vaddpd              -48(%rdx), %xmm3, %xmm4
+0x71	    vfmadd213pd         -32(%rdx), %xmm3, %xmm4
+0x77	    vmovapd             -16(%rdx), %xmm5
+0x7c	    vaddsd              %xmm3, %xmm5, %xmm5
+0x80	    vmulpd              %xmm3, %xmm5, %xmm5
+0x84	    vmulpd              %xmm5, %xmm4, %xmm4
+0x88	    vmovhlps            %xmm4, %xmm5, %xmm5
+0x8c	    vmulsd              %xmm5, %xmm4, %xmm4
+0x90	    je                  "0x7fff9793e100+0xc7"
+0x92	    vmovss              4(%rdx,%rax), %xmm5
+0x98	    vpsllq              $32, %xmm5, %xmm5
+0x9d	    vmulsd              112867(%rip), %xmm1, %xmm0
+0xa5	    vmulsd              112851(%rip), %xmm1, %xmm1
+0xad	    vaddsd              8(%rdx,%rax), %xmm0, %xmm0
+0xb3	    vaddsd              %xmm5, %xmm1, %xmm1
+0xb7	    vaddsd              %xmm4, %xmm0, %xmm0
+0xbb	    vaddsd              %xmm3, %xmm0, %xmm0
+0xbf	    vaddsd              %xmm1, %xmm0, %xmm0
+0xc3	    vzeroupper
+0xc6	    retq
+0xc7	    vmovss              4(%rdx,%rax), %xmm5
+0xcd	    vpsllq              $32, %xmm5, %xmm5
+0xd2	    vaddsd              %xmm5, %xmm3, %xmm3
+0xd6	    vaddsd              8(%rdx,%rax), %xmm4, %xmm0
+0xdc	    vaddsd              %xmm3, %xmm0, %xmm0
+0xe0	    vzeroupper
+0xe3	    retq

Solving the problem

Ideally we could just push a faster log implementation to the R core repo for Windows, but the fact of the matter is the R core team does not accept code changes that only improve performance (and fair play to them, R is the most cross platform software out their, and that wasn’t easy to do).  Also, realistically the current implementation is feasible for most peoples work. So the solution is to switch out the log function for a better one only when it matters, and if it matters that means most of the library is already written in compiled code. Where to get a cross-platform compatible log function? It seems the folks over at Julia ran into a similar problem with slow Log on different machines. They looked at two version in libm (e.g.  Libm Version #1 and Libm Version #2), and also some really crazy implementations such as this one available from Yepp. However, what they settled on is a very good function, and we can port those over directly into R.  An example is shown in the pull request here.