Avoid integer undefined behavior (resolves #241)

CHANGELOG.md

@@ -3,6 +3,14 @@ Change Log
 
 ---
 
+## Unreleased
+
+### Fixed
+
+- #241: Signed left shifts, integer overflow invoke undefined behavior.
+
+---
+
 ## 1.0.1 (2023-12-15)
 
 This patch release primarily addresses minor bug fixes and is needed to update

src/template/decode.c

@@ -16,16 +16,27 @@ _t1(inv_lift, Int)(Int* p, ptrdiff_t s)
 
   /*
   ** non-orthogonal transform
+  **
   **       ( 4  6 -4 -1) (x)
   ** 1/4 * ( 4  2  4  5) (y)
   **       ( 4 -2  4 -5) (z)
   **       ( 4 -6 -4  1) (w)
+  **
+  ** original lifted version, which invokes UB due to signed left shift and
+  ** integer overflow:
+  **
+  ** y += w >> 1; w -= y >> 1;
+  ** y += w; w <<= 1; w -= y;
+  ** z += x; x <<= 1; x -= z;
+  ** y += z; z <<= 1; z -= y;
+  ** w += x; x <<= 1; x -= w;
   */
+
   y += w >> 1; w -= y >> 1;
-  y += w; w <<= 1; w -= y;
-  z += x; x <<= 1; x -= z;
-  y += z; z <<= 1; z -= y;
-  w += x; x <<= 1; x -= w;
+  y += w; w -= y - w;
+  z += x; x -= z - x;
+  y += z; z -= y - z;
+  w += x; x -= w - x;
 
   p -= s; *p = w;
   p -= s; *p = z;

src/template/revdecode.c

@@ -6,27 +6,30 @@ static void _t2(rev_inv_xform, Int, DIMS)(Int* p);
 static void
 _t1(rev_inv_lift, Int)(Int* p, ptrdiff_t s)
 {
-  Int x, y, z, w;
+  UInt x, y, z, w;
   x = *p; p += s;
   y = *p; p += s;
   z = *p; p += s;
   w = *p; p += s;
 
   /*
   ** high-order Lorenzo transform (P4 Pascal matrix)
+  **
   ** ( 1  0  0  0) (x)
   ** ( 1  1  0  0) (y)
   ** ( 1  2  1  0) (z)
   ** ( 1  3  3  1) (w)
+  **
+  ** unsigned arithmetic is used to avoid integer overflow
   */
   w += z;
   z += y; w += z;
   y += x; z += y; w += z;
 
-  p -= s; *p = w;
-  p -= s; *p = z;
-  p -= s; *p = y;
-  p -= s; *p = x;
+  p -= s; *p = (Int)w;
+  p -= s; *p = (Int)z;
+  p -= s; *p = (Int)y;
+  p -= s; *p = (Int)x;
 }
 
 /* decode block of integers using reversible algorithm */

src/template/revencode.c

@@ -6,27 +6,30 @@ static void _t2(rev_fwd_xform, Int, DIMS)(Int* p);
 static void
 _t1(rev_fwd_lift, Int)(Int* p, ptrdiff_t s)
 {
-  Int x, y, z, w;
+  UInt x, y, z, w;
   x = *p; p += s;
   y = *p; p += s;
   z = *p; p += s;
   w = *p; p += s;
 
   /*
   ** high-order Lorenzo transform
+  **
   ** ( 1  0  0  0) (x)
   ** (-1  1  0  0) (y)
   ** ( 1 -2  1  0) (z)
   ** (-1  3 -3  1) (w)
+  **
+  ** unsigned arithmetic is used to avoid integer overflow
   */
   w -= z; z -= y; y -= x;
   w -= z; z -= y;
   w -= z;
 
-  p -= s; *p = w;
-  p -= s; *p = z;
-  p -= s; *p = y;
-  p -= s; *p = x;
+  p -= s; *p = (Int)w;
+  p -= s; *p = (Int)z;
+  p -= s; *p = (Int)y;
+  p -= s; *p = (Int)x;
 }
 
 /* return precision required to encode block reversibly */