Remove undefined behaviour in transpose_and_copy.

There was a -DUBSAN guard for use with undefined behaviour sanitizers,
but actually the only code change that's necessary is to change the
explicit writes to a memcpy.  This fixes unoptimised zig-cc builds
which use -fsanitize=undefined by default.

I've tested gcc7, gcc13, clang10, clang16, zig cc, icc and icl.  All
seem fine with this change.

NB: Zig cc is the only on to vectorise the most naive implementation,
and it does so very well, basically creating an SSE4 version of the
loop. That is the even fastest implementation for zig cc.

Author: jkbonfield

htscodecs/rANS_static32x16pr_sse4.c

@@ -807,19 +807,59 @@ unsigned char *rans_uncompress_O0_32x16_sse4(unsigned char *in,
 static inline void transpose_and_copy(uint8_t *out, int iN[32],
                                       uint8_t t[32][32]) {
     int z;
-#ifdef UBSAN
-    // Simplified version to avoid undefined behaviour sanitiser warnings.
+
+    // Simplified version from below.
+    // This is pretty good with zig cc, but slow on clang and very slow on
+    // gcc, even with -O3
+    /*
     for (z = 0; z < NX; z++) {
         int k;
         for (k = 0; k < 32; k++)
             out[iN[z]+k] = t[k][z];
         iN[z] += 32;
     }
-#else
-    // Unaligned access.  We know we can get away with this as this
-    // code is only ever executed on x86 platforms which permit this.
-    for (z = 0; z < NX; z+=4) {
-        *(uint64_t *)&out[iN[z]] =
+    */
+
+
+    // A better approach for clang and gcc can be had with some manual
+    // restructuring to attempt to do the two loops in explcit blocks.
+    // With gcc -O3 or -O2 -ftree-vectorize this is quite fast, as is clang
+    // and zig but neither beat the version below (or, for zig, the basic
+    // code above).
+    //
+    // It's left here incase we ever want to move to tidier code and
+    // to understand what auto-vectorises and what doesn't.
+    /*
+#define NZ 2
+#define NK 8
+    for (z = 0; z < 32; z+=NZ) {
+        for (int k = 0; k < 32; k+=NK) {
+            for (int z0 = z; z0 < z+NZ; z0++) {
+                uint8_t tmp[NK];// __attribute__((aligned(32)));
+                //uint8_t (*RESTRICT t0)[32] = &t[k];
+                for (int k0 = 0; k0 < NK; k0++)
+                    //tmp[k0] = t0[k0][z0];
+                    tmp[k0] = t[k+k0][z0];
+                memcpy(&out[iN[z0]+k], tmp, NK);
+            }
+        }
+        for (int z0 = z; z0 < z+NZ; z0++)
+            iN[z0] += 32;
+    }
+    */
+
+    // Manually unrolled code.
+    // This is fastest on gcc and clang and not far behind with zig cc.
+    // It also doesn't need aggressive gcc optimisation levels to be
+    // efficient.
+    //
+    // It works by constructing 64-bit ints and copying them with a single
+    // memory write.  The fixed size memcpys just boil down to a memory write,
+    // but unlike the earlier versions that did this direct, this isn't
+    // exploiting undefined behaviour.
+    for (z = 0; z < 32; z+=4) {
+        uint64_t i64;
+        i64 =
             ((uint64_t)(t[0][z])<< 0) +
             ((uint64_t)(t[1][z])<< 8) +
             ((uint64_t)(t[2][z])<<16) +
@@ -828,7 +868,8 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[5][z])<<40) +
             ((uint64_t)(t[6][z])<<48) +
             ((uint64_t)(t[7][z])<<56);
-        *(uint64_t *)&out[iN[z+1]] =
+        memcpy(&out[iN[z]], &i64, 8);
+        i64 =
             ((uint64_t)(t[0][z+1])<< 0) +
             ((uint64_t)(t[1][z+1])<< 8) +
             ((uint64_t)(t[2][z+1])<<16) +
@@ -837,7 +878,8 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[5][z+1])<<40) +
             ((uint64_t)(t[6][z+1])<<48) +
             ((uint64_t)(t[7][z+1])<<56);
