Fix fast_cast int64/uint64→double under -ffast-math

With -ffast-math (-fassociative-math), the compiler may reassociate
(xH-C)+xL into xH+(xL-C). Since xL<<C this causes catastrophic
cancellation, zeroing the result.

Add an asm register barrier between the sub and add to prevent
reassociation. The barrier forces the intermediate into a register
without emitting any instructions.

Fixes #1264

Author: DiamonDinoia

include/xsimd/arch/xsimd_avx2.hpp

@@ -528,6 +528,13 @@ namespace xsimd
                                                  0xFFFF, 0xFFFF, 0x0000, 0x0000, 0xFFFF, 0xFFFF, 0x0000, 0x0000);
                 __m256i xL = _mm256_or_si256(_mm256_and_si256(mask, x), _mm256_andnot_si256(mask, _mm256_castpd_si256(_mm256_set1_pd(0x0010000000000000)))); //  2^52
                 __m256d f = _mm256_sub_pd(_mm256_castsi256_pd(xH), _mm256_set1_pd(19342813118337666422669312.)); //  2^84 + 2^52
+                // With -ffast-math, the compiler may reassociate (xH-C)+xL into
+                // xH+(xL-C). Since xL<<C this causes catastrophic cancellation.
+                // The asm barrier forces f into a register before the add, blocking
+                // the reorder. It emits zero instructions.
+#if defined(__GNUC__)
+                __asm__ volatile("" : "+x"(f));
+#endif
                 return _mm256_add_pd(f, _mm256_castsi256_pd(xL));
             }
 
@@ -543,6 +550,10 @@ namespace xsimd
                                                  0xFFFF, 0xFFFF, 0xFFFF, 0x0000, 0xFFFF, 0xFFFF, 0xFFFF, 0x0000);
                 __m256i xL = _mm256_or_si256(_mm256_and_si256(mask, x), _mm256_andnot_si256(mask, _mm256_castpd_si256(_mm256_set1_pd(0x0010000000000000)))); //  2^52
                 __m256d f = _mm256_sub_pd(_mm256_castsi256_pd(xH), _mm256_set1_pd(442726361368656609280.)); //  3*2^67 + 2^52
+                // See above: prevent -ffast-math from reassociating (xH-C)+xL.
+#if defined(__GNUC__)
+                __asm__ volatile("" : "+x"(f));
+#endif
                 return _mm256_add_pd(f, _mm256_castsi256_pd(xL));
             }
         }

include/xsimd/arch/xsimd_sse4_1.hpp

@@ -53,6 +53,13 @@ namespace xsimd
                 xH = _mm_add_epi64(xH, _mm_castpd_si128(_mm_set1_pd(442721857769029238784.))); //  3*2^67
                 __m128i xL = _mm_blend_epi16(x, _mm_castpd_si128(_mm_set1_pd(0x0010000000000000)), 0x88); //  2^52
                 __m128d f = _mm_sub_pd(_mm_castsi128_pd(xH), _mm_set1_pd(442726361368656609280.)); //  3*2^67 + 2^52
+                // With -ffast-math, the compiler may reassociate (xH-C)+xL into
+                // xH+(xL-C). Since xL<<C this causes catastrophic cancellation.
+                // The asm barrier forces f into a register before the add, blocking
+                // the reorder. It emits zero instructions.
+#if defined(__GNUC__)
+                __asm__ volatile("" : "+x"(f));
+#endif
                 return _mm_add_pd(f, _mm_castsi128_pd(xL));
             }
 
@@ -64,6 +71,10 @@ namespace xsimd
                 xH = _mm_or_si128(xH, _mm_castpd_si128(_mm_set1_pd(19342813113834066795298816.))); //  2^84
                 __m128i xL = _mm_blend_epi16(x, _mm_castpd_si128(_mm_set1_pd(0x0010000000000000)), 0xcc); //  2^52
                 __m128d f = _mm_sub_pd(_mm_castsi128_pd(xH), _mm_set1_pd(19342813118337666422669312.)); //  2^84 + 2^52
+                // See above: prevent -ffast-math from reassociating (xH-C)+xL.
+#if defined(__GNUC__)
+                __asm__ volatile("" : "+x"(f));
+#endif
                 return _mm_add_pd(f, _mm_castsi128_pd(xL));
             }
         }