On Jun 8, 2018, at 5:47 AM, KHMan <keinhong@gmail.com> wrote:
Operation / Op_Count / Diff_Count
=================================
A+B 1,000,000,000 0
A-B 1,000,000,000 0
A*B 1,000,000,000 243,584
A/B 1,000,000,000 244,107
For A*B and A/B, this is roughly 1 in 4000.
All different results differ by 1 ULP. Nothing more than 1 ULP was seen.
I found another blog doing the same test, error also about 1:4000
Just some math behind the observation:
Double rounding occurs when both rounded in the same direction.
Difference between 53-bits and 64-bits are 11-bits.
--> chances of same direction rounding is about 50%
--> chances of double rounding = 50% / 2^11 = 1 / 4096
With just 1 op, of course the error is at most 1 ULP
Operation : *
n1 = 0.18987016622488206 0x1.84daa65360288p-3
n2 = 0.98339548253950804 0x1.f77f9cd91528bp-1
Here are the hex representations of the significands:
fpu8087 = 7e65b9c2a810e
binary128 = 7E65B9C2A810E5EAB7A33195161D
fpu64bit = 7e65b9c2a810f
Using this comparison, it appears that the binary128 result is nearer to fpu8087 than to fpu64bit. We can also verify this by getting a more accurate decimal version of the 64-bit binary representation (done by manual fiddling, so not definitive):
You do realize we have to work in 53-bits for comparison.
Getting more precise result with more bits does not show anything.
My guess is you set (n1, n2) as decimal string -> __float128.
Had you use the actual hexfloats, you should get this:
n1 * n2
= 0x1.84daa65360288p-3 * 0x1.f77f9cd91528bp-1
= 0x1.7e65b9c2a810e 800eb4c1577ec p-3
--> value should round-up, to 0x1.7e65b9c2a810fp-3
--> 53-bits rounding is doing the best it can with 2 doubles.