-        *(uint64_t *)&out[iN[z+2]] =
+        memcpy(&out[iN[z+1]], &i64, 8);
+        i64 =
             ((uint64_t)(t[0][z+2])<< 0) +
             ((uint64_t)(t[1][z+2])<< 8) +
             ((uint64_t)(t[2][z+2])<<16) +
@@ -846,7 +888,8 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[5][z+2])<<40) +
             ((uint64_t)(t[6][z+2])<<48) +
             ((uint64_t)(t[7][z+2])<<56);
-        *(uint64_t *)&out[iN[z+3]] =
+        memcpy(&out[iN[z+2]], &i64, 8);
+        i64 =
             ((uint64_t)(t[0][z+3])<< 0) +
             ((uint64_t)(t[1][z+3])<< 8) +
             ((uint64_t)(t[2][z+3])<<16) +
@@ -855,8 +898,9 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[5][z+3])<<40) +
             ((uint64_t)(t[6][z+3])<<48) +
             ((uint64_t)(t[7][z+3])<<56);
+        memcpy(&out[iN[z+3]], &i64, 8);
 
-        *(uint64_t *)&out[iN[z]+8] =
+        i64 =
             ((uint64_t)(t[8+0][z])<< 0) +
             ((uint64_t)(t[8+1][z])<< 8) +
             ((uint64_t)(t[8+2][z])<<16) +
@@ -865,7 +909,8 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[8+5][z])<<40) +
             ((uint64_t)(t[8+6][z])<<48) +
             ((uint64_t)(t[8+7][z])<<56);
-        *(uint64_t *)&out[iN[z+1]+8] =
+        memcpy(&out[iN[z]+8], &i64, 8);
+        i64 =
             ((uint64_t)(t[8+0][z+1])<< 0) +
             ((uint64_t)(t[8+1][z+1])<< 8) +
             ((uint64_t)(t[8+2][z+1])<<16) +
@@ -874,7 +919,8 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[8+5][z+1])<<40) +
             ((uint64_t)(t[8+6][z+1])<<48) +
             ((uint64_t)(t[8+7][z+1])<<56);
-        *(uint64_t *)&out[iN[z+2]+8] =
+        memcpy(&out[iN[z+1]+8], &i64, 8);
+        i64 =
             ((uint64_t)(t[8+0][z+2])<< 0) +
             ((uint64_t)(t[8+1][z+2])<< 8) +
             ((uint64_t)(t[8+2][z+2])<<16) +
@@ -883,7 +929,8 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[8+5][z+2])<<40) +
             ((uint64_t)(t[8+6][z+2])<<48) +
             ((uint64_t)(t[8+7][z+2])<<56);
-        *(uint64_t *)&out[iN[z+3]+8] =
+        memcpy(&out[iN[z+2]+8], &i64, 8);
+        i64 =
             ((uint64_t)(t[8+0][z+3])<< 0) +
             ((uint64_t)(t[8+1][z+3])<< 8) +
             ((uint64_t)(t[8+2][z+3])<<16) +
@@ -892,8 +939,9 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[8+5][z+3])<<40) +
             ((uint64_t)(t[8+6][z+3])<<48) +
             ((uint64_t)(t[8+7][z+3])<<56);
+        memcpy(&out[iN[z+3]+8], &i64, 8);
 
-        *(uint64_t *)&out[iN[z]+16] =
+        i64 =
             ((uint64_t)(t[16+0][z])<< 0) +
             ((uint64_t)(t[16+1][z])<< 8) +
             ((uint64_t)(t[16+2][z])<<16) +
@@ -902,7 +950,8 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[16+5][z])<<40) +
             ((uint64_t)(t[16+6][z])<<48) +
             ((uint64_t)(t[16+7][z])<<56);
-        *(uint64_t *)&out[iN[z+1]+16] =
+        memcpy(&out[iN[z]+16], &i64, 8);
+        i64 =
             ((uint64_t)(t[16+0][z+1])<< 0) +
             ((uint64_t)(t[16+1][z+1])<< 8) +
             ((uint64_t)(t[16+2][z+1])<<16) +
