<div dir="ltr"><div class="gmail_quote"><div dir="ltr">On Sat, 4 Jul 2015 at 02:12 Jason Swails <<a href="mailto:jason.swails@gmail.com">jason.swails@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Fri, Jul 3, 2015 at 11:13 AM, Oscar Benjamin <span dir="ltr"><<a href="mailto:oscar.j.benjamin@gmail.com" target="_blank">oscar.j.benjamin@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><span>On 2 July 2015 at 18:29, Jason Swails <<a href="mailto:jason.swails@gmail.com" target="_blank">jason.swails@gmail.com</a>> wrote:</span><br>
<br>
Where is the 32 bit one looks like:<br>
<br>
$ objdump -d a.out.32 | less<br>
... </blockquote></div></div></div></blockquote><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
804843e: fildl -0x14(%ebp)<br>
8048441: fmull -0x10(%ebp)<br>
8048444: fnstcw -0x1a(%ebp)<br>
8048447: movzwl -0x1a(%ebp),%eax<br>
804844b: mov $0xc,%ah<br>
804844d: mov %ax,-0x1c(%ebp)<br>
8048451: fldcw -0x1c(%ebp)<br>
8048454: fistpl -0x20(%ebp)<br>
8048457: fldcw -0x1a(%ebp)<br>
804845a: mov -0x20(%ebp),%eax<br>
804845d: cmp -0x14(%ebp),%eax<br>
8048460: jne 8048477 <main+0x5c><br>
8048462: sub $0x8,%esp<br>
8048465: pushl -0x14(%ebp)<br>
8048468: push $0x8048520<br>
804846d: call 80482f0 <printf@plt><br>
8048472: add $0x10,%esp<br>
8048475: jmp 8048484 <main+0x69><br>
8048477: addl $0x1,-0x14(%ebp)<br>
804847b: cmpl $0xf423f,-0x14(%ebp)<br>
8048482: jle 804843e <main+0x23><br>
...<br>
<br>
So the 64 bit one is using SSE instructions and the 32-bit one is<br>
using x87. That could explain the difference you see at the C level<br>
but I don't see it on this CPU (/proc/cpuinfo says Intel(R) Core(TM)<br>
i5-3427U CPU @ 1.80GHz).<br></blockquote><div><br></div></div></div></div><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div><div class="gmail_default" style="color:rgb(0,0,0);display:inline">Hmm. Well that could explain why you don't get the same results as me. My CPU is a </div>AMD FX(tm)-6100 Six-Core Processor<div class="gmail_default" style="color:rgb(0,0,0);display:inline"> (from /proc/cpuinfo). My objdump looks the same as yours for the 64-bit version, but for 32-bit it looks like:</div></div></div></div></div></blockquote><div><br></div><div>So if we have different generated machine instructions it suggests a difference in the way it was compiled rather than in the hardware itself. (Although it could be that the compilers were changed because the hardware was inconsistent in this particular usage).<br></div><div><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><br><div><div class="gmail_default" style="color:rgb(0,0,0)">However, I have no experience looking at raw assembler, so I can't discern what it is I'm even looking at (nor do I know what explicit SSE instructions look like in assembler).</div></div></div></div></div></blockquote><div><br></div><div>The give away is that SSE instructions use the XMM registers so where you see %xmm0 etc data is being loaded ready for SSE instructions. I'll translate the important part of the 32 bit code below:<br></div><div><div><div class="gmail_default" style="color:rgb(0,0,0);display:inline"><br></div></div><div><div class="gmail_default" style="color:rgb(0,0,0);display:inline">...<br></div></div>
804843a: db 44 24 14 fildl 0x14(%esp)<div class="gmail_default" style="color:rgb(0,0,0);display:inline"></div><br><br></div><div>FILDL - F is for float and indicates that this is an x87 instruction. LD is for load and indicates that we want to load from memory to the x87 register stack. The I is for integer: we're loading from int format. I can't remember what the L suffix means.<br><br></div><div>This instruction will load the integer from the memory stack at offset 0x14 from the stack pointer (esp) pushing it onto the x87 stack. Since the x87 stack uses 80-bit extended-real precision format internally this instruction performs a format conversion. The conversion should be lossless as 80-bit FP uses a 64 bit mantissa so it should be able to represent any 64-bit signed or unsigned integer (we're only using 32 bits in this example).<br></div><div><br><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"> 804843e: dc 4c 24 18 fmull 0x18(%esp)<br><br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)">FMULL (x87 multiply) - This multiplies the top of the x87 stack (the number we previously stored) by the 64 bit float at 0x18 stack offset. In context that's the number 1-.5**53 which was presumably calculated earlier. The top of the x87 stack is overwritten with the result. This calculation is performed using 80-bit floating point with a 64 bit mantissa.<br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"><br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"> 8048442: dd 5c 24 08 fstpl 0x8(%esp)<br><br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)">FSTPL (x87 store and pop) - Pop the top of the x87 stack pushing to the memory stack. This is storing the result of the multiplication back to RAM (0x8 from stack) as a 64 bit float. Since it was internally stored as an 80-bit float there is a format conversion here that can lose precision. Rounding can occur according to the settings in the x87 control word.<br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"><br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"> 8048446: f2 0f 2c 44 24 08 cvttsd2si 0x8(%esp),%eax<br><br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)">This converts the 64 bit float stored at 0x8 from stack to a signed int, truncating if inexact, and stores the result in register eax. This operation is guaranteed to truncate (rather than floor etc) if it needs to round.<br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"><br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"> 804844c: 3b 44 24 14 cmp 0x14(%esp),%eax<br><br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)">This compares the result of the above processing (now stored in eax) with the int we started with (still at offset 0x14 from stack). The cmp instruction sets flags ready for the jne instruction below.<br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"><br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"> 8048450: 75 16 jne 8048468 <main+0x4b><br><br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)">JNE (jump if not equal) - this is the branching instruction. If the previously compared values are unequal we skip forward to 8048468 otherwise we execute the instructions below which call printf and jmp (break) out of the loop.<br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"><br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"> 8048452: 8b 44 24 14 mov 0x14(%esp),%eax</span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"> 8048456: 89 44 24 04 mov %eax,0x4(%esp)</span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"> 804845a: c7 04 24 10 85 04 08 movl $0x8048510,(%esp)</span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"> 8048461: e8 8a fe ff ff call 80482f0 <printf@plt></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"> 8048466: eb 0f jmp 8048477 <main+0x5a><br><br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)">This is where the jne goes. We increment the int and jump (jle) back to the top of the loop (unless the loop is finished).<br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"><br></span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"> 8048468: 83 44 24 14 01 addl $0x1,0x14(%esp)</span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"> 804846d: 81 7c 24 14 3f 42 0f cmpl $0xf423f,0x14(%esp)</span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"> 8048474: 00 </span></div><div class="gmail_default" style="color:rgb(0,0,0)"><span style="color:rgb(34,34,34)"> 8048475: 7e c3 jle 804843a <main+0x1d></span></div><div><div class="gmail_default" style="color:rgb(0,0,0);display:inline">...</div><div class="gmail_default" style="color:rgb(0,0,0);display:inline"></div> <br><br></div><div>My suspicion falls on the FSTPL instruction. The rounding behaviour there depends on the rounding flags set in the x87 control word which have not been explicitly controlled here. The calculation we want to perform cannot be performed exactly in 64 bit floating point:<br><br>>>> from fractions import Fraction as F<br>>>> x = 1-.5**53<br>>>> F(x)*F(2049) == F(x*2049)<br>False<br><br><div>This is because the true result cannot be represented as a 64 bit float. Instead the result comes out to the nearest available float:<br></div><br>>>> float(F(x)*F(2049)) == x*2049<br>True<br><br></div><div>This means that the x87 register will be storing a higher precision result in its 80 bit format. This result will have to be rounded by the FSTPL instruction.<br><br>If you look at the assembly output I showed you'll see the instructions FNSTCW/FLDCW (x87 store/load control word) which are used to manipulate the control word to tell the FPU how to perform this kind of rounding. The fact that we don't see it in your compiled output could indicate a C compiler bug which could in turn explain the different behaviour people see from Python.<br><br></div><div>To understand exactly why 2049 is the number where it fails consider that it is the smallest integer that requires 12 bits of mantissa in floating point format. The number 1-.5**53 has a mantissa that is 53 ones:<br><br>>>> x = 1-.5**53<br>>>> x.hex()<br>'0x1.fffffffffffffp-1'<br><br>When extended to 80 bit real-extended format with a 64 bit mantissa it will have 11 trailing zeros. So I think multiplication of x with any integer less than 2049 can be performed exactly by the FMULL instruction. I haven't fully considered what impact that would have but it seems reasonable that this is why 2049 is the first number that fails.<br></div><div><br>--<br></div><div>Oscar<br></div></div></div></div>