diff --git a/docs/language-reference.rst b/docs/language-reference.rst
--- a/docs/language-reference.rst
+++ b/docs/language-reference.rst
@@ -291,7 +291,7 @@
 
   def pair 't (x: t) (y: t) = (x, y)
 
-The application ``pair [1] [2,3]`` will fail at run-time.
+The application ``pair [1] [2,3]`` is ill-typed.
 
 To simplify the handling of in-place updates (see
 :ref:`in-place-updates`), the value returned by a function may not
diff --git a/futhark.cabal b/futhark.cabal
--- a/futhark.cabal
+++ b/futhark.cabal
@@ -1,6 +1,6 @@
 cabal-version: 2.4
 name:           futhark
-version:        0.22.5
+version:        0.22.6
 synopsis:       An optimising compiler for a functional, array-oriented language.
 
 description:    Futhark is a small programming language designed to be compiled to
diff --git a/prelude/math.fut b/prelude/math.fut
--- a/prelude/math.fut
+++ b/prelude/math.fut
@@ -332,8 +332,8 @@
   def set_bit (bit: i32) (x: t) (b: i32) =
     ((x & i32 (!(1 intrinsics.<< bit))) | i32 (b intrinsics.<< bit))
   def popc = intrinsics.popc8
-  def mul_hi a b = intrinsics.mul_hi8 (i8 a, i8 b)
-  def mad_hi a b c = intrinsics.mad_hi8 (i8 a, i8 b, i8 c)
+  def mul_hi a b = intrinsics.smul_hi8 (i8 a, i8 b)
+  def mad_hi a b c = intrinsics.smad_hi8 (i8 a, i8 b, i8 c)
   def clz = intrinsics.clz8
   def ctz = intrinsics.ctz8
 
@@ -405,8 +405,8 @@
   def set_bit (bit: i32) (x: t) (b: i32) =
     ((x & i32 (!(1 intrinsics.<< bit))) | i32 (b intrinsics.<< bit))
   def popc = intrinsics.popc16
-  def mul_hi a b = intrinsics.mul_hi16 (i16 a, i16 b)
-  def mad_hi a b c = intrinsics.mad_hi16 (i16 a, i16 b, i16 c)
+  def mul_hi a b = intrinsics.smul_hi16 (i16 a, i16 b)
+  def mad_hi a b c = intrinsics.smad_hi16 (i16 a, i16 b, i16 c)
   def clz = intrinsics.clz16
   def ctz = intrinsics.ctz16
 
@@ -481,8 +481,8 @@
   def set_bit (bit: i32) (x: t) (b: i32) =
     ((x & i32 (!(1 intrinsics.<< bit))) | i32 (b intrinsics.<< bit))
   def popc = intrinsics.popc32
-  def mul_hi a b = intrinsics.mul_hi32 (i32 a, i32 b)
-  def mad_hi a b c = intrinsics.mad_hi32 (i32 a, i32 b, i32 c)
+  def mul_hi a b = intrinsics.smul_hi32 (i32 a, i32 b)
+  def mad_hi a b c = intrinsics.smad_hi32 (i32 a, i32 b, i32 c)
   def clz = intrinsics.clz32
   def ctz = intrinsics.ctz32
 
@@ -557,8 +557,8 @@
   def set_bit (bit: i32) (x: t) (b: i32) =
     ((x & i32 (!(1 intrinsics.<< bit))) | intrinsics.zext_i32_i64 (b intrinsics.<< bit))
   def popc = intrinsics.popc64
-  def mul_hi a b = intrinsics.mul_hi64 (i64 a, i64 b)
-  def mad_hi a b c = intrinsics.mad_hi64 (i64 a, i64 b, i64 c)
+  def mul_hi a b = intrinsics.smul_hi64 (i64 a, i64 b)
+  def mad_hi a b c = intrinsics.smad_hi64 (i64 a, i64 b, i64 c)
   def clz = intrinsics.clz64
   def ctz = intrinsics.ctz64
 
@@ -633,8 +633,8 @@
   def set_bit (bit: i32) (x: t) (b: i32) =
     ((x & i32 (!(1 intrinsics.<< bit))) | i32 (b intrinsics.<< bit))
   def popc x = intrinsics.popc8 (sign x)
-  def mul_hi a b = unsign (intrinsics.mul_hi8 (sign a, sign b))
-  def mad_hi a b c = unsign (intrinsics.mad_hi8 (sign a, sign b, sign c))
+  def mul_hi a b = unsign (intrinsics.umul_hi8 (sign a, sign b))
+  def mad_hi a b c = unsign (intrinsics.umad_hi8 (sign a, sign b, sign c))
   def clz x = intrinsics.clz8 (sign x)
   def ctz x = intrinsics.ctz8 (sign x)
 
@@ -709,8 +709,8 @@
   def set_bit (bit: i32) (x: t) (b: i32) =
     ((x & i32 (!(1 intrinsics.<< bit))) | i32 (b intrinsics.<< bit))
   def popc x = intrinsics.popc16 (sign x)
-  def mul_hi a b = unsign (intrinsics.mul_hi16 (sign a, sign b))
-  def mad_hi a b c = unsign (intrinsics.mad_hi16 (sign a, sign b, sign c))
+  def mul_hi a b = unsign (intrinsics.umul_hi16 (sign a, sign b))
+  def mad_hi a b c = unsign (intrinsics.umad_hi16 (sign a, sign b, sign c))
   def clz x = intrinsics.clz16 (sign x)
   def ctz x = intrinsics.ctz16 (sign x)
 
@@ -785,8 +785,8 @@
   def set_bit (bit: i32) (x: t) (b: i32) =
     ((x & i32 (!(1 intrinsics.<< bit))) | i32 (b intrinsics.<< bit))
   def popc x = intrinsics.popc32 (sign x)
-  def mul_hi a b = unsign (intrinsics.mul_hi32 (sign a, sign b))
-  def mad_hi a b c = unsign (intrinsics.mad_hi32 (sign a, sign b, sign c))
+  def mul_hi a b = unsign (intrinsics.umul_hi32 (sign a, sign b))
+  def mad_hi a b c = unsign (intrinsics.umad_hi32 (sign a, sign b, sign c))
   def clz x = intrinsics.clz32 (sign x)
   def ctz x = intrinsics.ctz32 (sign x)
 
@@ -861,8 +861,8 @@
   def set_bit (bit: i32) (x: t) (b: i32) =
     ((x & i32 (!(1 intrinsics.<< bit))) | i32 (b intrinsics.<< bit))
   def popc x = intrinsics.popc64 (sign x)
-  def mul_hi a b = unsign (intrinsics.mul_hi64 (sign a, sign b))
-  def mad_hi a b c = unsign (intrinsics.mad_hi64 (sign a, sign b, sign c))
+  def mul_hi a b = unsign (intrinsics.umul_hi64 (sign a, sign b))
+  def mad_hi a b c = unsign (intrinsics.umad_hi64 (sign a, sign b, sign c))
   def clz x = intrinsics.clz64 (sign x)
   def ctz x = intrinsics.ctz64 (sign x)
 
diff --git a/rts/c/scalar.h b/rts/c/scalar.h
--- a/rts/c/scalar.h
+++ b/rts/c/scalar.h
@@ -1289,70 +1289,28 @@
 #endif
 
 #if defined(__OPENCL_VERSION__)
-static uint8_t futrts_mul_hi8(uint8_t a, uint8_t b) {
-  return mul_hi(a, b);
-}
-
-static uint16_t futrts_mul_hi16(uint16_t a, uint16_t b) {
-  return mul_hi(a, b);
-}
-
-static uint32_t futrts_mul_hi32(uint32_t a, uint32_t b) {
-  return mul_hi(a, b);
-}
-
-static uint64_t futrts_mul_hi64(uint64_t a, uint64_t b) {
-  return mul_hi(a, b);
-}
-
+static uint8_t  futrts_umul_hi8 ( uint8_t a,  uint8_t b) { return mul_hi(a, b); }
+static uint16_t futrts_umul_hi16(uint16_t a, uint16_t b) { return mul_hi(a, b); }
+static uint32_t futrts_umul_hi32(uint32_t a, uint32_t b) { return mul_hi(a, b); }
+static uint64_t futrts_umul_hi64(uint64_t a, uint64_t b) { return mul_hi(a, b); }
+static uint8_t  futrts_smul_hi8 ( int8_t a,  int8_t b) { return mul_hi(a, b); }
+static uint16_t futrts_smul_hi16(int16_t a, int16_t b) { return mul_hi(a, b); }
+static uint32_t futrts_smul_hi32(int32_t a, int32_t b) { return mul_hi(a, b); }
+static uint64_t futrts_smul_hi64(int64_t a, int64_t b) { return mul_hi(a, b); }
 #elif defined(__CUDA_ARCH__)
-
-static uint8_t futrts_mul_hi8(uint8_t a, uint8_t b) {
-  uint16_t aa = a;
-  uint16_t bb = b;
-
-  return aa * bb >> 8;
-}
-
-static uint16_t futrts_mul_hi16(uint16_t a, uint16_t b) {
-  uint32_t aa = a;
-  uint32_t bb = b;
-
-  return aa * bb >> 16;
-}
-
-static uint32_t futrts_mul_hi32(uint32_t a, uint32_t b) {
-  return mulhi(a, b);
-}
-
-static uint64_t futrts_mul_hi64(uint64_t a, uint64_t b) {
-  return mul64hi(a, b);
-}
-
+static  uint8_t futrts_umul_hi8(uint8_t a, uint8_t b) { return ((uint16_t)a) * ((uint16_t)b) >> 8; }
+static uint16_t futrts_umul_hi16(uint16_t a, uint16_t b) { return ((uint32_t)a) * ((uint32_t)b) >> 16; }
+static uint32_t futrts_umul_hi32(uint32_t a, uint32_t b) { return __umulhi(a, b); }
+static uint64_t futrts_umul_hi64(uint64_t a, uint64_t b) { return __umul64hi(a, b); }
+static  uint8_t futrts_smul_hi8 ( int8_t a, int8_t b) { return ((int16_t)a) * ((int16_t)b) >> 8; }
+static uint16_t futrts_smul_hi16(int16_t a, int16_t b) { return ((int32_t)a) * ((int32_t)b) >> 16; }
+static uint32_t futrts_smul_hi32(int32_t a, int32_t b) { return __mulhi(a, b); }
+static uint64_t futrts_smul_hi64(int64_t a, int64_t b) { return __mul64hi(a, b); }
 #elif ISPC