@@ -911,7 +960,8 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[16+5][z+1])<<40) +
             ((uint64_t)(t[16+6][z+1])<<48) +
             ((uint64_t)(t[16+7][z+1])<<56);
-        *(uint64_t *)&out[iN[z+2]+16] =
+        memcpy(&out[iN[z+1]+16], &i64, 8);
+        i64 =
             ((uint64_t)(t[16+0][z+2])<< 0) +
             ((uint64_t)(t[16+1][z+2])<< 8) +
             ((uint64_t)(t[16+2][z+2])<<16) +
@@ -920,7 +970,8 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[16+5][z+2])<<40) +
             ((uint64_t)(t[16+6][z+2])<<48) +
             ((uint64_t)(t[16+7][z+2])<<56);
-        *(uint64_t *)&out[iN[z+3]+16] =
+        memcpy(&out[iN[z+2]+16], &i64, 8);
+        i64 =
             ((uint64_t)(t[16+0][z+3])<< 0) +
             ((uint64_t)(t[16+1][z+3])<< 8) +
             ((uint64_t)(t[16+2][z+3])<<16) +
@@ -929,8 +980,9 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[16+5][z+3])<<40) +
             ((uint64_t)(t[16+6][z+3])<<48) +
             ((uint64_t)(t[16+7][z+3])<<56);
+        memcpy(&out[iN[z+3]+16], &i64, 8);
 
-        *(uint64_t *)&out[iN[z]+24] =
+        i64 =
             ((uint64_t)(t[24+0][z])<< 0) +
             ((uint64_t)(t[24+1][z])<< 8) +
             ((uint64_t)(t[24+2][z])<<16) +
@@ -939,7 +991,8 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[24+5][z])<<40) +
             ((uint64_t)(t[24+6][z])<<48) +
             ((uint64_t)(t[24+7][z])<<56);
-        *(uint64_t *)&out[iN[z+1]+24] =
+        memcpy(&out[iN[z]+24], &i64, 8);
+        i64 =
             ((uint64_t)(t[24+0][z+1])<< 0) +
             ((uint64_t)(t[24+1][z+1])<< 8) +
             ((uint64_t)(t[24+2][z+1])<<16) +
@@ -948,7 +1001,8 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[24+5][z+1])<<40) +
             ((uint64_t)(t[24+6][z+1])<<48) +
             ((uint64_t)(t[24+7][z+1])<<56);
-        *(uint64_t *)&out[iN[z+2]+24] =
+        memcpy(&out[iN[z+1]+24], &i64, 8);
+        i64 =
             ((uint64_t)(t[24+0][z+2])<< 0) +
             ((uint64_t)(t[24+1][z+2])<< 8) +
             ((uint64_t)(t[24+2][z+2])<<16) +
@@ -957,7 +1011,8 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[24+5][z+2])<<40) +
             ((uint64_t)(t[24+6][z+2])<<48) +
             ((uint64_t)(t[24+7][z+2])<<56);
-        *(uint64_t *)&out[iN[z+3]+24] =
+        memcpy(&out[iN[z+2]+24], &i64, 8);
+        i64 =
             ((uint64_t)(t[24+0][z+3])<< 0) +
             ((uint64_t)(t[24+1][z+3])<< 8) +
             ((uint64_t)(t[24+2][z+3])<<16) +
@@ -966,13 +1021,13 @@ static inline void transpose_and_copy(uint8_t *out, int iN[32],
             ((uint64_t)(t[24+5][z+3])<<40) +
             ((uint64_t)(t[24+6][z+3])<<48) +
             ((uint64_t)(t[24+7][z+3])<<56);
+        memcpy(&out[iN[z+3]+24], &i64, 8);
 
         iN[z+0] += 32;
         iN[z+1] += 32;
         iN[z+2] += 32;
         iN[z+3] += 32;
     }
-#endif
 }
 
 unsigned char *rans_uncompress_O1_32x16_sse4(unsigned char *in,