-
-static uint8_t futrts_mul_hi8(uint8_t a, uint8_t b) {
-  uint16_t aa = a;
-  uint16_t bb = b;
-
-  return aa * bb >> 8;
-}
-
-static uint16_t futrts_mul_hi16(uint16_t a, uint16_t b) {
-  uint32_t aa = a;
-  uint32_t bb = b;
-
-  return aa * bb >> 16;
-}
-
-static uint32_t futrts_mul_hi32(uint32_t a, uint32_t b) {
-  uint64_t aa = a;
-  uint64_t bb = b;
-
-  return aa * bb >> 32;
-}
-
-static uint64_t futrts_mul_hi64(uint64_t a, uint64_t b) {
+static uint8_t futrts_umul_hi8(uint8_t a, uint8_t b) { return ((uint16_t)a) * ((uint16_t)b) >> 8; }
+static uint16_t futrts_umul_hi16(uint16_t a, uint16_t b) { return ((uint32_t)a) * ((uint32_t)b) >> 16; }
+static uint32_t futrts_umul_hi32(uint32_t a, uint32_t b) { return ((uint64_t)a) * ((uint64_t)b) >> 32; }
+static uint64_t futrts_umul_hi64(uint64_t a, uint64_t b) {
   uint64_t ah = a >> 32;
   uint64_t al = a & 0xffffffff;
   uint64_t bh = b >> 32;
@@ -1369,78 +1327,69 @@
   uint64_t p2l = p2 & 0xffffffff;
   uint64_t p3l = p3 & 0xffffffff;
 
-  uint64_t l = p1h + p2l  + p3l;
+  uint64_t l = p1h + p2l + p3l;
   uint64_t m = (p2 >> 32) + (p3 >> 32);
   uint64_t h = (l >> 32) + m + p4;
 
   return h;
 }
-
-#else // Not OpenCL, ISPC, or CUDA, but plain C.
-
-static uint8_t futrts_mul_hi8(uint8_t a, uint8_t b) {
-  uint16_t aa = a;
-  uint16_t bb = b;
-
-  return aa * bb >> 8;
-}
+static  int8_t futrts_smul_hi8 ( int8_t a,  int8_t b) { return ((uint16_t)a) * ((uint16_t)b) >> 8; }
+static int16_t futrts_smul_hi16(int16_t a, int16_t b) { return ((uint32_t)a) * ((uint32_t)b) >> 16; }
+static int32_t futrts_smul_hi32(int32_t a, int32_t b) { return ((uint64_t)a) * ((uint64_t)b) >> 32; }
+static int64_t futrts_smul_hi64(int64_t a, int64_t b) {
+  uint64_t ah = a >> 32;
+  uint64_t al = a & 0xffffffff;
+  uint64_t bh = b >> 32;
+  uint64_t bl = b & 0xffffffff;
 
-static uint16_t futrts_mul_hi16(uint16_t a, uint16_t b) {
-  uint32_t aa = a;
-  uint32_t bb = b;
+  uint64_t p1 =  al * bl;
+  int64_t  p2 = al * bh;
+  int64_t  p3 = ah * bl;
+  uint64_t p4 =  ah * bh;
 
-  return aa * bb >> 16;
-}
+  uint64_t p1h = p1 >> 32;
+  uint64_t p2h = p2 >> 32;
+  uint64_t p3h = p3 >> 32;
+  uint64_t p2l = p2 & 0xffffffff;
+  uint64_t p3l = p3 & 0xffffffff;
 
-static uint32_t futrts_mul_hi32(uint32_t a, uint32_t b) {
-  uint64_t aa = a;
-  uint64_t bb = b;
+  uint64_t l = p1h + p2l + p3l;
+  uint64_t m = (p2 >> 32) + (p3 >> 32);
+  uint64_t h = (l >> 32) + m + p4;
 
-  return aa * bb >> 32;
+  return h;
 }
 
-static uint64_t futrts_mul_hi64(uint64_t a, uint64_t b) {
-  __uint128_t aa = a;
-  __uint128_t bb = b;
-
-  return aa * bb >> 64;
-}
+#else // Not OpenCL, ISPC, or CUDA, but plain C.
+static uint8_t futrts_umul_hi8(uint8_t a, uint8_t b) { return ((uint16_t)a) * ((uint16_t)b) >> 8; }
+static uint16_t futrts_umul_hi16(uint16_t a, uint16_t b) { return ((uint32_t)a) * ((uint32_t)b) >> 16; }
+static uint32_t futrts_umul_hi32(uint32_t a, uint32_t b) { return ((uint64_t)a) * ((uint64_t)b) >> 32; }
+static uint64_t futrts_umul_hi64(uint64_t a, uint64_t b) { return ((__uint128_t)a) * ((__uint128_t)b) >> 64; }
+static int8_t futrts_smul_hi8(int8_t a, int8_t b) { return ((int16_t)a) * ((int16_t)b) >> 8; }
+static int16_t futrts_smul_hi16(int16_t a, int16_t b) { return ((int32_t)a) * ((int32_t)b) >> 16; }
+static int32_t futrts_smul_hi32(int32_t a, int32_t b) { return ((int64_t)a) * ((int64_t)b) >> 32; }
+static int64_t futrts_smul_hi64(int64_t a, int64_t b) { return ((__int128_t)a) * ((__int128_t)b) >> 64; }
 #endif
 
 #if defined(__OPENCL_VERSION__)
-static uint8_t futrts_mad_hi8(uint8_t a, uint8_t b, uint8_t c) {
-  return mad_hi(a, b, c);
-}
-
-static uint16_t futrts_mad_hi16(uint16_t a, uint16_t b, uint16_t c) {
-  return mad_hi(a, b, c);
-}
-
-static uint32_t futrts_mad_hi32(uint32_t a, uint32_t b, uint32_t c) {
-  return mad_hi(a, b, c);
-}
-
-static uint64_t futrts_mad_hi64(uint64_t a, uint64_t b, uint64_t c) {
-  return mad_hi(a, b, c);
-}
-
+static  uint8_t futrts_umad_hi8 ( uint8_t a,  uint8_t b,  uint8_t c) { return mad_hi(a, b, c); }
+static uint16_t futrts_umad_hi16(uint16_t a, uint16_t b, uint16_t c) { return mad_hi(a, b, c); }
+static uint32_t futrts_umad_hi32(uint32_t a, uint32_t b, uint32_t c) { return mad_hi(a, b, c); }
+static uint64_t futrts_umad_hi64(uint64_t a, uint64_t b, uint64_t c) { return mad_hi(a, b, c); }
+static  uint8_t futrts_smad_hi8( int8_t a,  int8_t b,   int8_t c) { return mad_hi(a, b, c); }
+static uint16_t futrts_smad_hi16(int16_t a, int16_t b, int16_t c) { return mad_hi(a, b, c); }
+static uint32_t futrts_smad_hi32(int32_t a, int32_t b, int32_t c) { return mad_hi(a, b, c); }
+static uint64_t futrts_smad_hi64(int64_t a, int64_t b, int64_t c) { return mad_hi(a, b, c); }
 #else // Not OpenCL
 
-static uint8_t futrts_mad_hi8(uint8_t a, uint8_t b, uint8_t c) {
-  return futrts_mul_hi8(a, b) + c;
-}
-
-static uint16_t futrts_mad_hi16(uint16_t a, uint16_t b, uint16_t c) {
-  return futrts_mul_hi16(a, b) + c;
-}
-
-static uint32_t futrts_mad_hi32(uint32_t a, uint32_t b, uint32_t c) {
-  return futrts_mul_hi32(a, b) + c;
-}
-
-static uint64_t futrts_mad_hi64(uint64_t a, uint64_t b, uint64_t c) {
-  return futrts_mul_hi64(a, b) + c;
-}
+static  uint8_t futrts_umad_hi8( uint8_t a,  uint8_t b,  uint8_t c) { return futrts_umul_hi8(a, b) + c; }
+static uint16_t futrts_umad_hi16(uint16_t a, uint16_t b, uint16_t c) { return futrts_umul_hi16(a, b) + c; }
+static uint32_t futrts_umad_hi32(uint32_t a, uint32_t b, uint32_t c) { return futrts_umul_hi32(a, b) + c; }
+static uint64_t futrts_umad_hi64(uint64_t a, uint64_t b, uint64_t c) { return futrts_umul_hi64(a, b) + c; }
+static  uint8_t futrts_smad_hi8 ( int8_t a,  int8_t b,  int8_t c) { return futrts_smul_hi8(a, b) + c; }
+static uint16_t futrts_smad_hi16(int16_t a, int16_t b, int16_t c) { return futrts_smul_hi16(a, b) + c; }
+static uint32_t futrts_smad_hi32(int32_t a, int32_t b, int32_t c) { return futrts_smul_hi32(a, b) + c; }
+static uint64_t futrts_smad_hi64(int64_t a, int64_t b, int64_t c) { return futrts_smul_hi64(a, b) + c; }
 #endif
 
 #if defined(__OPENCL_VERSION__)
diff --git a/rts/python/scalar.py b/rts/python/scalar.py
--- a/rts/python/scalar.py
+++ b/rts/python/scalar.py
@@ -394,38 +394,38 @@
 def fpconv_f64_f32(x):
   return np.float32(x)
 
-def futhark_mul_hi8(a, b):
-  a = np.uint64(np.uint8(a))
-  b = np.uint64(np.uint8(b))
-  return np.int8((a*b) >> np.uint64(8))
+def futhark_umul_hi8(a, b):
+  return np.int8((np.uint64(np.uint8(a))*np.uint64(np.uint8(b))) >> np.uint64(8))
 
-def futhark_mul_hi16(a, b):
-  a = np.uint64(np.uint16(a))
-  b = np.uint64(np.uint16(b))
-  return np.int16((a*b) >> np.uint64(16))
+def futhark_umul_hi16(a, b):
+  return np.int16((np.uint64(np.uint16(a))*np.uint64(np.uint16(b))) >> np.uint64(16))
 
-def futhark_mul_hi32(a, b):
-  a = np.uint64(np.uint32(a))
-  b = np.uint64(np.uint32(b))
-  return np.int32((a*b) >> np.uint64(32))
+def futhark_umul_hi32(a, b):
+  return np.int32((np.uint64(np.uint32(a))*np.uint64(np.uint32(b))) >> np.uint64(32))
 
-# This one is done with arbitrary-precision integers.
-def futhark_mul_hi64(a, b):
-  a = int(np.uint64(a))
-  b = int(np.uint64(b))
-  return np.int64(np.uint64(a*b >> 64))
+def futhark_umul_hi64(a, b):
+  return np.int64(np.uint64(int(np.uint64(a))*int(np.uint64(b)) >> 64))
 
-def futhark_mad_hi8(a, b, c):
-  return futhark_mul_hi8(a,b) + c
+def futhark_smul_hi8(a, b):
+  return np.int8((np.int64(a)*np.int64(b)) >> np.int64(8))
 
-def futhark_mad_hi16(a, b, c):
-  return futhark_mul_hi16(a,b) + c
+def futhark_smul_hi16(a, b):
+  return np.int16((np.int64(a)*np.int64(b)) >> np.int64(16))
 
-def futhark_mad_hi32(a, b, c):
-  return futhark_mul_hi32(a,b) + c
+def futhark_smul_hi32(a, b):
+  return np.int32((np.int64(a)*np.int64(b)) >> np.int64(32))
 
-def futhark_mad_hi64(a, b, c):
-  return futhark_mul_hi64(a,b) + c
+def futhark_smul_hi64(a, b):
+  return np.int64(int(a)*int(b) >> 64)
+
+def futhark_umad_hi8(a, b, c):  return futhark_umul_hi8(a,b) + c
+def futhark_umad_hi16(a, b, c): return futhark_umul_hi16(a,b) + c
+def futhark_umad_hi32(a, b, c): return futhark_umul_hi32(a,b) + c
+def futhark_umad_hi64(a, b, c): return futhark_umul_hi64(a,b) + c
+def futhark_smad_hi8(a, b, c):  return futhark_smul_hi8(a,b) + c
+def futhark_smad_hi16(a, b, c): return futhark_smul_hi16(a,b) + c
+def futhark_smad_hi32(a, b, c): return futhark_smul_hi32(a,b) + c
+def futhark_smad_hi64(a, b, c): return futhark_smul_hi64(a,b) + c
 
 def futhark_log64(x):
   return np.float64(np.log(x))
diff --git a/src/Futhark/AD/Derivatives.hs b/src/Futhark/AD/Derivatives.hs
--- a/src/Futhark/AD/Derivatives.hs
+++ b/src/Futhark/AD/Derivatives.hs
@@ -346,14 +346,22 @@
 pdBuiltin "erfc64" [z] =
   Just [untyped $ negate $ (2 / sqrt pi) * exp (negate (isF64 z * isF64 z))]
 -- More problematic derivatives follow below.
-pdBuiltin "mul_hi8" [x, y] = Just [y, x]
-pdBuiltin "mul_hi16" [x, y] = Just [y, x]
-pdBuiltin "mul_hi32" [x, y] = Just [y, x]
-pdBuiltin "mul_hi64" [x, y] = Just [y, x]
-pdBuiltin "mad_hi8" [a, b, _c] = Just [b, a, iConst Int8 1]
-pdBuiltin "mad_hi16" [a, b, _c] = Just [b, a, iConst Int16 1]
-pdBuiltin "mad_hi32" [a, b, _c] = Just [b, a, iConst Int32 1]
-pdBuiltin "mad_hi64" [a, b, _c] = Just [b, a, iConst Int64 1]
+pdBuiltin "umul_hi8" [x, y] = Just [y, x]
+pdBuiltin "umul_hi16" [x, y] = Just [y, x]
+pdBuiltin "umul_hi32" [x, y] = Just [y, x]
+pdBuiltin "umul_hi64" [x, y] = Just [y, x]
+pdBuiltin "umad_hi8" [a, b, _c] = Just [b, a, iConst Int8 1]
+pdBuiltin "umad_hi16" [a, b, _c] = Just [b, a, iConst Int16 1]
+pdBuiltin "umad_hi32" [a, b, _c] = Just [b, a, iConst Int32 1]
+pdBuiltin "umad_hi64" [a, b, _c] = Just [b, a, iConst Int64 1]
+pdBuiltin "smul_hi8" [x, y] = Just [y, x]
+pdBuiltin "smul_hi16" [x, y] = Just [y, x]
+pdBuiltin "smul_hi32" [x, y] = Just [y, x]
+pdBuiltin "smul_hi64" [x, y] = Just [y, x]
+pdBuiltin "smad_hi8" [a, b, _c] = Just [b, a, iConst Int8 1]
+pdBuiltin "smad_hi16" [a, b, _c] = Just [b, a, iConst Int16 1]
+pdBuiltin "smad_hi32" [a, b, _c] = Just [b, a, iConst Int32 1]
+pdBuiltin "smad_hi64" [a, b, _c] = Just [b, a, iConst Int64 1]
 pdBuiltin "isnan16" [_] = Just [untyped false]
 pdBuiltin "isnan32" [_] = Just [untyped false]
 pdBuiltin "isnan64" [_] = Just [untyped false]
diff --git a/src/Futhark/CLI/REPL.hs b/src/Futhark/CLI/REPL.hs
--- a/src/Futhark/CLI/REPL.hs
+++ b/src/Futhark/CLI/REPL.hs
@@ -20,7 +20,7 @@
 import Futhark.MonadFreshNames
 import Futhark.Util (fancyTerminal)
 import Futhark.Util.Options
-import Futhark.Util.Pretty (AnsiStyle, Color (..), Doc, annotate, bgColorDull, bold, brackets, color, docText, docTextForHandle, hardline, pretty, putDoc, putDocLn, (<+>))
+import Futhark.Util.Pretty (AnsiStyle, Color (..), Doc, align, annotate, bgColorDull, bold, brackets, color, docText, docTextForHandle, hardline, pretty, putDoc, putDocLn, unAnnotate, (<+>))
 import Futhark.Version
 import Language.Futhark
 import Language.Futhark.Interpreter qualified as I
@@ -324,7 +324,7 @@
     intOp (I.ExtOpError err) =
       pure $ Left err
     intOp (I.ExtOpTrace w v c) = do
-      liftIO $ putStrLn $ w ++ ": " ++ v
+      liftIO $ putDocLn $ pretty w <> ":" <+> unAnnotate v
       c
     intOp (I.ExtOpBreak w why callstack c) = do
       s <- get
@@ -376,7 +376,7 @@
   where
     intOp (I.ExtOpError err) = pure $ Left err
     intOp (I.ExtOpTrace w v c) = do
-      liftIO $ putStrLn $ w ++ ": " ++ v
+      liftIO $ putDocLn $ pretty w <> ":" <+> align (unAnnotate v)
       c
     intOp (I.ExtOpBreak _ _ _ c) = c
 
diff --git a/src/Futhark/CLI/Run.hs b/src/Futhark/CLI/Run.hs
--- a/src/Futhark/CLI/Run.hs
+++ b/src/Futhark/CLI/Run.hs
@@ -14,7 +14,7 @@
 import Futhark.Pipeline
 import Futhark.Util (toPOSIX)
 import Futhark.Util.Options
-import Futhark.Util.Pretty (AnsiStyle, Doc, hPutDoc)
+import Futhark.Util.Pretty (AnsiStyle, Doc, align, hPutDoc, hPutDocLn, pretty, unAnnotate, (<+>))
 import Language.Futhark
 import Language.Futhark.Interpreter qualified as I
 import Language.Futhark.Semantic qualified as T
@@ -154,6 +154,6 @@
   where
     intOp (I.ExtOpError err) = pure $ Left err
     intOp (I.ExtOpTrace w v c) = do
-      liftIO $ hPutStrLn stderr $ w ++ ": " ++ v
+      liftIO $ hPutDocLn stderr $ pretty w <> ":" <+> align (unAnnotate v)
       c
     intOp (I.ExtOpBreak _ _ _ c) = c
diff --git a/src/Futhark/CodeGen/Backends/MulticoreISPC.hs b/src/Futhark/CodeGen/Backends/MulticoreISPC.hs
--- a/src/Futhark/CodeGen/Backends/MulticoreISPC.hs
+++ b/src/Futhark/CodeGen/Backends/MulticoreISPC.hs
@@ -551,14 +551,9 @@
 compileCode (Write dest (Count idx) elemtype DefaultSpace _ elemexp)
   | isConstExp (untyped idx) = do
       dest' <- GC.rawMem dest
-      idxexp <- compileExp (untyped idx)
-      varis <- mapM getVariability (namesToList $ freeIn idx)
-      let quals = if all (== Uniform) varis then [C.ctyquals|$tyqual:uniform|] else []
-      tmp <- newVName "tmp_idx"
-      -- Disambiguate the variability of the constant index
-      GC.decl [C.cdecl|$tyquals:quals typename int64_t $id:tmp = $exp:idxexp;|]
+      idxexp <- compileExp $ constFoldPrimExp $ untyped idx
       deref <-
-        GC.derefPointer dest' [C.cexp|$id:tmp|]
+        GC.derefPointer dest' [C.cexp|($tyquals:([varying]) typename int64_t)$exp:idxexp|]
           <$> getMemType dest elemtype
       elemexp' <- toStorage elemtype <$> compileExp elemexp
       GC.stm [C.cstm|$exp:deref = $exp:elemexp';|]
diff --git a/src/Futhark/CodeGen/ImpGen.hs b/src/Futhark/CodeGen/ImpGen.hs
--- a/src/Futhark/CodeGen/ImpGen.hs
+++ b/src/Futhark/CodeGen/ImpGen.hs
@@ -806,9 +806,9 @@
   Exp rep ->
   ImpM rep r op ()
 defCompileExp pat (Match ses cases defbody _) =
-  foldl f (compileBody pat defbody) cases
+  foldr f (compileBody pat defbody) cases
   where
-    f rest (Case vs body) = sIf (caseMatch ses vs) (compileBody pat body) rest
+    f (Case vs body) = sIf (caseMatch ses vs) (compileBody pat body)
 defCompileExp pat (Apply fname args _ _) = do
   dest <- destinationFromPat pat
   targets <- funcallTargets dest
diff --git a/src/Futhark/CodeGen/ImpGen/Multicore/Base.hs b/src/Futhark/CodeGen/ImpGen/Multicore/Base.hs
--- a/src/Futhark/CodeGen/ImpGen/Multicore/Base.hs
+++ b/src/Futhark/CodeGen/ImpGen/Multicore/Base.hs
@@ -7,7 +7,6 @@
     MulticoreGen,
     decideScheduling,
     decideScheduling',
-    groupResultArrays,
     renameSegBinOp,
     freeParams,
     renameHistOpLambda,
@@ -149,18 +148,6 @@
   ts <- mapM lookupType free
   concat <$> zipWithM toParam free ts
 
--- | Arrays for storing group results shared between threads
-groupResultArrays ::
-  String ->
-  SubExp ->
-  [SegBinOp MCMem] ->
-  MulticoreGen [[VName]]
-groupResultArrays s num_threads reds =
-  forM reds $ \(SegBinOp _ lam _ shape) ->
-    forM (lambdaReturnType lam) $ \t -> do
-      let full_shape = Shape [num_threads] <> shape <> arrayShape t
-      sAllocArray s (elemType t) full_shape DefaultSpace
-
 isLoadBalanced :: Imp.MCCode -> Bool
 isLoadBalanced (a Imp.:>>: b) = isLoadBalanced a && isLoadBalanced b
 isLoadBalanced (Imp.For _ _ a) = isLoadBalanced a
@@ -173,16 +160,10 @@
 isLoadBalanced (Imp.Op (Imp.ISPCKernel a _)) = isLoadBalanced a
 isLoadBalanced _ = True
 
-segBinOpComm' :: [SegBinOp rep] -> Commutativity
-segBinOpComm' = mconcat . map segBinOpComm
-
 decideScheduling' :: SegOp () rep -> Imp.MCCode -> Imp.Scheduling
 decideScheduling' SegHist {} _ = Imp.Static
 decideScheduling' SegScan {} _ = Imp.Static
-decideScheduling' (SegRed _ _ reds _ _) code =
-  case segBinOpComm' reds of
-    Commutative -> decideScheduling code
-    Noncommutative -> Imp.Static
+decideScheduling' SegRed {} _ = Imp.Static
 decideScheduling' SegMap {} code = decideScheduling code
 
 decideScheduling :: Imp.MCCode -> Imp.Scheduling
diff --git a/src/Futhark/CodeGen/ImpGen/Multicore/SegRed.hs b/src/Futhark/CodeGen/ImpGen/Multicore/SegRed.hs
--- a/src/Futhark/CodeGen/ImpGen/Multicore/SegRed.hs
+++ b/src/Futhark/CodeGen/ImpGen/Multicore/SegRed.hs
@@ -80,6 +80,18 @@
   let op' = op {segBinOpLambda = lambda'}
   pure slug {slugOp = op'}
 
+-- | Arrays for storing group results shared between threads
+groupResultArrays ::
+  String ->
+  SubExp ->
+  [SegBinOp MCMem] ->
+  MulticoreGen [[VName]]
+groupResultArrays s num_threads reds =
+  forM reds $ \(SegBinOp _ lam _ shape) ->
+    forM (lambdaReturnType lam) $ \t -> do
+      let full_shape = Shape [num_threads] <> shape <> arrayShape t
+      sAllocArray s (elemType t) full_shape DefaultSpace
+
 nonsegmentedReduction ::
   Pat LetDecMem ->
   SegSpace ->
@@ -216,7 +228,6 @@
     forM_ (zip3 all_red_res' slugs slug_local_accs) $ \(red_res, slug, local_accs) ->
       getNestLoop typ (slugShape slug) $ \vec_is -> do
         let lamtypes = lambdaReturnType $ segBinOpLambda $ slugOp slug
-        -- Load accum params
         getRedLoop typ i $ \uni -> do
           sComment "Load accum params" $
             forM_ (zip3 (accParams slug) local_accs lamtypes) $
diff --git a/src/Futhark/IR/Prop/Types.hs b/src/Futhark/IR/Prop/Types.hs
--- a/src/Futhark/IR/Prop/Types.hs
+++ b/src/Futhark/IR/Prop/Types.hs
@@ -4,6 +4,7 @@
     arrayRank,
     arrayShape,
     setArrayShape,
+    isEmptyArray,
     existential,
     uniqueness,
     unique,
@@ -116,6 +117,13 @@
   newshape ->
   TypeBase newshape u
 setArrayShape t ds = modifyArrayShape (const ds) t
+
+-- | If the array is statically an empty array (meaning any dimension
+-- is a static zero), return the element type and the shape.
+isEmptyArray :: Type -> Maybe (PrimType, Shape)
+isEmptyArray (Array pt (Shape ds) _)
+  | intConst Int64 0 `elem` ds = Just (pt, Shape ds)
+isEmptyArray _ = Nothing
 
 -- | True if the given type has a dimension that is existentially sized.
 existential :: ExtType -> Bool
diff --git a/src/Futhark/IR/SOACS/Simplify.hs b/src/Futhark/IR/SOACS/Simplify.hs
--- a/src/Futhark/IR/SOACS/Simplify.hs
+++ b/src/Futhark/IR/SOACS/Simplify.hs
@@ -400,14 +400,21 @@
 removeUnusedSOACInput _ pat aux op
   | Just (Screma w arrs form :: SOAC rep) <- asSOAC op,
     ScremaForm scan reduce map_lam <- form,
-    (used, unused) <- partition (usedInput map_lam) (zip (lambdaParams map_lam) arrs),
-    not (null unused) = Simplify $ do
-      let (used_params, used_arrs) = unzip used
-          map_lam' = map_lam {lambdaParams = used_params}
-      auxing aux $ letBind pat $ Op $ soacOp $ Screma w used_arrs (ScremaForm scan reduce map_lam')
+    Just (used_arrs, map_lam') <- remove map_lam arrs =
+      Simplify . auxing aux . letBind pat . Op $
+        soacOp (Screma w used_arrs (ScremaForm scan reduce map_lam'))
+  | Just (Scatter w arrs map_lam dests :: SOAC rep) <- asSOAC op,
+    Just (used_arrs, map_lam') <- remove map_lam arrs =
+      Simplify . auxing aux . letBind pat . Op $
+        soacOp (Scatter w used_arrs map_lam' dests)
   where
     used_in_body map_lam = freeIn $ lambdaBody map_lam
     usedInput map_lam (param, _) = paramName param `nameIn` used_in_body map_lam
+    remove map_lam arrs =
+      let (used, unused) = partition (usedInput map_lam) (zip (lambdaParams map_lam) arrs)
+          (used_params, used_arrs) = unzip used
+          map_lam' = map_lam {lambdaParams = used_params}
+       in if null unused then Nothing else Just (used_arrs, map_lam')
 removeUnusedSOACInput _ _ _ _ = Skip
 
 removeDeadMapping :: BottomUpRuleOp (Wise SOACS)
diff --git a/src/Futhark/Internalise/Defunctionalise.hs b/src/Futhark/Internalise/Defunctionalise.hs
--- a/src/Futhark/Internalise/Defunctionalise.hs
+++ b/src/Futhark/Internalise/Defunctionalise.hs
@@ -842,7 +842,7 @@
           params_for_rettype = params ++ svParams sv1 ++ svParams sv2
           svParams (LambdaSV sv_pat _ _ _) = [sv_pat]
           svParams _ = []
-          rettype = buildRetType closure_env params_for_rettype (unRetType e0_t) $ typeOf e0'
+          lifted_rettype = buildRetType closure_env params_for_rettype (unRetType e0_t) $ typeOf e0'
 
           already_bound =
             globals
@@ -872,7 +872,7 @@
       fname <- newNameFromString $ liftedName (0 :: Int) e1
       liftValDec
         fname
-        (RetType [] $ toStruct rettype)
+        (RetType [] $ toStruct lifted_rettype)
         (dims ++ more_dims ++ missing_dims)
         params'
         e0'
@@ -886,39 +886,30 @@
               ( Info
                   ( Scalar . Arrow mempty Unnamed t1 . RetType [] $
                       Scalar . Arrow mempty Unnamed t2 $
-                        RetType [] rettype
+                        RetType [] lifted_rettype
                   )
               )
               loc
 
-          -- FIXME: what if this application returns both a function
-          -- and a value?
-          callret
-            | orderZero ret = AppRes ret ext
-            | otherwise = AppRes rettype ext
+          callret = AppRes (combineTypeShapes ret lifted_rettype) ext
 
+          innercallret =
+            AppRes
+              (Scalar $ Arrow mempty Unnamed t2 $ RetType [] lifted_rettype)
+              []
+
       pure
-        ( Parens
-            ( AppExp
-                ( Apply
-                    ( AppExp
-                        (Apply fname'' e1' (Info (Observe, Nothing)) loc)
-                        ( Info $
-                            AppRes
-                              ( Scalar $
-                                  Arrow mempty Unnamed t2 $
-                                    RetType [] rettype
-                              )
-                              []
-                        )
-                    )
-                    e2'
-                    d
-                    loc
+        ( AppExp
+            ( Apply
+                ( AppExp
+                    (Apply fname'' e1' (Info (Observe, Nothing)) loc)
+                    (Info innercallret)
                 )
-                (Info callret)
+                e2'
+                d
+                loc
             )
-            mempty,
+            (Info callret),
           sv
         )
 
@@ -1182,6 +1173,7 @@
   | otherwise =
       error $ "matchPatSV: missing constructor in type: " ++ prettyString c1
 matchPatSV pat (Dynamic t) = matchPatSV pat $ svFromType t
+matchPatSV pat (HoleSV t _) = matchPatSV pat $ svFromType t
 matchPatSV pat sv =
   error $
     "Tried to match pattern "
@@ -1235,6 +1227,7 @@
 updatePat (PatConstr c1 _ ps loc) (Dynamic t) =
   PatConstr c1 (Info t) ps loc
 updatePat pat (Dynamic t) = updatePat pat (svFromType t)
+updatePat pat (HoleSV t _) = updatePat pat (svFromType t)
 updatePat pat sv =
   error $
     "Tried to update pattern "
diff --git a/src/Futhark/Optimise/Fusion.hs b/src/Futhark/Optimise/Fusion.hs
--- a/src/Futhark/Optimise/Fusion.hs
+++ b/src/Futhark/Optimise/Fusion.hs
@@ -179,29 +179,27 @@
   | otherwise = pure dg
 
 hFuseContexts :: DepContext -> DepContext -> FusionM (Maybe DepContext)
-hFuseContexts
-  c1@(_, _, nodeT1, _)
-  c2@(_, _, nodeT2, _) = do
-    fres <- hFuseNodeT nodeT1 nodeT2
-    case fres of
-      Just nodeT -> pure $ Just (mergedContext nodeT c1 c2)
-      Nothing -> pure Nothing
+hFuseContexts c1 c2 = do
+  let (_, _, nodeT1, _) = c1
+      (_, _, nodeT2, _) = c2
+  fres <- hFuseNodeT nodeT1 nodeT2
+  case fres of
+    Just nodeT -> pure $ Just (mergedContext nodeT c1 c2)
+    Nothing -> pure Nothing
 
 vFuseContexts :: [EdgeT] -> [VName] -> DepContext -> DepContext -> FusionM (Maybe DepContext)
-vFuseContexts
-  edgs
-  infusable
-  c1@(i1, n1, nodeT1, o1)
-  c2@(_i2, n2, nodeT2, o2) = do
-    fres <-
-      vFuseNodeT
-        edgs
-        infusable
-        (nodeT1, map fst $ filter ((/=) n2 . snd) i1, map fst o1)
-        (nodeT2, map fst $ filter ((/=) n1 . snd) o2)
-    case fres of
-      Just nodeT -> pure $ Just (mergedContext nodeT c1 c2)
-      Nothing -> pure Nothing
+vFuseContexts edgs infusable c1 c2 = do
+  let (i1, n1, nodeT1, o1) = c1
+      (_i2, n2, nodeT2, o2) = c2
+  fres <-
+    vFuseNodeT
+      edgs
+      infusable
+      (nodeT1, map fst $ filter ((/=) n2 . snd) i1, map fst o1)
+      (nodeT2, map fst $ filter ((/=) n1 . snd) o2)
+  case fres of
+    Just nodeT -> pure $ Just (mergedContext nodeT c1 c2)
+    Nothing -> pure Nothing
 
 makeCopiesOfFusedExcept ::
   (LocalScope SOACS m, MonadFreshNames m) =>
@@ -284,7 +282,7 @@
     ok <- okToFuseProducer soac1
     r <-
       if ok && ots1 == mempty
-        then TF.attemptFusion preserve (patNames pats1) soac1 ker
+        then TF.attemptFusion TF.Vertical preserve (patNames pats1) soac1 ker
         else pure Nothing
     case r of
       Just ker' -> do
@@ -319,7 +317,7 @@
                 TF.fsOutNames = patNames pats2
               }
           preserve = namesFromList $ patNames pats1
-      r <- TF.attemptFusion preserve (patNames pats1) soac1 ker
+      r <- TF.attemptFusion TF.Horizontal preserve (patNames pats1) soac1 ker
       case r of
         Just ker' -> do
           let pats2' =
diff --git a/src/Futhark/Optimise/Fusion/TryFusion.hs b/src/Futhark/Optimise/Fusion/TryFusion.hs
--- a/src/Futhark/Optimise/Fusion/TryFusion.hs
+++ b/src/Futhark/Optimise/Fusion/TryFusion.hs
@@ -8,6 +8,7 @@
 -- various ways.
 module Futhark.Optimise.Fusion.TryFusion
   ( FusedSOAC (..),
+    Mode (..),
     attemptFusion,
   )
 where
@@ -92,17 +93,18 @@
 setInputs inps ker = ker {fsSOAC = inps `SOAC.setInputs` fsSOAC ker}
 
 tryOptimizeSOAC ::
+  Mode ->
   Names ->
   [VName] ->
   SOAC ->
   FusedSOAC ->
   TryFusion FusedSOAC
-tryOptimizeSOAC unfus_nms outVars soac ker = do
+tryOptimizeSOAC mode unfus_nms outVars soac ker = do
   (soac', ots) <- optimizeSOAC Nothing soac mempty
   let ker' = map (addInitialTransformIfRelevant ots) (inputs ker) `setInputs` ker
       outIdents = zipWith Ident outVars $ SOAC.typeOf soac'
       ker'' = fixInputTypes outIdents ker'
-  applyFusionRules unfus_nms outVars soac' ker''
+  applyFusionRules mode unfus_nms outVars soac' ker''
   where
     addInitialTransformIfRelevant ots inp
       | SOAC.inputArray inp `elem` outVars =
@@ -111,32 +113,34 @@
           inp
 
 tryOptimizeKernel ::
+  Mode ->
   Names ->
   [VName] ->
   SOAC ->
   FusedSOAC ->
   TryFusion FusedSOAC
-tryOptimizeKernel unfus_nms outVars soac ker = do
+tryOptimizeKernel mode unfus_nms outVars soac ker = do
   ker' <- optimizeKernel (Just outVars) ker
-  applyFusionRules unfus_nms outVars soac ker'
+  applyFusionRules mode unfus_nms outVars soac ker'
 
 tryExposeInputs ::
+  Mode ->
   Names ->
   [VName] ->
   SOAC ->
   FusedSOAC ->
   TryFusion FusedSOAC
-tryExposeInputs unfus_nms outVars soac ker = do
+tryExposeInputs mode unfus_nms outVars soac ker = do
   (ker', ots) <- exposeInputs outVars ker
   if SOAC.nullTransforms ots
-    then fuseSOACwithKer unfus_nms outVars soac ker'
+    then fuseSOACwithKer mode unfus_nms outVars soac ker'
     else do
       guard $ unfus_nms == mempty
       (soac', ots') <- pullOutputTransforms soac ots
       let outIdents = zipWith Ident outVars $ SOAC.typeOf soac'
           ker'' = fixInputTypes outIdents ker'
       if SOAC.nullTransforms ots'
-        then applyFusionRules unfus_nms outVars soac' ker''
+        then applyFusionRules mode unfus_nms outVars soac' ker''
         else fail "tryExposeInputs could not pull SOAC transforms"
 
 fixInputTypes :: [Ident] -> FusedSOAC -> FusedSOAC
@@ -151,20 +155,27 @@
     fixInputType inp = inp
 
 applyFusionRules ::
+  Mode ->
   Names ->
   [VName] ->
   SOAC ->
   FusedSOAC ->
   TryFusion FusedSOAC
-applyFusionRules unfus_nms outVars soac ker =
-  tryOptimizeSOAC unfus_nms outVars soac ker
-    <|> tryOptimizeKernel unfus_nms outVars soac ker
-    <|> fuseSOACwithKer unfus_nms outVars soac ker
-    <|> tryExposeInputs unfus_nms outVars soac ker
+applyFusionRules mode unfus_nms outVars soac ker =
+  tryOptimizeSOAC mode unfus_nms outVars soac ker
+    <|> tryOptimizeKernel mode unfus_nms outVars soac ker
+    <|> fuseSOACwithKer mode unfus_nms outVars soac ker
+    <|> tryExposeInputs mode unfus_nms outVars soac ker
 
+-- | Whether we are doing horizontal or vertical fusion.  Note that
+-- vertical also includes "diagonal" fusion, where some producer
+-- results are also produced by the final SOAC.
+data Mode = Horizontal | Vertical
+
 -- | Attempt fusing the producer into the consumer.
 attemptFusion ::
   (HasScope SOACS m, MonadFreshNames m) =>
+  Mode ->
   -- | Outputs of the producer that should still be output by the
   -- fusion result (corresponding to "diagonal fusion").
   Names ->
@@ -173,9 +184,9 @@
   SOAC ->
   FusedSOAC ->
   m (Maybe FusedSOAC)
-attemptFusion unfus_nms outVars soac ker = do
+attemptFusion mode unfus_nms outVars soac ker = do
   scope <- askScope
-  tryFusion (applyFusionRules unfus_nms outVars soac ker) scope
+  tryFusion (applyFusionRules mode unfus_nms outVars soac ker) scope
 
 -- | Check that the consumer does not use any scan or reduce results.
 scremaFusionOK :: ([VName], [VName]) -> FusedSOAC -> Bool
@@ -190,17 +201,17 @@
 
 -- | The brain of this module: Fusing a SOAC with a Kernel.
 fuseSOACwithKer ::
+  Mode ->
   Names ->
   [VName] ->
   SOAC ->
   FusedSOAC ->
   TryFusion FusedSOAC
-fuseSOACwithKer unfus_set outVars soac_p ker = do
+fuseSOACwithKer mode unfus_set outVars soac_p ker = do
   -- We are fusing soac_p into soac_c, i.e, the output of soac_p is going
   -- into soac_c.
   let soac_c = fsSOAC ker
       inp_p_arr = SOAC.inputs soac_p
-      horizFuse = unfus_set /= mempty
       inp_c_arr = SOAC.inputs soac_c
       lam_p = SOAC.lambda soac_p
       lam_c = SOAC.lambda soac_c
@@ -242,13 +253,14 @@
             res_lam' = res_lam {lambdaReturnType = lambdaReturnType res_lam ++ extra_rtps}
          in (extra_nms, res_lam', new_inp)
 
-  when (horizFuse && not (SOAC.nullTransforms $ fsOutputTransform ker)) $
-    fail "Horizontal fusion is invalid in the presence of output transforms."
-
-  case (soac_c, soac_p) of
+  case (soac_c, soac_p, mode) of
     _ | SOAC.width soac_p /= SOAC.width soac_c -> fail "SOAC widths must match."
+    (_, _, Horizontal)
+      | not (SOAC.nullTransforms $ fsOutputTransform ker) ->
+          fail "Horizontal fusion is invalid in the presence of output transforms."
     ( SOAC.Screma _ (ScremaForm scans_c reds_c _) _,
-      SOAC.Screma _ (ScremaForm scans_p reds_p _) _
+      SOAC.Screma _ (ScremaForm scans_p reds_p _) _,
+      _
       )
         | scremaFusionOK (splitAt (Futhark.scanResults scans_p + Futhark.redResults reds_p) outVars) ker -> do
             let red_nes_p = concatMap redNeutral reds_p
@@ -296,7 +308,8 @@
     -- checking that the Scatter is not writing to any array used in
     -- the Map.
     ( SOAC.Scatter _len _lam _ivs dests,
-      SOAC.Screma _ form _
+      SOAC.Screma _ form _,
+      _
       )
         | isJust $ isMapSOAC form,
           -- 1. all arrays produced by the map are ONLY used (consumed)
@@ -314,7 +327,8 @@
     -- checking that the Hist is not writing to any array used in
     -- the Map.
     ( SOAC.Hist _ ops _ _,
-      SOAC.Screma _ form _
+      SOAC.Screma _ form _,
+      _
       )
         | isJust $ isMapSOAC form,
           -- 1. all arrays produced by the map are ONLY used (consumed)
@@ -328,66 +342,70 @@
 
     -- Hist-Hist fusion
     ( SOAC.Hist _ ops_c _ _,
-      SOAC.Hist _ ops_p _ _
-      )
-        | horizFuse -> do
-            let p_num_buckets = length ops_p
-                c_num_buckets = length ops_c
-                (body_p, body_c) = (lambdaBody lam_p, lambdaBody lam_c)
-                body' =
-                  Body
-                    { bodyDec = bodyDec body_p, -- body_p and body_c have the same decorations
-                      bodyStms = bodyStms body_p <> bodyStms body_c,
-                      bodyResult =
-                        take c_num_buckets (bodyResult body_c)
-                          ++ take p_num_buckets (bodyResult body_p)
-                          ++ drop c_num_buckets (bodyResult body_c)
-                          ++ drop p_num_buckets (bodyResult body_p)
-                    }
-                lam' =
-                  Lambda
-                    { lambdaParams = lambdaParams lam_c ++ lambdaParams lam_p,
-                      lambdaBody = body',
-                      lambdaReturnType =
-                        replicate (c_num_buckets + p_num_buckets) (Prim int64)
-                          ++ drop c_num_buckets (lambdaReturnType lam_c)
-                          ++ drop p_num_buckets (lambdaReturnType lam_p)
-                    }
-            success (fsOutNames ker ++ returned_outvars) $
-              SOAC.Hist w (ops_c <> ops_p) lam' (inp_c_arr <> inp_p_arr)
+      SOAC.Hist _ ops_p _ _,
+      Horizontal
+      ) -> do
+        let p_num_buckets = length ops_p
+            c_num_buckets = length ops_c
+            (body_p, body_c) = (lambdaBody lam_p, lambdaBody lam_c)
+            body' =
+              Body
+                { bodyDec = bodyDec body_p, -- body_p and body_c have the same decorations
+                  bodyStms = bodyStms body_p <> bodyStms body_c,
+                  bodyResult =
+                    take c_num_buckets (bodyResult body_c)
+                      ++ take p_num_buckets (bodyResult body_p)
+                      ++ drop c_num_buckets (bodyResult body_c)
+                      ++ drop p_num_buckets (bodyResult body_p)
+                }
+            lam' =
+              Lambda
+                { lambdaParams = lambdaParams lam_c ++ lambdaParams lam_p,
+                  lambdaBody = body',
+                  lambdaReturnType =
+                    replicate (c_num_buckets + p_num_buckets) (Prim int64)
+                      ++ drop c_num_buckets (lambdaReturnType lam_c)
+                      ++ drop p_num_buckets (lambdaReturnType lam_p)
+                }
+        success (fsOutNames ker ++ returned_outvars) $
+          SOAC.Hist w (ops_c <> ops_p) lam' (inp_c_arr <> inp_p_arr)
 
     -- Scatter-write fusion.
     ( SOAC.Scatter _len_c _lam_c ivs_c as_c,
-      SOAC.Scatter _len_p _lam_p ivs_p as_p
-      )
-        | horizFuse -> do
-            let zipW as_xs xs as_ys ys = xs_indices ++ ys_indices ++ xs_vals ++ ys_vals
-                  where
-                    (xs_indices, xs_vals) = splitScatterResults as_xs xs
-                    (ys_indices, ys_vals) = splitScatterResults as_ys ys
-            let (body_p, body_c) = (lambdaBody lam_p, lambdaBody lam_c)
-            let body' =
-                  Body
-                    { bodyDec = bodyDec body_p, -- body_p and body_c have the same decorations
-                      bodyStms = bodyStms body_p <> bodyStms body_c,
-                      bodyResult = zipW as_c (bodyResult body_c) as_p (bodyResult body_p)
-                    }
-            let lam' =
-                  Lambda
-                    { lambdaParams = lambdaParams lam_c ++ lambdaParams lam_p,
-                      lambdaBody = body',
-                      lambdaReturnType = zipW as_c (lambdaReturnType lam_c) as_p (lambdaReturnType lam_p)
-                    }
-            success (fsOutNames ker ++ returned_outvars) $
-              SOAC.Scatter w lam' (ivs_c ++ ivs_p) (as_c ++ as_p)
-    (SOAC.Scatter {}, _) ->
-      fail "Cannot fuse a write with anything else than a write or a map"
-    (_, SOAC.Scatter {}) ->
-      fail "Cannot fuse a write with anything else than a write or a map"
+      SOAC.Scatter
+        _len_p
+        _lam_p
+        ivs_p
+        as_p,
+      Horizontal
+      ) -> do
+        let zipW as_xs xs as_ys ys = xs_indices ++ ys_indices ++ xs_vals ++ ys_vals
+              where
+                (xs_indices, xs_vals) = splitScatterResults as_xs xs
+                (ys_indices, ys_vals) = splitScatterResults as_ys ys
+        let (body_p, body_c) = (lambdaBody lam_p, lambdaBody lam_c)
+        let body' =
+              Body
+                { bodyDec = bodyDec body_p, -- body_p and body_c have the same decorations
+                  bodyStms = bodyStms body_p <> bodyStms body_c,
+                  bodyResult = zipW as_c (bodyResult body_c) as_p (bodyResult body_p)
+                }
+        let lam' =
+              Lambda
+                { lambdaParams = lambdaParams lam_c ++ lambdaParams lam_p,
+                  lambdaBody = body',
+                  lambdaReturnType = zipW as_c (lambdaReturnType lam_c) as_p (lambdaReturnType lam_p)
+                }
+        success (fsOutNames ker ++ returned_outvars) $
+          SOAC.Scatter w lam' (ivs_c ++ ivs_p) (as_c ++ as_p)
+    (SOAC.Scatter {}, _, _) ->
+      fail "Cannot fuse a scatter with anything else than a scatter or a map"
+    (_, SOAC.Scatter {}, _) ->
+      fail "Cannot fuse a scatter with anything else than a scatter or a map"
     ----------------------------
     -- Stream-Stream Fusions: --
     ----------------------------
-    (SOAC.Stream {}, SOAC.Stream {}) -> do
+    (SOAC.Stream {}, SOAC.Stream {}, _) -> do
       -- fuse two SEQUENTIAL streams
       (res_nms, res_stream) <- fuseStreamHelper (fsOutNames ker) unfus_set outVars outPairs soac_c soac_p
       success res_nms res_stream
@@ -400,18 +418,19 @@
     ---   we could run in an infinite recursion, i.e., repeatedly   ---
     ---   fusing map o scan into an infinity of Stream levels!      ---
     -------------------------------------------------------------------
-    (SOAC.Stream {}, _) -> do
+    (SOAC.Stream {}, _, _) -> do
       -- If this rule is matched then soac_p is NOT a stream.
       -- To fuse a stream kernel, we transform soac_p to a stream, which
       -- borrows the sequential/parallel property of the soac_c Stream,
       -- and recursively perform stream-stream fusion.
       (soac_p', newacc_ids) <- SOAC.soacToStream soac_p
       fuseSOACwithKer
+        mode
         (namesFromList (map identName newacc_ids) <> unfus_set)
         (map identName newacc_ids ++ outVars)
         soac_p'
         ker
-    (_, SOAC.Screma _ form _) | Just _ <- Futhark.isScanomapSOAC form -> do
+    (_, SOAC.Screma _ form _, _) | Just _ <- Futhark.isScanomapSOAC form -> do
       -- A Scan soac can be currently only fused as a (sequential) stream,
       -- hence it is first translated to a (sequential) Stream and then
       -- fusion with a kernel is attempted.
@@ -419,12 +438,13 @@
       if soac_p' /= soac_p
         then
           fuseSOACwithKer
+            mode
             (namesFromList (map identName newacc_ids) <> unfus_set)
             (map identName newacc_ids ++ outVars)
             soac_p'
             ker
         else fail "SOAC could not be turned into stream."
-    (_, SOAC.Stream {}) -> do
+    (_, SOAC.Stream {}, _) -> do
       -- If it reached this case then soac_c is NOT a Stream kernel,
       -- hence transform the kernel's soac to a stream and attempt
       -- stream-stream fusion recursivelly.
@@ -434,6 +454,7 @@
       if soac_c' /= soac_c
         then
           fuseSOACwithKer
+            mode
             (namesFromList (map identName newacc_ids) <> unfus_set)
             outVars
             soac_p
diff --git a/src/Futhark/Optimise/Simplify/Rules.hs b/src/Futhark/Optimise/Simplify/Rules.hs
--- a/src/Futhark/Optimise/Simplify/Rules.hs
+++ b/src/Futhark/Optimise/Simplify/Rules.hs
@@ -38,7 +38,8 @@
 topDownRules :: BuilderOps rep => [TopDownRule rep]
 topDownRules =
   [ RuleGeneric constantFoldPrimFun,
-    RuleGeneric withAccTopDown
+    RuleGeneric withAccTopDown,
+    RuleGeneric emptyArrayToScratch
   ]
 
 bottomUpRules :: (BuilderOps rep, TraverseOpStms rep) => [BottomUpRule rep]
@@ -111,6 +112,18 @@
     isConst (Constant v) = Just v
     isConst _ = Nothing
 constantFoldPrimFun _ _ = Skip
+
+-- | If an expression produces an array with a constant zero anywhere
+-- in its shape, just turn that into a Scratch.
+emptyArrayToScratch :: BuilderOps rep => TopDownRuleGeneric rep
+emptyArrayToScratch _ (Let pat@(Pat [pe]) aux e)
+  | Just (pt, shape) <- isEmptyArray $ patElemType pe,
+    not $ isScratch e =
+      Simplify $ auxing aux $ letBind pat $ BasicOp $ Scratch pt $ shapeDims shape
+  where
+    isScratch (BasicOp Scratch {}) = True
+    isScratch _ = False
+emptyArrayToScratch _ _ = Skip
 
 simplifyIndex :: BuilderOps rep => BottomUpRuleBasicOp rep
 simplifyIndex (vtable, used) pat@(Pat [pe]) (StmAux cs attrs _) (Index idd inds)
diff --git a/src/Futhark/Optimise/Simplify/Rules/BasicOp.hs b/src/Futhark/Optimise/Simplify/Rules/BasicOp.hs
--- a/src/Futhark/Optimise/Simplify/Rules/BasicOp.hs
+++ b/src/Futhark/Optimise/Simplify/Rules/BasicOp.hs
@@ -10,6 +10,7 @@
 import Control.Monad
 import Data.List (find, foldl', isSuffixOf, sort)
 import Data.List.NonEmpty (NonEmpty (..))
+import Data.Maybe (isNothing)
 import Futhark.Analysis.PrimExp.Convert
 import Futhark.Analysis.SymbolTable qualified as ST
 import Futhark.Construct
@@ -49,11 +50,9 @@
   m VName
 fromConcatArg t (ArgArrayLit ses, cs) =
   certifying cs $ letExp "concat_lit" $ BasicOp $ ArrayLit ses $ rowType t
-fromConcatArg elem_type (ArgReplicate ws se, cs) = do
-  let elem_shape = arrayShape elem_type
-  certifying cs $ do
-    w <- letSubExp "concat_rep_w" =<< toExp (sum $ map pe64 ws)
-    letExp "concat_rep" $ BasicOp $ Replicate (setDim 0 elem_shape w) se
+fromConcatArg _ (ArgReplicate ws se, cs) = certifying cs $ do
+  w <- letSubExp "concat_rep_w" =<< toExp (sum $ map pe64 ws)
+  letExp "concat_rep" $ BasicOp $ Replicate (Shape [w]) se
 fromConcatArg _ (ArgVar v, _) =
   pure v
 
@@ -121,6 +120,12 @@
       Just (Concat j (y :| ys) _, v_cs) | j == i -> (y : ys, v_cs)
       _ -> ([v], mempty)
 
+-- Removing empty arrays from concatenations.
+simplifyConcat (vtable, _) pat aux (Concat i (x :| xs) new_d)
+  | Just ts <- mapM (`ST.lookupType` vtable) $ x : xs,
+    x' : xs' <- map fst $ filter (isNothing . isEmptyArray . snd) $ zip (x : xs) ts,
+    xs' /= xs =
+      Simplify $ auxing aux $ letBind pat $ BasicOp $ Concat i (x' :| xs') new_d
 -- Fusing arguments to the concat when possible.  Only done when
 -- concatenating along the outer dimension for now.
 simplifyConcat (vtable, _) pat aux (Concat 0 (x :| xs) outer_w)
@@ -187,7 +192,7 @@
           True
     arrayFrom _ =
       False
-ruleBasicOp vtable pat aux (Update _ dest is se)
+ruleBasicOp vtable pat aux (Update Unsafe dest is se)
   | Just dest_t <- ST.lookupType dest vtable,
     isFullSlice (arrayShape dest_t) is = Simplify . auxing aux $
       case se of
diff --git a/src/Futhark/Optimise/Simplify/Rules/Loop.hs b/src/Futhark/Optimise/Simplify/Rules/Loop.hs
--- a/src/Futhark/Optimise/Simplify/Rules/Loop.hs
+++ b/src/Futhark/Optimise/Simplify/Rules/Loop.hs
@@ -35,8 +35,7 @@
 
           resIsNecessary ((v, _), _) =
             usedAfterLoop v
-              || paramName v
-              `nameIn` necessaryForReturned
+              || (paramName v `nameIn` necessaryForReturned)
               || referencedInPat v
               || referencedInForm v
 
@@ -153,6 +152,8 @@
               p_init == resSubExp resExp,
               p_res == Var pat_name =
                 True
+            -- (3) It is a statically empty array.
+            | isJust $ isEmptyArray (paramType mergeParam) = True
             | otherwise = False
     checkInvariance
       (_pat_name, (mergeParam, mergeInit), resExp)
diff --git a/src/Futhark/Pass/ExtractKernels/Distribution.hs b/src/Futhark/Pass/ExtractKernels/Distribution.hs
--- a/src/Futhark/Pass/ExtractKernels/Distribution.hs
+++ b/src/Futhark/Pass/ExtractKernels/Distribution.hs
@@ -326,9 +326,12 @@
   let bound_by_stms = namesFromList $ M.keys $ scopeOf stms
       (inner_pat', inner_res', inner_identity_map, inner_expand_target) =
         removeIdentityMappingGeneral bound_by_stms inner_pat inner_res
+      free =
+        (foldMap freeIn stms <> freeIn (map resCerts inner_res))
+          `namesSubtract` bound_by_stms
    in ( DistributionBody
           { distributionTarget = Targets (inner_pat', inner_res') targets,
-            distributionFreeInBody = foldMap freeIn stms `namesSubtract` bound_by_stms,
+            distributionFreeInBody = free,
             distributionIdentityMap = inner_identity_map,
             distributionExpandTarget = inner_expand_target
           },
diff --git a/src/Futhark/Util/Pretty.hs b/src/Futhark/Util/Pretty.hs
--- a/src/Futhark/Util/Pretty.hs
+++ b/src/Futhark/Util/Pretty.hs
@@ -103,7 +103,10 @@
 -- | Convert a 'Doc' to text.  Thsi ignores any annotations (i.e. it
 -- will be non-coloured output).
 docText :: Doc a -> T.Text
-docText = Prettyprinter.Render.Text.renderStrict . layoutSmart defaultLayoutOptions
+docText = Prettyprinter.Render.Text.renderStrict . layouter
+  where
+    layouter =
+      layoutSmart defaultLayoutOptions {layoutPageWidth = Unbounded}
 
 -- | Prettyprint a value to a 'Text' on a single line.
 prettyTextOneLine :: Pretty a => a -> Text
diff --git a/src/Language/Futhark/Interpreter.hs b/src/Language/Futhark/Interpreter.hs
--- a/src/Language/Futhark/Interpreter.hs
+++ b/src/Language/Futhark/Interpreter.hs
@@ -76,7 +76,7 @@
     BreakNaN
 
 data ExtOp a
-  = ExtOpTrace String String a
+  = ExtOpTrace T.Text (Doc ()) a
   | ExtOpBreak Loc BreakReason (NE.NonEmpty StackFrame) a
   | ExtOpError InterpreterError
 
@@ -330,9 +330,9 @@
   liftF . ExtOpError . InterpreterError $
     "Error at\n" <> prettyStacktrace 0 ss <> s
 
-trace :: String -> Value -> EvalM ()
+trace :: T.Text -> Value -> EvalM ()
 trace w v = do
-  liftF $ ExtOpTrace w (T.unpack $ docText $ oneLine $ prettyValue v) ()
+  liftF $ ExtOpTrace w (prettyValue v) ()
 
 typeCheckerEnv :: Env -> T.Env
 typeCheckerEnv env =
@@ -995,11 +995,11 @@
   eval env e
 eval env (Attr (AttrAtom (AtomName "trace") _) e loc) = do
   v <- eval env e
-  trace (locStr (locOf loc)) v
+  trace (locText (locOf loc)) v
   pure v
 eval env (Attr (AttrComp "trace" [AttrAtom (AtomName tag) _] _) e _) = do
   v <- eval env e
-  trace (nameToString tag) v
+  trace (nameToText tag) v
   pure v
 eval env (Attr _ e _) =
   eval env e
diff --git a/src/Language/Futhark/Pretty.hs b/src/Language/Futhark/Pretty.hs
--- a/src/Language/Futhark/Pretty.hs
+++ b/src/Language/Futhark/Pretty.hs
@@ -308,6 +308,9 @@
           "@" <> parens (align $ pretty t')
     _ -> mempty
 
+prettyAttr :: Pretty a => a -> Doc ann
+prettyAttr attr = "#[" <> pretty attr <> "]"
+
 prettyExp :: (Eq vn, IsName vn, Annot f) => Int -> ExpBase f vn -> Doc a
 prettyExp _ (Var name t _) = pretty name <> prettyInst t
 prettyExp _ (Hole t _) = "???" <> prettyInst t
@@ -366,7 +369,7 @@
   parens $ "." <> brackets (commasep (map pretty idxs))
 prettyExp _ (Constr n cs _ _) = "#" <> pretty n <+> sep (map pretty cs)
 prettyExp _ (Attr attr e _) =
-  "#[" <> pretty attr <> "]" </> prettyExp (-1) e
+  prettyAttr attr </> prettyExp (-1) e
 prettyExp i (AppExp e _) = prettyAppExp i e
 
 instance (Eq vn, IsName vn, Annot f) => Pretty (ExpBase f vn) where
@@ -467,7 +470,7 @@
 
 instance (Eq vn, IsName vn, Annot f) => Pretty (ValBindBase f vn) where
   pretty (ValBind entry name retdecl rettype tparams args body _ attrs _) =
-    mconcat (map ((<> line) . pretty) attrs)
+    mconcat (map ((<> line) . prettyAttr) attrs)
       <> fun
       <+> align
         ( sep
diff --git a/src/Language/Futhark/Primitive.hs b/src/Language/Futhark/Primitive.hs
--- a/src/Language/Futhark/Primitive.hs
+++ b/src/Language/Futhark/Primitive.hs
@@ -1294,78 +1294,22 @@
       i16 "popc16" $ IntValue . Int32Value . fromIntegral . popCount,
       i32 "popc32" $ IntValue . Int32Value . fromIntegral . popCount,
       i64 "popc64" $ IntValue . Int32Value . fromIntegral . popCount,
-      ( "mad_hi8",
-        ( [IntType Int8, IntType Int8, IntType Int8],
-          IntType Int8,
-          \case
-            [IntValue (Int8Value a), IntValue (Int8Value b), IntValue (Int8Value c)] ->
-              Just $ IntValue . Int8Value $ mad_hi8 (Int8Value a) (Int8Value b) c
-            _ -> Nothing
-        )
-      ),
-      ( "mad_hi16",
-        ( [IntType Int16, IntType Int16, IntType Int16],
-          IntType Int16,
-          \case
-            [IntValue (Int16Value a), IntValue (Int16Value b), IntValue (Int16Value c)] ->
-              Just $ IntValue . Int16Value $ mad_hi16 (Int16Value a) (Int16Value b) c
-            _ -> Nothing
-        )
-      ),
-      ( "mad_hi32",
-        ( [IntType Int32, IntType Int32, IntType Int32],
-          IntType Int32,
-          \case
-            [IntValue (Int32Value a), IntValue (Int32Value b), IntValue (Int32Value c)] ->
-              Just $ IntValue . Int32Value $ mad_hi32 (Int32Value a) (Int32Value b) c
-            _ -> Nothing
-        )
-      ),
-      ( "mad_hi64",
-        ( [IntType Int64, IntType Int64, IntType Int64],
-          IntType Int64,
-          \case
-            [IntValue (Int64Value a), IntValue (Int64Value b), IntValue (Int64Value c)] ->
-              Just $ IntValue . Int64Value $ mad_hi64 (Int64Value a) (Int64Value b) c
-            _ -> Nothing
-        )
-      ),
-      ( "mul_hi8",
-        ( [IntType Int8, IntType Int8],
-          IntType Int8,
-          \case
-            [IntValue (Int8Value a), IntValue (Int8Value b)] ->
-              Just $ IntValue . Int8Value $ mul_hi8 (Int8Value a) (Int8Value b)
-            _ -> Nothing
-        )
-      ),
-      ( "mul_hi16",
-        ( [IntType Int16, IntType Int16],
-          IntType Int16,
-          \case
-            [IntValue (Int16Value a), IntValue (Int16Value b)] ->
-              Just $ IntValue . Int16Value $ mul_hi16 (Int16Value a) (Int16Value b)
-            _ -> Nothing
-        )
-      ),
-      ( "mul_hi32",
-        ( [IntType Int32, IntType Int32],
-          IntType Int32,
-          \case
-            [IntValue (Int32Value a), IntValue (Int32Value b)] ->
-              Just $ IntValue . Int32Value $ mul_hi32 (Int32Value a) (Int32Value b)
-            _ -> Nothing
-        )
-      ),
-      ( "mul_hi64",
-        ( [IntType Int64, IntType Int64],
-          IntType Int64,
-          \case
-            [IntValue (Int64Value a), IntValue (Int64Value b)] ->
-              Just $ IntValue . Int64Value $ mul_hi64 (Int64Value a) (Int64Value b)
-            _ -> Nothing
-        )
-      ),
+      i8_3 "umad_hi8" umad_hi8,
+      i16_3 "umad_hi16" umad_hi16,
+      i32_3 "umad_hi32" umad_hi32,
+      i64_3 "umad_hi64" umad_hi64,
+      i8_2 "umul_hi8" umul_hi8,
+      i16_2 "umul_hi16" umul_hi16,
+      i32_2 "umul_hi32" umul_hi32,
+      i64_2 "umul_hi64" umul_hi64,
+      i8_3 "smad_hi8" smad_hi8,
+      i16_3 "smad_hi16" smad_hi16,
+      i32_3 "smad_hi32" smad_hi32,
+      i64_3 "smad_hi64" smad_hi64,
+      i8_2 "smul_hi8" smul_hi8,
+      i16_2 "smul_hi16" smul_hi16,
+      i32_2 "smul_hi32" smul_hi32,
+      i64_2 "smul_hi64" smul_hi64,
       --
       ( "atan2_16",
         ( [FloatType Float16, FloatType Float16],
@@ -1542,63 +1486,33 @@
     f16 s f = (s, ([FloatType Float16], FloatType Float16, f16PrimFun f))
     f32 s f = (s, ([FloatType Float32], FloatType Float32, f32PrimFun f))
     f64 s f = (s, ([FloatType Float64], FloatType Float64, f64PrimFun f))
-    f16_2 s f =
-      ( s,
-        ( [FloatType Float16, FloatType Float16],
-          FloatType Float16,
-          f16PrimFun2 f
-        )
-      )
-    f32_2 s f =
-      ( s,
-        ( [FloatType Float32, FloatType Float32],
-          FloatType Float32,
-          f32PrimFun2 f
-        )
-      )
-    f64_2 s f =
-      ( s,
-        ( [FloatType Float64, FloatType Float64],
-          FloatType Float64,
-          f64PrimFun2 f
-        )
-      )
-    f16_3 s f =
-      ( s,
-        ( [FloatType Float16, FloatType Float16, FloatType Float16],
-          FloatType Float16,
-          f16PrimFun3 f
-        )
-      )
-    f32_3 s f =
-      ( s,
-        ( [FloatType Float32, FloatType Float32, FloatType Float32],
-          FloatType Float32,
-          f32PrimFun3 f
-        )
-      )
-    f64_3 s f =
-      ( s,
-        ( [FloatType Float64, FloatType Float64, FloatType Float64],
-          FloatType Float64,
-          f64PrimFun3 f
-        )
-      )
+    t_2 t s f = (s, ([t, t], t, f))
+    t_3 t s f = (s, ([t, t, t], t, f))
+    f16_2 s f = t_2 (FloatType Float16) s (f16PrimFun2 f)
+    f32_2 s f = t_2 (FloatType Float32) s (f32PrimFun2 f)
+    f64_2 s f = t_2 (FloatType Float64) s (f64PrimFun2 f)
+    f16_3 s f = t_3 (FloatType Float16) s (f16PrimFun3 f)
+    f32_3 s f = t_3 (FloatType Float32) s (f32PrimFun3 f)
+    f64_3 s f = t_3 (FloatType Float64) s (f64PrimFun3 f)
+    i8_2 s f = t_2 (IntType Int8) s (i8PrimFun2 f)
+    i16_2 s f = t_2 (IntType Int16) s (i16PrimFun2 f)
+    i32_2 s f = t_2 (IntType Int32) s (i32PrimFun2 f)
+    i64_2 s f = t_2 (IntType Int64) s (i64PrimFun2 f)
+    i8_3 s f = t_3 (IntType Int8) s (i8PrimFun3 f)
+    i16_3 s f = t_3 (IntType Int16) s (i16PrimFun3 f)
+    i32_3 s f = t_3 (IntType Int32) s (i32PrimFun3 f)
+    i64_3 s f = t_3 (IntType Int64) s (i64PrimFun3 f)
 
-    i8PrimFun f [IntValue (Int8Value x)] =
-      Just $ f x
+    i8PrimFun f [IntValue (Int8Value x)] = Just $ f x
     i8PrimFun _ _ = Nothing
 
-    i16PrimFun f [IntValue (Int16Value x)] =
-      Just $ f x
+    i16PrimFun f [IntValue (Int16Value x)] = Just $ f x
     i16PrimFun _ _ = Nothing
 
-    i32PrimFun f [IntValue (Int32Value x)] =
-      Just $ f x
+    i32PrimFun f [IntValue (Int32Value x)] = Just $ f x
     i32PrimFun _ _ = Nothing
 
-    i64PrimFun f [IntValue (Int64Value x)] =
-      Just $ f x
+    i64PrimFun f [IntValue (Int64Value x)] = Just $ f x
     i64PrimFun _ _ = Nothing
 
     f16PrimFun f [FloatValue (Float16Value x)] =
@@ -1664,6 +1578,54 @@
         Just $ FloatValue $ Float64Value $ f a b c
     f64PrimFun3 _ _ = Nothing
 
+    i8PrimFun2
+      f
+      [IntValue (Int8Value a), IntValue (Int8Value b)] =
+        Just $ IntValue $ Int8Value $ f a b
+    i8PrimFun2 _ _ = Nothing
+
+    i16PrimFun2
+      f
+      [IntValue (Int16Value a), IntValue (Int16Value b)] =
+        Just $ IntValue $ Int16Value $ f a b
+    i16PrimFun2 _ _ = Nothing
+
+    i32PrimFun2
+      f
+      [IntValue (Int32Value a), IntValue (Int32Value b)] =
+        Just $ IntValue $ Int32Value $ f a b
+    i32PrimFun2 _ _ = Nothing
+
+    i64PrimFun2
+      f
+      [IntValue (Int64Value a), IntValue (Int64Value b)] =
+        Just $ IntValue $ Int64Value $ f a b
+    i64PrimFun2 _ _ = Nothing
+
+    i8PrimFun3
+      f
+      [IntValue (Int8Value a), IntValue (Int8Value b), IntValue (Int8Value c)] =
+        Just $ IntValue $ Int8Value $ f a b c
+    i8PrimFun3 _ _ = Nothing
+
+    i16PrimFun3
+      f
+      [IntValue (Int16Value a), IntValue (Int16Value b), IntValue (Int16Value c)] =
+        Just $ IntValue $ Int16Value $ f a b c
+    i16PrimFun3 _ _ = Nothing
+
+    i32PrimFun3
+      f
+      [IntValue (Int32Value a), IntValue (Int32Value b), IntValue (Int32Value c)] =
+        Just $ IntValue $ Int32Value $ f a b c
+    i32PrimFun3 _ _ = Nothing
+
+    i64PrimFun3
+      f
+      [IntValue (Int64Value a), IntValue (Int64Value b), IntValue (Int64Value c)] =
+        Just $ IntValue $ Int64Value $ f a b c
+    i64PrimFun3 _ _ = Nothing
+
 -- | Is the given value kind of zero?
 zeroIsh :: PrimValue -> Bool
 zeroIsh (IntValue k) = zeroIshInt k
@@ -1878,38 +1840,74 @@
 prettySigned True (IntType it) = 'u' : drop 1 (prettyString it)
 prettySigned _ t = prettyString t
 
-mul_hi8 :: IntValue -> IntValue -> Int8
-mul_hi8 a b =
-  let a' = intToWord64 a
-      b' = intToWord64 b
+umul_hi8 :: Int8 -> Int8 -> Int8
+umul_hi8 a b =
+  let a' = fromIntegral (fromIntegral a :: Word8) :: Word64
+      b' = fromIntegral (fromIntegral b :: Word8) :: Word64
    in fromIntegral (shiftR (a' * b') 8)
 
-mul_hi16 :: IntValue -> IntValue -> Int16
-mul_hi16 a b =
-  let a' = intToWord64 a
-      b' = intToWord64 b
+umul_hi16 :: Int16 -> Int16 -> Int16
+umul_hi16 a b =
+  let a' = fromIntegral (fromIntegral a :: Word16) :: Word64
+      b' = fromIntegral (fromIntegral b :: Word16) :: Word64
    in fromIntegral (shiftR (a' * b') 16)
 
-mul_hi32 :: IntValue -> IntValue -> Int32
-mul_hi32 a b =
-  let a' = intToWord64 a
-      b' = intToWord64 b
+umul_hi32 :: Int32 -> Int32 -> Int32
+umul_hi32 a b =
+  let a' = fromIntegral (fromIntegral a :: Word32) :: Word64
+      b' = fromIntegral (fromIntegral b :: Word32) :: Word64
    in fromIntegral (shiftR (a' * b') 32)
 
-mul_hi64 :: IntValue -> IntValue -> Int64
-mul_hi64 a b =
-  let a' = (toInteger . intToWord64) a
-      b' = (toInteger . intToWord64) b
+umul_hi64 :: Int64 -> Int64 -> Int64
+umul_hi64 a b =
+  let a' = toInteger (fromIntegral a :: Word64)
+      b' = toInteger (fromIntegral b :: Word64)
    in fromIntegral (shiftR (a' * b') 64)
 
-mad_hi8 :: IntValue -> IntValue -> Int8 -> Int8
-mad_hi8 a b c = mul_hi8 a b + c
+umad_hi8 :: Int8 -> Int8 -> Int8 -> Int8
+umad_hi8 a b c = umul_hi8 a b + c
 
-mad_hi16 :: IntValue -> IntValue -> Int16 -> Int16
-mad_hi16 a b c = mul_hi16 a b + c
+umad_hi16 :: Int16 -> Int16 -> Int16 -> Int16
+umad_hi16 a b c = umul_hi16 a b + c
 
-mad_hi32 :: IntValue -> IntValue -> Int32 -> Int32
-mad_hi32 a b c = mul_hi32 a b + c
+umad_hi32 :: Int32 -> Int32 -> Int32 -> Int32
+umad_hi32 a b c = umul_hi32 a b + c
 
-mad_hi64 :: IntValue -> IntValue -> Int64 -> Int64
-mad_hi64 a b c = mul_hi64 a b + c
+umad_hi64 :: Int64 -> Int64 -> Int64 -> Int64
+umad_hi64 a b c = umul_hi64 a b + c
+
+smul_hi8 :: Int8 -> Int8 -> Int8
+smul_hi8 a b =
+  let a' = fromIntegral a :: Int64
+      b' = fromIntegral b :: Int64
+   in fromIntegral (shiftR (a' * b') 8)
+
+smul_hi16 :: Int16 -> Int16 -> Int16
+smul_hi16 a b =
+  let a' = fromIntegral a :: Int64
+      b' = fromIntegral b :: Int64
+   in fromIntegral (shiftR (a' * b') 16)
+
+smul_hi32 :: Int32 -> Int32 -> Int32
+smul_hi32 a b =
+  let a' = fromIntegral a :: Int64
+      b' = fromIntegral b :: Int64
+   in fromIntegral (shiftR (a' * b') 32)
+
+smul_hi64 :: Int64 -> Int64 -> Int64
+smul_hi64 a b =
+  let a' = toInteger a
+      b' = toInteger b
+   in fromIntegral (shiftR (a' * b') 64)
+
+smad_hi8 :: Int8 -> Int8 -> Int8 -> Int8
+smad_hi8 a b c = smul_hi8 a b + c
+
+smad_hi16 :: Int16 -> Int16 -> Int16 -> Int16
+smad_hi16 a b c = smul_hi16 a b + c
+
+smad_hi32 :: Int32 -> Int32 -> Int32 -> Int32
+smad_hi32 a b c = smul_hi32 a b + c
+
+smad_hi64 :: Int64 -> Int64 -> Int64 -> Int64
+smad_hi64 a b c = smul_hi64 a b + c
diff --git a/src/Language/Futhark/TypeChecker/Unify.hs b/src/Language/Futhark/TypeChecker/Unify.hs
--- a/src/Language/Futhark/TypeChecker/Unify.hs
+++ b/src/Language/Futhark/TypeChecker/Unify.hs
@@ -274,7 +274,7 @@
       aNote $
         prettyName v
           <+> "="
-          <+> mconcat (map ppConstr (M.toList cs))
+          <+> hsep (map ppConstr (M.toList cs))
           <+> "..."
     note (Overloaded ts _) =
       aNote $ prettyName v <+> "must be one of" <+> mconcat (punctuate ", " (map pretty ts))
@@ -351,6 +351,14 @@
 rigidConstraint UnknowableSize {} = True
 rigidConstraint _ = False
 
+unsharedConstructorsMsg :: M.Map Name t -> M.Map Name t -> Doc a
+unsharedConstructorsMsg cs1 cs2 =
+  "Unshared constructors:" <+> commasep (map (("#" <>) . pretty) missing) <> "."
+  where
+    missing =
+      filter (`notElem` M.keys cs1) (M.keys cs2)
+        ++ filter (`notElem` M.keys cs2) (M.keys cs1)
+
 -- | Instantiate existential context in return type.
 instantiateEmptyArrayDims ::
   MonadUnify m =>
@@ -535,12 +543,8 @@
           )
             | M.keys cs == M.keys arg_cs ->
                 unifySharedConstructors onDims usage bound bcs cs arg_cs
-            | otherwise -> do
-                let missing =
-                      filter (`notElem` M.keys arg_cs) (M.keys cs)
-                        ++ filter (`notElem` M.keys cs) (M.keys arg_cs)
-                unifyError usage mempty bcs $
-                  "Unshared constructors:" <+> commasep (map (("#" <>) . pretty) missing) <> "."
+            | otherwise ->
+                unifyError usage mempty bcs $ unsharedConstructorsMsg arg_cs cs
         _
           | t1' == t2' -> pure ()
           | otherwise -> failure
@@ -775,6 +779,8 @@
               modifyConstraints $
                 M.insert vn (lvl, Constraint (RetType ext tp') usage)
               unifySharedConstructors onDims usage bound bcs required_cs ts
+          | otherwise ->
+              unsharedConstructors required_cs ts =<< typeVarNotes vn
         Scalar (TypeVar _ _ (QualName [] v) []) -> do
           case M.lookup v constraints of
             Just (_, HasConstrs _ v_cs _) ->
@@ -794,21 +800,27 @@
             combineConstrs (_, HasConstrs l1 cs1 usage1) (_, HasConstrs l2 cs2 _) =
               (lvl, HasConstrs (l1 `min` l2) (M.union cs1 cs2) usage1)
             combineConstrs hasCs _ = hasCs
-        _ -> noSumType mempty
+        _ -> noSumType =<< typeVarNotes vn
     _ -> link
   where
+    unsharedConstructors cs1 cs2 notes =
+      unifyError
+        usage
+        notes
+        bcs
+        (unsharedConstructorsMsg cs1 cs2)
     noSumType notes =
       unifyError
         usage
         notes
         bcs
-        "Cannot unify a sum type with a non-sum type"
+        "Cannot unify a sum type with a non-sum type."
     noRecordType notes =
       unifyError
         usage
         notes
         bcs
-        "Cannot unify a record type with a non-record type"
+        "Cannot unify a record type with a non-record type."
 
 linkVarToDim ::
   MonadUnify m =>