SSE4.1 and AVX2 highbd_blend_a64_d16_mask

Add SSE4.1 and AVX2 SIMD optimised implementations of
highbd_blend_a64_d16_mask.

Unit tests speed-ups is approximately 2x to 2.2x for SSE4.1 and 2.8x to
3.6x for AVX2.

On 15 frames of crowd_run_360 at 10/12bd cpu=0 speed-up is about 3-4%.
At cpu=1 speed-up is about 3%.

I also corrected a comment at the top of blend_a64_mask.c which had not
been updated to reflect the rename of d32 blend functions to d16.

Change-Id: I789932a792c78e1bc64eabbf7143ecb0e67532da
diff --git a/aom_dsp/aom_dsp_rtcd_defs.pl b/aom_dsp/aom_dsp_rtcd_defs.pl
index 51e0621..704e358 100755
--- a/aom_dsp/aom_dsp_rtcd_defs.pl
+++ b/aom_dsp/aom_dsp_rtcd_defs.pl
@@ -544,7 +544,6 @@
 #
 add_proto qw/void aom_lowbd_blend_a64_d16_mask/, "uint8_t *dst, uint32_t dst_stride, const CONV_BUF_TYPE *src0, uint32_t src0_stride, const CONV_BUF_TYPE *src1, uint32_t src1_stride, const uint8_t *mask, uint32_t mask_stride, int w, int h, int subx, int suby, ConvolveParams *conv_params";
 specialize qw/aom_lowbd_blend_a64_d16_mask sse4_1 avx2 neon/;
-add_proto qw/void aom_highbd_blend_a64_d16_mask/, "uint8_t *dst, uint32_t dst_stride, const CONV_BUF_TYPE *src0, uint32_t src0_stride, const CONV_BUF_TYPE *src1, uint32_t src1_stride, const uint8_t *mask, uint32_t mask_stride, int w, int h, int subx, int suby, ConvolveParams *conv_params, const int bd";
 add_proto qw/void aom_blend_a64_mask/, "uint8_t *dst, uint32_t dst_stride, const uint8_t *src0, uint32_t src0_stride, const uint8_t *src1, uint32_t src1_stride, const uint8_t *mask, uint32_t mask_stride, int w, int h, int subx, int suby";
 add_proto qw/void aom_blend_a64_hmask/, "uint8_t *dst, uint32_t dst_stride, const uint8_t *src0, uint32_t src0_stride, const uint8_t *src1, uint32_t src1_stride, const uint8_t *mask, int w, int h";
 add_proto qw/void aom_blend_a64_vmask/, "uint8_t *dst, uint32_t dst_stride, const uint8_t *src0, uint32_t src0_stride, const uint8_t *src1, uint32_t src1_stride, const uint8_t *mask, int w, int h";
@@ -555,9 +554,11 @@
 add_proto qw/void aom_highbd_blend_a64_mask/, "uint8_t *dst, uint32_t dst_stride, const uint8_t *src0, uint32_t src0_stride, const uint8_t *src1, uint32_t src1_stride, const uint8_t *mask, uint32_t mask_stride, int w, int h, int subx, int suby, int bd";
 add_proto qw/void aom_highbd_blend_a64_hmask/, "uint8_t *dst, uint32_t dst_stride, const uint8_t *src0, uint32_t src0_stride, const uint8_t *src1, uint32_t src1_stride, const uint8_t *mask, int w, int h, int bd";
 add_proto qw/void aom_highbd_blend_a64_vmask/, "uint8_t *dst, uint32_t dst_stride, const uint8_t *src0, uint32_t src0_stride, const uint8_t *src1, uint32_t src1_stride, const uint8_t *mask, int w, int h, int bd";
+add_proto qw/void aom_highbd_blend_a64_d16_mask/, "uint8_t *dst, uint32_t dst_stride, const CONV_BUF_TYPE *src0, uint32_t src0_stride, const CONV_BUF_TYPE *src1, uint32_t src1_stride, const uint8_t *mask, uint32_t mask_stride, int w, int h, int subx, int suby, ConvolveParams *conv_params, const int bd";
 specialize "aom_highbd_blend_a64_mask", qw/sse4_1/;
 specialize "aom_highbd_blend_a64_hmask", qw/sse4_1/;
 specialize "aom_highbd_blend_a64_vmask", qw/sse4_1/;
+specialize "aom_highbd_blend_a64_d16_mask", qw/sse4_1 avx2/;
 
 if (aom_config("CONFIG_AV1_ENCODER") eq "yes") {
   #
diff --git a/aom_dsp/blend_a64_mask.c b/aom_dsp/blend_a64_mask.c
index 992cc5c..79956c3 100644
--- a/aom_dsp/blend_a64_mask.c
+++ b/aom_dsp/blend_a64_mask.c
@@ -22,7 +22,7 @@
 // as described for AOM_BLEND_A64 in aom_dsp/blend.h. src0 or src1 can
 // be the same as dst, or dst can be different from both sources.
 
-// NOTE(david.barker): The input and output of aom_blend_a64_d32_mask_c() are
+// NOTE(david.barker): The input and output of aom_blend_a64_d16_mask_c() are
 // in a higher intermediate precision, and will later be rounded down to pixel
 // precision.
 // Thus, in order to avoid double-rounding, we want to use normal right shifts
@@ -30,7 +30,7 @@
 // This works because of the identity:
 // ROUND_POWER_OF_TWO(x >> y, z) == ROUND_POWER_OF_TWO(x, y+z)
 //
-// In contrast, the output of the non-d32 functions will not be further rounded,
+// In contrast, the output of the non-d16 functions will not be further rounded,
 // so we *should* use ROUND_POWER_OF_TWO there.
 
 void aom_lowbd_blend_a64_d16_mask_c(
diff --git a/aom_dsp/x86/blend_a64_mask_avx2.c b/aom_dsp/x86/blend_a64_mask_avx2.c
index 67fb4d3..057f615 100644
--- a/aom_dsp/x86/blend_a64_mask_avx2.c
+++ b/aom_dsp/x86/blend_a64_mask_avx2.c
@@ -898,3 +898,475 @@
     }
   }
 }
+
+//////////////////////////////////////////////////////////////////////////////
+// aom_highbd_blend_a64_d16_mask_avx2()
+//////////////////////////////////////////////////////////////////////////////
+
+static INLINE void highbd_blend_a64_d16_mask_w4_avx2(
+    uint16_t *dst, int dst_stride, const CONV_BUF_TYPE *src0, int src0_stride,
+    const CONV_BUF_TYPE *src1, int src1_stride, const __m256i *mask0,
+    const __m256i *round_offset, int shift, const __m256i *clip_low,
+    const __m256i *clip_high, const __m256i *mask_max) {
+  // Load 4x u16 pixels from each of 4 rows from each source
+  const __m256i s0 = _mm256_set_epi64x(*(uint64_t *)(src0 + 3 * src0_stride),
+                                       *(uint64_t *)(src0 + 2 * src0_stride),
+                                       *(uint64_t *)(src0 + 1 * src0_stride),
+                                       *(uint64_t *)(src0 + 0 * src0_stride));
+  const __m256i s1 = _mm256_set_epi64x(*(uint64_t *)(src1 + 3 * src1_stride),
+                                       *(uint64_t *)(src1 + 2 * src1_stride),
+                                       *(uint64_t *)(src1 + 1 * src1_stride),
+                                       *(uint64_t *)(src1 + 0 * src1_stride));
+  // Generate the inverse mask
+  const __m256i mask1 = _mm256_sub_epi16(*mask_max, *mask0);
+
+  // Multiply each mask by the respective source
+  const __m256i mul0_highs = _mm256_mulhi_epu16(*mask0, s0);
+  const __m256i mul0_lows = _mm256_mullo_epi16(*mask0, s0);
+  const __m256i mul0h = _mm256_unpackhi_epi16(mul0_lows, mul0_highs);
+  const __m256i mul0l = _mm256_unpacklo_epi16(mul0_lows, mul0_highs);
+  // Note that AVX2 unpack orders 64-bit words as [3 1] [2 0] to keep within
+  // lanes Later, packs does the same again which cancels this out with no need
+  // for a permute.  The intermediate values being reordered makes no difference
+
+  const __m256i mul1_highs = _mm256_mulhi_epu16(mask1, s1);
+  const __m256i mul1_lows = _mm256_mullo_epi16(mask1, s1);
+  const __m256i mul1h = _mm256_unpackhi_epi16(mul1_lows, mul1_highs);
+  const __m256i mul1l = _mm256_unpacklo_epi16(mul1_lows, mul1_highs);
+
+  const __m256i sumh = _mm256_add_epi32(mul0h, mul1h);
+  const __m256i suml = _mm256_add_epi32(mul0l, mul1l);
+
+  const __m256i roundh =
+      _mm256_srai_epi32(_mm256_sub_epi32(sumh, *round_offset), shift);
+  const __m256i roundl =
+      _mm256_srai_epi32(_mm256_sub_epi32(suml, *round_offset), shift);
+
+  const __m256i pack = _mm256_packs_epi32(roundl, roundh);
+  const __m256i clip =
+      _mm256_min_epi16(_mm256_max_epi16(pack, *clip_low), *clip_high);
+
+  // _mm256_extract_epi64 doesn't exist on x86, so do it the old-fashioned way:
+  const __m128i cliph = _mm256_extracti128_si256(clip, 1);
+  xx_storel_64(dst + 3 * dst_stride, _mm_srli_si128(cliph, 8));
+  xx_storel_64(dst + 2 * dst_stride, cliph);
+  const __m128i clipl = _mm256_castsi256_si128(clip);
+  xx_storel_64(dst + 1 * dst_stride, _mm_srli_si128(clipl, 8));
+  xx_storel_64(dst + 0 * dst_stride, clipl);
+}
+
+static INLINE void highbd_blend_a64_d16_mask_subw0_subh0_w4_avx2(
+    uint16_t *dst, uint32_t dst_stride, const CONV_BUF_TYPE *src0,
+    uint32_t src0_stride, const CONV_BUF_TYPE *src1, uint32_t src1_stride,
+    const uint8_t *mask, uint32_t mask_stride, int h,
+    const __m256i *round_offset, int shift, const __m256i *clip_low,
+    const __m256i *clip_high, const __m256i *mask_max) {
+  do {
+    // Load 8x u8 pixels from each of 4 rows of the mask, pad each to u16
+    const __m128i mask08 = _mm_set_epi32(*(uint32_t *)(mask + 3 * mask_stride),
+                                         *(uint32_t *)(mask + 2 * mask_stride),
+                                         *(uint32_t *)(mask + 1 * mask_stride),
+                                         *(uint32_t *)(mask + 0 * mask_stride));
+    const __m256i mask0 = _mm256_cvtepu8_epi16(mask08);
+
+    highbd_blend_a64_d16_mask_w4_avx2(dst, dst_stride, src0, src0_stride, src1,
+                                      src1_stride, &mask0, round_offset, shift,
+                                      clip_low, clip_high, mask_max);
+
+    dst += dst_stride * 4;
+    src0 += src0_stride * 4;
+    src1 += src1_stride * 4;
+    mask += mask_stride * 4;
+  } while (h -= 4);
+}
+
+static INLINE void highbd_blend_a64_d16_mask_subw1_subh1_w4_avx2(
+    uint16_t *dst, uint32_t dst_stride, const CONV_BUF_TYPE *src0,
+    uint32_t src0_stride, const CONV_BUF_TYPE *src1, uint32_t src1_stride,
+    const uint8_t *mask, uint32_t mask_stride, int h,
+    const __m256i *round_offset, int shift, const __m256i *clip_low,
+    const __m256i *clip_high, const __m256i *mask_max) {
+  const __m256i one_b = _mm256_set1_epi8(1);
+  const __m256i two_w = _mm256_set1_epi16(2);
+  do {
+    // Load 8 pixels from each of 8 rows of mask,
+    // (saturating) add together rows then use madd to add adjacent pixels
+    // Finally, divide each value by 4 (with rounding)
+    const __m256i m0246 =
+        _mm256_set_epi64x(*(uint64_t *)(mask + 6 * mask_stride),
+                          *(uint64_t *)(mask + 4 * mask_stride),
+                          *(uint64_t *)(mask + 2 * mask_stride),
+                          *(uint64_t *)(mask + 0 * mask_stride));
+    const __m256i m1357 =
+        _mm256_set_epi64x(*(uint64_t *)(mask + 7 * mask_stride),
+                          *(uint64_t *)(mask + 5 * mask_stride),
+                          *(uint64_t *)(mask + 3 * mask_stride),
+                          *(uint64_t *)(mask + 1 * mask_stride));
+    const __m256i addrows = _mm256_adds_epu8(m0246, m1357);
+    const __m256i adjacent = _mm256_maddubs_epi16(addrows, one_b);
+    const __m256i mask0 =
+        _mm256_srli_epi16(_mm256_add_epi16(adjacent, two_w), 2);
+
+    highbd_blend_a64_d16_mask_w4_avx2(dst, dst_stride, src0, src0_stride, src1,
+                                      src1_stride, &mask0, round_offset, shift,
+                                      clip_low, clip_high, mask_max);
+
+    dst += dst_stride * 4;
+    src0 += src0_stride * 4;
+    src1 += src1_stride * 4;
+    mask += mask_stride * 8;
+  } while (h -= 4);
+}
+
+static INLINE void highbd_blend_a64_d16_mask_w8_avx2(
+    uint16_t *dst, int dst_stride, const CONV_BUF_TYPE *src0, int src0_stride,
+    const CONV_BUF_TYPE *src1, int src1_stride, const __m256i *mask0a,
+    const __m256i *mask0b, const __m256i *round_offset, int shift,
+    const __m256i *clip_low, const __m256i *clip_high,
+    const __m256i *mask_max) {
+  // Load 8x u16 pixels from each of 4 rows from each source
+  const __m256i s0a =
+      yy_loadu2_128(src0 + 0 * src0_stride, src0 + 1 * src0_stride);
+  const __m256i s0b =
+      yy_loadu2_128(src0 + 2 * src0_stride, src0 + 3 * src0_stride);
+  const __m256i s1a =
+      yy_loadu2_128(src1 + 0 * src1_stride, src1 + 1 * src1_stride);
+  const __m256i s1b =
+      yy_loadu2_128(src1 + 2 * src1_stride, src1 + 3 * src1_stride);
+
+  // Generate inverse masks
+  const __m256i mask1a = _mm256_sub_epi16(*mask_max, *mask0a);
+  const __m256i mask1b = _mm256_sub_epi16(*mask_max, *mask0b);
+
+  // Multiply sources by respective masks
+  const __m256i mul0a_highs = _mm256_mulhi_epu16(*mask0a, s0a);
+  const __m256i mul0a_lows = _mm256_mullo_epi16(*mask0a, s0a);
+  const __m256i mul0ah = _mm256_unpackhi_epi16(mul0a_lows, mul0a_highs);
+  const __m256i mul0al = _mm256_unpacklo_epi16(mul0a_lows, mul0a_highs);
+  // Note that AVX2 unpack orders 64-bit words as [3 1] [2 0] to keep within
+  // lanes Later, packs does the same again which cancels this out with no need
+  // for a permute.  The intermediate values being reordered makes no difference
+
+  const __m256i mul1a_highs = _mm256_mulhi_epu16(mask1a, s1a);
+  const __m256i mul1a_lows = _mm256_mullo_epi16(mask1a, s1a);
+  const __m256i mul1ah = _mm256_unpackhi_epi16(mul1a_lows, mul1a_highs);
+  const __m256i mul1al = _mm256_unpacklo_epi16(mul1a_lows, mul1a_highs);
+
+  const __m256i sumah = _mm256_add_epi32(mul0ah, mul1ah);
+  const __m256i sumal = _mm256_add_epi32(mul0al, mul1al);
+
+  const __m256i mul0b_highs = _mm256_mulhi_epu16(*mask0b, s0b);
+  const __m256i mul0b_lows = _mm256_mullo_epi16(*mask0b, s0b);
+  const __m256i mul0bh = _mm256_unpackhi_epi16(mul0b_lows, mul0b_highs);
+  const __m256i mul0bl = _mm256_unpacklo_epi16(mul0b_lows, mul0b_highs);
+
+  const __m256i mul1b_highs = _mm256_mulhi_epu16(mask1b, s1b);
+  const __m256i mul1b_lows = _mm256_mullo_epi16(mask1b, s1b);
+  const __m256i mul1bh = _mm256_unpackhi_epi16(mul1b_lows, mul1b_highs);
+  const __m256i mul1bl = _mm256_unpacklo_epi16(mul1b_lows, mul1b_highs);
+
+  const __m256i sumbh = _mm256_add_epi32(mul0bh, mul1bh);
+  const __m256i sumbl = _mm256_add_epi32(mul0bl, mul1bl);
+
+  // Divide down each result, with rounding
+  const __m256i roundah =
+      _mm256_srai_epi32(_mm256_sub_epi32(sumah, *round_offset), shift);
+  const __m256i roundal =
+      _mm256_srai_epi32(_mm256_sub_epi32(sumal, *round_offset), shift);
+  const __m256i roundbh =
+      _mm256_srai_epi32(_mm256_sub_epi32(sumbh, *round_offset), shift);
+  const __m256i roundbl =
+      _mm256_srai_epi32(_mm256_sub_epi32(sumbl, *round_offset), shift);
+
+  // Pack each i32 down to an i16 with saturation, then clip to valid range
+  const __m256i packa = _mm256_packs_epi32(roundal, roundah);
+  const __m256i clipa =
+      _mm256_min_epi16(_mm256_max_epi16(packa, *clip_low), *clip_high);
+  const __m256i packb = _mm256_packs_epi32(roundbl, roundbh);
+  const __m256i clipb =
+      _mm256_min_epi16(_mm256_max_epi16(packb, *clip_low), *clip_high);
+
+  // Store 8x u16 pixels to each of 4 rows in the destination
+  yy_storeu2_128(dst + 0 * dst_stride, dst + 1 * dst_stride, clipa);
+  yy_storeu2_128(dst + 2 * dst_stride, dst + 3 * dst_stride, clipb);
+}
+
+static INLINE void highbd_blend_a64_d16_mask_subw0_subh0_w8_avx2(
+    uint16_t *dst, int dst_stride, const CONV_BUF_TYPE *src0, int src0_stride,
+    const CONV_BUF_TYPE *src1, int src1_stride, const uint8_t *mask,
+    int mask_stride, int h, const __m256i *round_offset, int shift,
+    const __m256i *clip_low, const __m256i *clip_high,
+    const __m256i *mask_max) {
+  do {
+    // Load 8x u8 pixels from each of 4 rows in the mask
+    const __m128i mask0a8 =
+        _mm_set_epi64x(*(uint64_t *)mask, *(uint64_t *)(mask + mask_stride));
+    const __m128i mask0b8 =
+        _mm_set_epi64x(*(uint64_t *)(mask + 2 * mask_stride),
+                       *(uint64_t *)(mask + 3 * mask_stride));
+    const __m256i mask0a = _mm256_cvtepu8_epi16(mask0a8);
+    const __m256i mask0b = _mm256_cvtepu8_epi16(mask0b8);
+
+    highbd_blend_a64_d16_mask_w8_avx2(
+        dst, dst_stride, src0, src0_stride, src1, src1_stride, &mask0a, &mask0b,
+        round_offset, shift, clip_low, clip_high, mask_max);
+
+    dst += dst_stride * 4;
+    src0 += src0_stride * 4;
+    src1 += src1_stride * 4;
+    mask += mask_stride * 4;
+  } while (h -= 4);
+}
+
+static INLINE void highbd_blend_a64_d16_mask_subw1_subh1_w8_avx2(
+    uint16_t *dst, int dst_stride, const CONV_BUF_TYPE *src0, int src0_stride,
+    const CONV_BUF_TYPE *src1, int src1_stride, const uint8_t *mask,
+    int mask_stride, int h, const __m256i *round_offset, int shift,
+    const __m256i *clip_low, const __m256i *clip_high,
+    const __m256i *mask_max) {
+  const __m256i one_b = _mm256_set1_epi8(1);
+  const __m256i two_w = _mm256_set1_epi16(2);
+  do {
+    // Load 16x u8 pixels from each of 8 rows in the mask,
+    // (saturating) add together rows then use madd to add adjacent pixels
+    // Finally, divide each value by 4 (with rounding)
+    const __m256i m02 =
+        yy_loadu2_128(mask + 0 * mask_stride, mask + 2 * mask_stride);
+    const __m256i m13 =
+        yy_loadu2_128(mask + 1 * mask_stride, mask + 3 * mask_stride);
+    const __m256i m0123 =
+        _mm256_maddubs_epi16(_mm256_adds_epu8(m02, m13), one_b);
+    const __m256i mask_0a =
+        _mm256_srli_epi16(_mm256_add_epi16(m0123, two_w), 2);
+    const __m256i m46 =
+        yy_loadu2_128(mask + 4 * mask_stride, mask + 6 * mask_stride);
+    const __m256i m57 =
+        yy_loadu2_128(mask + 5 * mask_stride, mask + 7 * mask_stride);
+    const __m256i m4567 =
+        _mm256_maddubs_epi16(_mm256_adds_epu8(m46, m57), one_b);
+    const __m256i mask_0b =
+        _mm256_srli_epi16(_mm256_add_epi16(m4567, two_w), 2);
+
+    highbd_blend_a64_d16_mask_w8_avx2(
+        dst, dst_stride, src0, src0_stride, src1, src1_stride, &mask_0a,
+        &mask_0b, round_offset, shift, clip_low, clip_high, mask_max);
+
+    dst += dst_stride * 4;
+    src0 += src0_stride * 4;
+    src1 += src1_stride * 4;
+    mask += mask_stride * 8;
+  } while (h -= 4);
+}
+
+static INLINE void highbd_blend_a64_d16_mask_w16_avx2(
+    uint16_t *dst, int dst_stride, const CONV_BUF_TYPE *src0, int src0_stride,
+    const CONV_BUF_TYPE *src1, int src1_stride, const __m256i *mask0a,
+    const __m256i *mask0b, const __m256i *round_offset, int shift,
+    const __m256i *clip_low, const __m256i *clip_high,
+    const __m256i *mask_max) {
+  // Load 16x pixels from each of 2 rows from each source
+  const __m256i s0a = yy_loadu_256(src0);
+  const __m256i s0b = yy_loadu_256(src0 + src0_stride);
+  const __m256i s1a = yy_loadu_256(src1);
+  const __m256i s1b = yy_loadu_256(src1 + src1_stride);
+
+  // Calculate inverse masks
+  const __m256i mask1a = _mm256_sub_epi16(*mask_max, *mask0a);
+  const __m256i mask1b = _mm256_sub_epi16(*mask_max, *mask0b);
+
+  // Multiply each source by appropriate mask
+  const __m256i mul0a_highs = _mm256_mulhi_epu16(*mask0a, s0a);
+  const __m256i mul0a_lows = _mm256_mullo_epi16(*mask0a, s0a);
+  const __m256i mul0ah = _mm256_unpackhi_epi16(mul0a_lows, mul0a_highs);
+  const __m256i mul0al = _mm256_unpacklo_epi16(mul0a_lows, mul0a_highs);
+  // Note that AVX2 unpack orders 64-bit words as [3 1] [2 0] to keep within
+  // lanes Later, packs does the same again which cancels this out with no need
+  // for a permute.  The intermediate values being reordered makes no difference
+
+  const __m256i mul1a_highs = _mm256_mulhi_epu16(mask1a, s1a);
+  const __m256i mul1a_lows = _mm256_mullo_epi16(mask1a, s1a);
+  const __m256i mul1ah = _mm256_unpackhi_epi16(mul1a_lows, mul1a_highs);
+  const __m256i mul1al = _mm256_unpacklo_epi16(mul1a_lows, mul1a_highs);
+
+  const __m256i mulah = _mm256_add_epi32(mul0ah, mul1ah);
+  const __m256i mulal = _mm256_add_epi32(mul0al, mul1al);
+
+  const __m256i mul0b_highs = _mm256_mulhi_epu16(*mask0b, s0b);
+  const __m256i mul0b_lows = _mm256_mullo_epi16(*mask0b, s0b);
+  const __m256i mul0bh = _mm256_unpackhi_epi16(mul0b_lows, mul0b_highs);
+  const __m256i mul0bl = _mm256_unpacklo_epi16(mul0b_lows, mul0b_highs);
+
+  const __m256i mul1b_highs = _mm256_mulhi_epu16(mask1b, s1b);
+  const __m256i mul1b_lows = _mm256_mullo_epi16(mask1b, s1b);
+  const __m256i mul1bh = _mm256_unpackhi_epi16(mul1b_lows, mul1b_highs);
+  const __m256i mul1bl = _mm256_unpacklo_epi16(mul1b_lows, mul1b_highs);
+
+  const __m256i mulbh = _mm256_add_epi32(mul0bh, mul1bh);
+  const __m256i mulbl = _mm256_add_epi32(mul0bl, mul1bl);
+
+  const __m256i resah =
+      _mm256_srai_epi32(_mm256_sub_epi32(mulah, *round_offset), shift);
+  const __m256i resal =
+      _mm256_srai_epi32(_mm256_sub_epi32(mulal, *round_offset), shift);
+  const __m256i resbh =
+      _mm256_srai_epi32(_mm256_sub_epi32(mulbh, *round_offset), shift);
+  const __m256i resbl =
+      _mm256_srai_epi32(_mm256_sub_epi32(mulbl, *round_offset), shift);
+
+  // Signed saturating pack from i32 to i16:
+  const __m256i packa = _mm256_packs_epi32(resal, resah);
+  const __m256i packb = _mm256_packs_epi32(resbl, resbh);
+
+  // Clip the values to the valid range
+  const __m256i clipa =
+      _mm256_min_epi16(_mm256_max_epi16(packa, *clip_low), *clip_high);
+  const __m256i clipb =
+      _mm256_min_epi16(_mm256_max_epi16(packb, *clip_low), *clip_high);
+
+  // Store 16 pixels
+  yy_storeu_256(dst, clipa);
+  yy_storeu_256(dst + dst_stride, clipb);
+}
+
+static INLINE void highbd_blend_a64_d16_mask_subw0_subh0_w16_avx2(
+    uint16_t *dst, int dst_stride, const CONV_BUF_TYPE *src0, int src0_stride,
+    const CONV_BUF_TYPE *src1, int src1_stride, const uint8_t *mask,
+    int mask_stride, int h, int w, const __m256i *round_offset, int shift,
+    const __m256i *clip_low, const __m256i *clip_high,
+    const __m256i *mask_max) {
+  for (int i = 0; i < h; i += 2) {
+    for (int j = 0; j < w; j += 16) {
+      // Load 16x u8 alpha-mask values from each of two rows and pad to u16
+      const __m128i masks_a8 = xx_loadu_128(mask + j);
+      const __m128i masks_b8 = xx_loadu_128(mask + mask_stride + j);
+      const __m256i mask0a = _mm256_cvtepu8_epi16(masks_a8);
+      const __m256i mask0b = _mm256_cvtepu8_epi16(masks_b8);
+
+      highbd_blend_a64_d16_mask_w16_avx2(
+          dst + j, dst_stride, src0 + j, src0_stride, src1 + j, src1_stride,
+          &mask0a, &mask0b, round_offset, shift, clip_low, clip_high, mask_max);
+    }
+    dst += dst_stride * 2;
+    src0 += src0_stride * 2;
+    src1 += src1_stride * 2;
+    mask += mask_stride * 2;
+  }
+}
+
+static INLINE void highbd_blend_a64_d16_mask_subw1_subh1_w16_avx2(
+    uint16_t *dst, int dst_stride, const CONV_BUF_TYPE *src0, int src0_stride,
+    const CONV_BUF_TYPE *src1, int src1_stride, const uint8_t *mask,
+    int mask_stride, int h, int w, const __m256i *round_offset, int shift,
+    const __m256i *clip_low, const __m256i *clip_high,
+    const __m256i *mask_max) {
+  const __m256i one_b = _mm256_set1_epi8(1);
+  const __m256i two_w = _mm256_set1_epi16(2);
+  for (int i = 0; i < h; i += 2) {
+    for (int j = 0; j < w; j += 16) {
+      // Load 32x u8 alpha-mask values from each of four rows
+      // (saturating) add pairs of rows, then use madd to add adjacent values
+      // Finally, divide down each result with rounding
+      const __m256i m0 = yy_loadu_256(mask + 0 * mask_stride + 2 * j);
+      const __m256i m1 = yy_loadu_256(mask + 1 * mask_stride + 2 * j);
+      const __m256i m2 = yy_loadu_256(mask + 2 * mask_stride + 2 * j);
+      const __m256i m3 = yy_loadu_256(mask + 3 * mask_stride + 2 * j);
+
+      const __m256i m01_8 = _mm256_adds_epu8(m0, m1);
+      const __m256i m23_8 = _mm256_adds_epu8(m2, m3);
+
+      const __m256i m01 = _mm256_maddubs_epi16(m01_8, one_b);
+      const __m256i m23 = _mm256_maddubs_epi16(m23_8, one_b);
+
+      const __m256i mask0a = _mm256_srli_epi16(_mm256_add_epi16(m01, two_w), 2);
+      const __m256i mask0b = _mm256_srli_epi16(_mm256_add_epi16(m23, two_w), 2);
+
+      highbd_blend_a64_d16_mask_w16_avx2(
+          dst + j, dst_stride, src0 + j, src0_stride, src1 + j, src1_stride,
+          &mask0a, &mask0b, round_offset, shift, clip_low, clip_high, mask_max);
+    }
+    dst += dst_stride * 2;
+    src0 += src0_stride * 2;
+    src1 += src1_stride * 2;
+    mask += mask_stride * 4;
+  }
+}
+
+void aom_highbd_blend_a64_d16_mask_avx2(
+    uint8_t *dst8, uint32_t dst_stride, const CONV_BUF_TYPE *src0,
+    uint32_t src0_stride, const CONV_BUF_TYPE *src1, uint32_t src1_stride,
+    const uint8_t *mask, uint32_t mask_stride, int w, int h, int subw, int subh,
+    ConvolveParams *conv_params, const int bd) {
+  uint16_t *dst = CONVERT_TO_SHORTPTR(dst8);
+  const int round_bits =
+      2 * FILTER_BITS - conv_params->round_0 - conv_params->round_1;
+  const int32_t round_offset =
+      ((1 << (round_bits + bd)) + (1 << (round_bits + bd - 1)) -
+       (1 << (round_bits - 1)))
+      << AOM_BLEND_A64_ROUND_BITS;
+  const __m256i v_round_offset = _mm256_set1_epi32(round_offset);
+  const int shift = round_bits + AOM_BLEND_A64_ROUND_BITS;
+
+  const __m256i clip_low = _mm256_set1_epi16(0);
+  const __m256i clip_high = _mm256_set1_epi16((1 << bd) - 1);
+  const __m256i mask_max = _mm256_set1_epi16(AOM_BLEND_A64_MAX_ALPHA);
+
+  assert(IMPLIES((void *)src0 == dst, src0_stride == dst_stride));
+  assert(IMPLIES((void *)src1 == dst, src1_stride == dst_stride));
+
+  assert(h >= 4);
+  assert(w >= 4);
+  assert(IS_POWER_OF_TWO(h));
+  assert(IS_POWER_OF_TWO(w));
+
+  if (subw == 0 && subh == 0) {
+    switch (w) {
+      case 4:
+        highbd_blend_a64_d16_mask_subw0_subh0_w4_avx2(
+            dst, dst_stride, src0, src0_stride, src1, src1_stride, mask,
+            mask_stride, h, &v_round_offset, shift, &clip_low, &clip_high,
+            &mask_max);
+        break;
+      case 8:
+        highbd_blend_a64_d16_mask_subw0_subh0_w8_avx2(
+            dst, dst_stride, src0, src0_stride, src1, src1_stride, mask,
+            mask_stride, h, &v_round_offset, shift, &clip_low, &clip_high,
+            &mask_max);
+        break;
+      default:  // >= 16
+        highbd_blend_a64_d16_mask_subw0_subh0_w16_avx2(
+            dst, dst_stride, src0, src0_stride, src1, src1_stride, mask,
+            mask_stride, h, w, &v_round_offset, shift, &clip_low, &clip_high,
+            &mask_max);
+        break;
+    }
+
+  } else if (subw == 1 && subh == 1) {
+    switch (w) {
+      case 4:
+        highbd_blend_a64_d16_mask_subw1_subh1_w4_avx2(
+            dst, dst_stride, src0, src0_stride, src1, src1_stride, mask,
+            mask_stride, h, &v_round_offset, shift, &clip_low, &clip_high,
+            &mask_max);
+        break;
+      case 8:
+        highbd_blend_a64_d16_mask_subw1_subh1_w8_avx2(
+            dst, dst_stride, src0, src0_stride, src1, src1_stride, mask,
+            mask_stride, h, &v_round_offset, shift, &clip_low, &clip_high,
+            &mask_max);
+        break;
+      default:  // >= 16
+        highbd_blend_a64_d16_mask_subw1_subh1_w16_avx2(
+            dst, dst_stride, src0, src0_stride, src1, src1_stride, mask,
+            mask_stride, h, w, &v_round_offset, shift, &clip_low, &clip_high,
+            &mask_max);
+        break;
+    }
+  } else {
+    // Sub-sampling in only one axis doesn't seem to happen very much, so fall
+    // back to the vanilla C implementation instead of having all the optimised
+    // code for these.
+    aom_highbd_blend_a64_d16_mask_c(dst8, dst_stride, src0, src0_stride, src1,
+                                    src1_stride, mask, mask_stride, w, h, subw,
+                                    subh, conv_params, bd);
+  }
+}
diff --git a/aom_dsp/x86/blend_a64_mask_sse4.c b/aom_dsp/x86/blend_a64_mask_sse4.c
index 9d6b4c2..b7a2468 100644
--- a/aom_dsp/x86/blend_a64_mask_sse4.c
+++ b/aom_dsp/x86/blend_a64_mask_sse4.c
@@ -1107,3 +1107,452 @@
     }
   }
 }
+
+//////////////////////////////////////////////////////////////////////////////
+// aom_highbd_blend_a64_d16_mask_sse4_1()
+//////////////////////////////////////////////////////////////////////////////
+
+static INLINE void highbd_blend_a64_d16_mask_w4_sse4_1(
+    uint16_t *dst, int dst_stride, const CONV_BUF_TYPE *src0, int src0_stride,
+    const CONV_BUF_TYPE *src1, int src1_stride, const __m128i *mask0a,
+    const __m128i *mask0b, const __m128i *round_offset, int shift,
+    const __m128i *clip_low, const __m128i *clip_high,
+    const __m128i *mask_max) {
+  // Load 4 pixels from each of 4 rows from each source
+  const __m128i s0a =
+      _mm_set_epi64x(*(uint64_t *)src0, *(uint64_t *)(src0 + src0_stride));
+  const __m128i s0b = _mm_set_epi64x(*(uint64_t *)(src0 + 2 * src0_stride),
+                                     *(uint64_t *)(src0 + 3 * src0_stride));
+  const __m128i s1a =
+      _mm_set_epi64x(*(uint64_t *)(src1), *(uint64_t *)(src1 + src1_stride));
+  const __m128i s1b = _mm_set_epi64x(*(uint64_t *)(src1 + 2 * src1_stride),
+                                     *(uint64_t *)(src1 + 3 * src1_stride));
+
+  // Generate the inverse masks
+  const __m128i mask1a = _mm_sub_epi16(*mask_max, *mask0a);
+  const __m128i mask1b = _mm_sub_epi16(*mask_max, *mask0b);
+
+  // Multiply each mask by the respective source
+  const __m128i mul0a_highs = _mm_mulhi_epu16(*mask0a, s0a);
+  const __m128i mul0a_lows = _mm_mullo_epi16(*mask0a, s0a);
+  const __m128i mul0ah = _mm_unpackhi_epi16(mul0a_lows, mul0a_highs);
+  const __m128i mul0al = _mm_unpacklo_epi16(mul0a_lows, mul0a_highs);
+  const __m128i mul1a_highs = _mm_mulhi_epu16(mask1a, s1a);
+  const __m128i mul1a_lows = _mm_mullo_epi16(mask1a, s1a);
+  const __m128i mul1ah = _mm_unpackhi_epi16(mul1a_lows, mul1a_highs);
+  const __m128i mul1al = _mm_unpacklo_epi16(mul1a_lows, mul1a_highs);
+
+  const __m128i mul0b_highs = _mm_mulhi_epu16(*mask0b, s0b);
+  const __m128i mul0b_lows = _mm_mullo_epi16(*mask0b, s0b);
+  const __m128i mul0bh = _mm_unpackhi_epi16(mul0b_lows, mul0b_highs);
+  const __m128i mul0bl = _mm_unpacklo_epi16(mul0b_lows, mul0b_highs);
+  const __m128i mul1b_highs = _mm_mulhi_epu16(mask1b, s1b);
+  const __m128i mul1b_lows = _mm_mullo_epi16(mask1b, s1b);
+  const __m128i mul1bh = _mm_unpackhi_epi16(mul1b_lows, mul1b_highs);
+  const __m128i mul1bl = _mm_unpacklo_epi16(mul1b_lows, mul1b_highs);
+
+  const __m128i sumah = _mm_add_epi32(mul0ah, mul1ah);
+  const __m128i sumal = _mm_add_epi32(mul0al, mul1al);
+  const __m128i sumbh = _mm_add_epi32(mul0bh, mul1bh);
+  const __m128i sumbl = _mm_add_epi32(mul0bl, mul1bl);
+
+  const __m128i roundah =
+      _mm_srai_epi32(_mm_sub_epi32(sumah, *round_offset), shift);
+  const __m128i roundbh =
+      _mm_srai_epi32(_mm_sub_epi32(sumbh, *round_offset), shift);
+  const __m128i roundal =
+      _mm_srai_epi32(_mm_sub_epi32(sumal, *round_offset), shift);
+  const __m128i roundbl =
+      _mm_srai_epi32(_mm_sub_epi32(sumbl, *round_offset), shift);
+
+  const __m128i packa = _mm_packs_epi32(roundal, roundah);
+  const __m128i packb = _mm_packs_epi32(roundbl, roundbh);
+
+  const __m128i clipa =
+      _mm_min_epi16(_mm_max_epi16(packa, *clip_low), *clip_high);
+  const __m128i clipb =
+      _mm_min_epi16(_mm_max_epi16(packb, *clip_low), *clip_high);
+
+  xx_storel_64(dst, _mm_srli_si128(clipa, 8));
+  xx_storel_64(dst + dst_stride, clipa);
+  xx_storel_64(dst + 2 * dst_stride, _mm_srli_si128(clipb, 8));
+  xx_storel_64(dst + 3 * dst_stride, clipb);
+}
+
+static INLINE void highbd_blend_a64_d16_mask_subw0_subh0_w4_sse4_1(
+    uint16_t *dst, uint32_t dst_stride, const CONV_BUF_TYPE *src0,
+    uint32_t src0_stride, const CONV_BUF_TYPE *src1, uint32_t src1_stride,
+    const uint8_t *mask, uint32_t mask_stride, int h,
+    const __m128i *round_offset, int shift, const __m128i *clip_low,
+    const __m128i *clip_high, const __m128i *mask_max) {
+  do {
+    const __m128i mask0a8 = _mm_set_epi32(0, 0, *(uint32_t *)mask,
+                                          *(uint32_t *)(mask + mask_stride));
+    const __m128i mask0b8 =
+        _mm_set_epi32(0, 0, *(uint32_t *)(mask + 2 * mask_stride),
+                      *(uint32_t *)(mask + 3 * mask_stride));
+    const __m128i mask0a = _mm_cvtepu8_epi16(mask0a8);
+    const __m128i mask0b = _mm_cvtepu8_epi16(mask0b8);
+
+    highbd_blend_a64_d16_mask_w4_sse4_1(
+        dst, dst_stride, src0, src0_stride, src1, src1_stride, &mask0a, &mask0b,
+        round_offset, shift, clip_low, clip_high, mask_max);
+
+    dst += dst_stride * 4;
+    src0 += src0_stride * 4;
+    src1 += src1_stride * 4;
+    mask += mask_stride * 4;
+  } while (h -= 4);
+}
+
+static INLINE void highbd_blend_a64_d16_mask_subw1_subh1_w4_sse4_1(
+    uint16_t *dst, uint32_t dst_stride, const CONV_BUF_TYPE *src0,
+    uint32_t src0_stride, const CONV_BUF_TYPE *src1, uint32_t src1_stride,
+    const uint8_t *mask, uint32_t mask_stride, int h,
+    const __m128i *round_offset, int shift, const __m128i *clip_low,
+    const __m128i *clip_high, const __m128i *mask_max) {
+  const __m128i one_b = _mm_set1_epi8(1);
+  const __m128i two_w = _mm_set1_epi16(2);
+  do {
+    // Load 8 pixels from each of 8 rows of mask,
+    // (saturating) add together rows then use madd to add adjacent pixels
+    // Finally, divide each value by 4 (with rounding)
+    const __m128i m02 = _mm_set_epi64x(*(uint64_t *)(mask),
+                                       *(uint64_t *)(mask + 2 * mask_stride));
+    const __m128i m13 = _mm_set_epi64x(*(uint64_t *)(mask + mask_stride),
+                                       *(uint64_t *)(mask + 3 * mask_stride));
+    const __m128i m0123 = _mm_maddubs_epi16(_mm_adds_epu8(m02, m13), one_b);
+    const __m128i mask_0a = _mm_srli_epi16(_mm_add_epi16(m0123, two_w), 2);
+    const __m128i m46 = _mm_set_epi64x(*(uint64_t *)(mask + 4 * mask_stride),
+                                       *(uint64_t *)(mask + 6 * mask_stride));
+    const __m128i m57 = _mm_set_epi64x(*(uint64_t *)(mask + 5 * mask_stride),
+                                       *(uint64_t *)(mask + 7 * mask_stride));
+    const __m128i m4567 = _mm_maddubs_epi16(_mm_adds_epu8(m46, m57), one_b);
+    const __m128i mask_0b = _mm_srli_epi16(_mm_add_epi16(m4567, two_w), 2);
+
+    highbd_blend_a64_d16_mask_w4_sse4_1(
+        dst, dst_stride, src0, src0_stride, src1, src1_stride, &mask_0a,
+        &mask_0b, round_offset, shift, clip_low, clip_high, mask_max);
+
+    dst += dst_stride * 4;
+    src0 += src0_stride * 4;
+    src1 += src1_stride * 4;
+    mask += mask_stride * 8;
+  } while (h -= 4);
+}
+
+static INLINE void highbd_blend_a64_d16_mask_w8_sse4_1(
+    uint16_t *dst, int dst_stride, const CONV_BUF_TYPE *src0, int src0_stride,
+    const CONV_BUF_TYPE *src1, int src1_stride, const __m128i *mask0a,
+    const __m128i *mask0b, const __m128i *round_offset, int shift,
+    const __m128i *clip_low, const __m128i *clip_high,
+    const __m128i *max_mask) {
+  // Load 8x pixels from each of 2 rows from each source
+  const __m128i s0a = xx_loadu_128(src0);
+  const __m128i s0b = xx_loadu_128(src0 + src0_stride);
+  const __m128i s1a = xx_loadu_128(src1);
+  const __m128i s1b = xx_loadu_128(src1 + src1_stride);
+
+  // Generate inverse masks
+  const __m128i mask1a = _mm_sub_epi16(*max_mask, *mask0a);
+  const __m128i mask1b = _mm_sub_epi16(*max_mask, *mask0b);
+
+  // Multiply sources by respective masks
+  const __m128i mul0a_highs = _mm_mulhi_epu16(*mask0a, s0a);
+  const __m128i mul0a_lows = _mm_mullo_epi16(*mask0a, s0a);
+  const __m128i mul0ah = _mm_unpackhi_epi16(mul0a_lows, mul0a_highs);
+  const __m128i mul0al = _mm_unpacklo_epi16(mul0a_lows, mul0a_highs);
+
+  const __m128i mul1a_highs = _mm_mulhi_epu16(mask1a, s1a);
+  const __m128i mul1a_lows = _mm_mullo_epi16(mask1a, s1a);
+  const __m128i mul1ah = _mm_unpackhi_epi16(mul1a_lows, mul1a_highs);
+  const __m128i mul1al = _mm_unpacklo_epi16(mul1a_lows, mul1a_highs);
+
+  const __m128i sumah = _mm_add_epi32(mul0ah, mul1ah);
+  const __m128i sumal = _mm_add_epi32(mul0al, mul1al);
+
+  const __m128i mul0b_highs = _mm_mulhi_epu16(*mask0b, s0b);
+  const __m128i mul0b_lows = _mm_mullo_epi16(*mask0b, s0b);
+  const __m128i mul0bh = _mm_unpackhi_epi16(mul0b_lows, mul0b_highs);
+  const __m128i mul0bl = _mm_unpacklo_epi16(mul0b_lows, mul0b_highs);
+
+  const __m128i mul1b_highs = _mm_mulhi_epu16(mask1b, s1b);
+  const __m128i mul1b_lows = _mm_mullo_epi16(mask1b, s1b);
+  const __m128i mul1bh = _mm_unpackhi_epi16(mul1b_lows, mul1b_highs);
+  const __m128i mul1bl = _mm_unpacklo_epi16(mul1b_lows, mul1b_highs);
+
+  const __m128i sumbh = _mm_add_epi32(mul0bh, mul1bh);
+  const __m128i sumbl = _mm_add_epi32(mul0bl, mul1bl);
+
+  const __m128i roundah =
+      _mm_srai_epi32(_mm_sub_epi32(sumah, *round_offset), shift);
+  const __m128i roundal =
+      _mm_srai_epi32(_mm_sub_epi32(sumal, *round_offset), shift);
+  const __m128i roundbh =
+      _mm_srai_epi32(_mm_sub_epi32(sumbh, *round_offset), shift);
+  const __m128i roundbl =
+      _mm_srai_epi32(_mm_sub_epi32(sumbl, *round_offset), shift);
+
+  const __m128i packa = _mm_packs_epi32(roundal, roundah);
+  const __m128i clipa =
+      _mm_min_epi16(_mm_max_epi16(packa, *clip_low), *clip_high);
+  const __m128i packb = _mm_packs_epi32(roundbl, roundbh);
+  const __m128i clipb =
+      _mm_min_epi16(_mm_max_epi16(packb, *clip_low), *clip_high);
+
+  xx_storeu_128(dst, clipa);
+  xx_storeu_128(dst + dst_stride, clipb);
+}
+
+static INLINE void highbd_blend_a64_d16_mask_subw0_subh0_w8_sse4_1(
+    uint16_t *dst, uint32_t dst_stride, const CONV_BUF_TYPE *src0,
+    uint32_t src0_stride, const CONV_BUF_TYPE *src1, uint32_t src1_stride,
+    const uint8_t *mask, uint32_t mask_stride, int h,
+    const __m128i *round_offset, int shift, const __m128i *clip_low,
+    const __m128i *clip_high, const __m128i *max_mask) {
+  do {
+    const __m128i mask0a = _mm_cvtepu8_epi16(xx_loadl_64(mask));
+    const __m128i mask0b = _mm_cvtepu8_epi16(xx_loadl_64(mask + mask_stride));
+    highbd_blend_a64_d16_mask_w8_sse4_1(
+        dst, dst_stride, src0, src0_stride, src1, src1_stride, &mask0a, &mask0b,
+        round_offset, shift, clip_low, clip_high, max_mask);
+
+    dst += dst_stride * 2;
+    src0 += src0_stride * 2;
+    src1 += src1_stride * 2;
+    mask += mask_stride * 2;
+  } while (h -= 2);
+}
+
+static INLINE void highbd_blend_a64_d16_mask_subw1_subh1_w8_sse4_1(
+    uint16_t *dst, uint32_t dst_stride, const CONV_BUF_TYPE *src0,
+    uint32_t src0_stride, const CONV_BUF_TYPE *src1, uint32_t src1_stride,
+    const uint8_t *mask, uint32_t mask_stride, int h,
+    const __m128i *round_offset, int shift, const __m128i *clip_low,
+    const __m128i *clip_high, const __m128i *max_mask) {
+  const __m128i one_b = _mm_set1_epi8(1);
+  const __m128i two_w = _mm_set1_epi16(2);
+  do {
+    const __m128i mask_thisrowa = xx_loadu_128(mask);
+    const __m128i mask_nextrowa = xx_loadu_128(mask + mask_stride);
+    const __m128i mask_thisrowb = xx_loadu_128(mask + 2 * mask_stride);
+    const __m128i mask_nextrowb = xx_loadu_128(mask + 3 * mask_stride);
+    const __m128i mask_bothrowsa = _mm_adds_epu8(mask_thisrowa, mask_nextrowa);
+    const __m128i mask_bothrowsb = _mm_adds_epu8(mask_thisrowb, mask_nextrowb);
+    const __m128i mask_16a = _mm_maddubs_epi16(mask_bothrowsa, one_b);
+    const __m128i mask_16b = _mm_maddubs_epi16(mask_bothrowsb, one_b);
+    const __m128i mask_sa = _mm_srli_epi16(_mm_add_epi16(mask_16a, two_w), 2);
+    const __m128i mask_sb = _mm_srli_epi16(_mm_add_epi16(mask_16b, two_w), 2);
+
+    highbd_blend_a64_d16_mask_w8_sse4_1(
+        dst, dst_stride, src0, src0_stride, src1, src1_stride, &mask_sa,
+        &mask_sb, round_offset, shift, clip_low, clip_high, max_mask);
+
+    dst += dst_stride * 2;
+    src0 += src0_stride * 2;
+    src1 += src1_stride * 2;
+    mask += mask_stride * 4;
+  } while (h -= 2);
+}
+
+static INLINE void highbd_blend_a64_d16_mask_w16_sse4_1(
+    uint16_t *dst, const CONV_BUF_TYPE *src0, const CONV_BUF_TYPE *src1,
+    const __m128i *round_offset, int shift, const __m128i *mask0l,
+    const __m128i *mask0h, const __m128i *clip_low, const __m128i *clip_high,
+    const __m128i *mask_max) {
+  // Load 16x u16 pixels for this row from each src
+  const __m128i s0l = xx_loadu_128(src0);
+  const __m128i s0h = xx_loadu_128(src0 + 8);
+  const __m128i s1l = xx_loadu_128(src1);
+  const __m128i s1h = xx_loadu_128(src1 + 8);
+
+  // Calculate inverse masks
+  const __m128i mask1h = _mm_sub_epi16(*mask_max, *mask0h);
+  const __m128i mask1l = _mm_sub_epi16(*mask_max, *mask0l);
+
+  const __m128i mul0_highs = _mm_mulhi_epu16(*mask0h, s0h);
+  const __m128i mul0_lows = _mm_mullo_epi16(*mask0h, s0h);
+  const __m128i mul0h = _mm_unpackhi_epi16(mul0_lows, mul0_highs);
+  const __m128i mul0l = _mm_unpacklo_epi16(mul0_lows, mul0_highs);
+
+  const __m128i mul1_highs = _mm_mulhi_epu16(mask1h, s1h);
+  const __m128i mul1_lows = _mm_mullo_epi16(mask1h, s1h);
+  const __m128i mul1h = _mm_unpackhi_epi16(mul1_lows, mul1_highs);
+  const __m128i mul1l = _mm_unpacklo_epi16(mul1_lows, mul1_highs);
+
+  const __m128i mulhh = _mm_add_epi32(mul0h, mul1h);
+  const __m128i mulhl = _mm_add_epi32(mul0l, mul1l);
+
+  const __m128i mul2_highs = _mm_mulhi_epu16(*mask0l, s0l);
+  const __m128i mul2_lows = _mm_mullo_epi16(*mask0l, s0l);
+  const __m128i mul2h = _mm_unpackhi_epi16(mul2_lows, mul2_highs);
+  const __m128i mul2l = _mm_unpacklo_epi16(mul2_lows, mul2_highs);
+
+  const __m128i mul3_highs = _mm_mulhi_epu16(mask1l, s1l);
+  const __m128i mul3_lows = _mm_mullo_epi16(mask1l, s1l);
+  const __m128i mul3h = _mm_unpackhi_epi16(mul3_lows, mul3_highs);
+  const __m128i mul3l = _mm_unpacklo_epi16(mul3_lows, mul3_highs);
+
+  const __m128i mullh = _mm_add_epi32(mul2h, mul3h);
+  const __m128i mulll = _mm_add_epi32(mul2l, mul3l);
+
+  const __m128i reshh =
+      _mm_srai_epi32(_mm_sub_epi32(mulhh, *round_offset), shift);
+  const __m128i reshl =
+      _mm_srai_epi32(_mm_sub_epi32(mulhl, *round_offset), shift);
+  const __m128i reslh =
+      _mm_srai_epi32(_mm_sub_epi32(mullh, *round_offset), shift);
+  const __m128i resll =
+      _mm_srai_epi32(_mm_sub_epi32(mulll, *round_offset), shift);
+
+  // Signed saturating pack from i32 to i16:
+  const __m128i packh = _mm_packs_epi32(reshl, reshh);
+  const __m128i packl = _mm_packs_epi32(resll, reslh);
+
+  // Clip the values to the valid range
+  const __m128i cliph =
+      _mm_min_epi16(_mm_max_epi16(packh, *clip_low), *clip_high);
+  const __m128i clipl =
+      _mm_min_epi16(_mm_max_epi16(packl, *clip_low), *clip_high);
+
+  // Store 16 pixels
+  xx_storeu_128(dst, clipl);
+  xx_storeu_128(dst + 8, cliph);
+}
+
+static INLINE void highbd_blend_a64_d16_mask_subw0_subh0_w16_sse4_1(
+    uint16_t *dst, uint32_t dst_stride, const CONV_BUF_TYPE *src0,
+    uint32_t src0_stride, const CONV_BUF_TYPE *src1, uint32_t src1_stride,
+    const uint8_t *mask, uint32_t mask_stride, int h, int w,
+    const __m128i *round_offset, int shift, const __m128i *clip_low,
+    const __m128i *clip_high, const __m128i *mask_max) {
+  for (int i = 0; i < h; i++) {
+    for (int j = 0; j < w; j += 16) {
+      // Load 16x u8 alpha-mask values and pad to u16
+      const __m128i masks_u8 = xx_loadu_128(mask + j);
+      const __m128i mask0l = _mm_cvtepu8_epi16(masks_u8);
+      const __m128i mask0h = _mm_cvtepu8_epi16(_mm_srli_si128(masks_u8, 8));
+
+      highbd_blend_a64_d16_mask_w16_sse4_1(
+          dst + j, src0 + j, src1 + j, round_offset, shift, &mask0l, &mask0h,
+          clip_low, clip_high, mask_max);
+    }
+    dst += dst_stride;
+    src0 += src0_stride;
+    src1 += src1_stride;
+    mask += mask_stride;
+  }
+}
+
+static INLINE void highbd_blend_a64_d16_mask_subw1_subh1_w16_sse4_1(
+    uint16_t *dst, uint32_t dst_stride, const CONV_BUF_TYPE *src0,
+    uint32_t src0_stride, const CONV_BUF_TYPE *src1, uint32_t src1_stride,
+    const uint8_t *mask, uint32_t mask_stride, int h, int w,
+    const __m128i *round_offset, int shift, const __m128i *clip_low,
+    const __m128i *clip_high, const __m128i *mask_max) {
+  const __m128i one_b = _mm_set1_epi8(1);
+  const __m128i two_w = _mm_set1_epi16(2);
+  for (int i = 0; i < h; i++) {
+    for (int j = 0; j < w; j += 16) {
+      const __m128i m_i00 = xx_loadu_128(mask + 2 * j);
+      const __m128i m_i01 = xx_loadu_128(mask + 2 * j + 16);
+      const __m128i m_i10 = xx_loadu_128(mask + mask_stride + 2 * j);
+      const __m128i m_i11 = xx_loadu_128(mask + mask_stride + 2 * j + 16);
+
+      const __m128i m0_ac = _mm_adds_epu8(m_i00, m_i10);
+      const __m128i m1_ac = _mm_adds_epu8(m_i01, m_i11);
+      const __m128i m0_acbd = _mm_maddubs_epi16(m0_ac, one_b);
+      const __m128i m1_acbd = _mm_maddubs_epi16(m1_ac, one_b);
+      const __m128i mask_l = _mm_srli_epi16(_mm_add_epi16(m0_acbd, two_w), 2);
+      const __m128i mask_h = _mm_srli_epi16(_mm_add_epi16(m1_acbd, two_w), 2);
+
+      highbd_blend_a64_d16_mask_w16_sse4_1(
+          dst + j, src0 + j, src1 + j, round_offset, shift, &mask_l, &mask_h,
+          clip_low, clip_high, mask_max);
+    }
+    dst += dst_stride;
+    src0 += src0_stride;
+    src1 += src1_stride;
+    mask += mask_stride * 2;
+  }
+}
+
+void aom_highbd_blend_a64_d16_mask_sse4_1(
+    uint8_t *dst8, uint32_t dst_stride, const CONV_BUF_TYPE *src0,
+    uint32_t src0_stride, const CONV_BUF_TYPE *src1, uint32_t src1_stride,
+    const uint8_t *mask, uint32_t mask_stride, int w, int h, int subw, int subh,
+    ConvolveParams *conv_params, const int bd) {
+  uint16_t *dst = CONVERT_TO_SHORTPTR(dst8);
+  const int round_bits =
+      2 * FILTER_BITS - conv_params->round_0 - conv_params->round_1;
+  const int32_t round_offset =
+      ((1 << (round_bits + bd)) + (1 << (round_bits + bd - 1)) -
+       (1 << (round_bits - 1)))
+      << AOM_BLEND_A64_ROUND_BITS;
+  const __m128i v_round_offset = _mm_set1_epi32(round_offset);
+  const int shift = round_bits + AOM_BLEND_A64_ROUND_BITS;
+
+  const __m128i clip_low = _mm_set1_epi16(0);
+  const __m128i clip_high = _mm_set1_epi16((1 << bd) - 1);
+  const __m128i mask_max = _mm_set1_epi16(AOM_BLEND_A64_MAX_ALPHA);
+
+  assert(IMPLIES((void *)src0 == dst, src0_stride == dst_stride));
+  assert(IMPLIES((void *)src1 == dst, src1_stride == dst_stride));
+
+  assert(h >= 4);
+  assert(w >= 4);
+  assert(IS_POWER_OF_TWO(h));
+  assert(IS_POWER_OF_TWO(w));
+
+  if (subw == 0 && subh == 0) {
+    switch (w) {
+      case 4:
+        highbd_blend_a64_d16_mask_subw0_subh0_w4_sse4_1(
+            dst, dst_stride, src0, src0_stride, src1, src1_stride, mask,
+            mask_stride, h, &v_round_offset, shift, &clip_low, &clip_high,
+            &mask_max);
+        break;
+      case 8:
+        highbd_blend_a64_d16_mask_subw0_subh0_w8_sse4_1(
+            dst, dst_stride, src0, src0_stride, src1, src1_stride, mask,
+            mask_stride, h, &v_round_offset, shift, &clip_low, &clip_high,
+            &mask_max);
+        break;
+      default:  // >=16
+        highbd_blend_a64_d16_mask_subw0_subh0_w16_sse4_1(
+            dst, dst_stride, src0, src0_stride, src1, src1_stride, mask,
+            mask_stride, h, w, &v_round_offset, shift, &clip_low, &clip_high,
+            &mask_max);
+        break;
+    }
+
+  } else if (subw == 1 && subh == 1) {
+    switch (w) {
+      case 4:
+        highbd_blend_a64_d16_mask_subw1_subh1_w4_sse4_1(
+            dst, dst_stride, src0, src0_stride, src1, src1_stride, mask,
+            mask_stride, h, &v_round_offset, shift, &clip_low, &clip_high,
+            &mask_max);
+        break;
+      case 8:
+        highbd_blend_a64_d16_mask_subw1_subh1_w8_sse4_1(
+            dst, dst_stride, src0, src0_stride, src1, src1_stride, mask,
+            mask_stride, h, &v_round_offset, shift, &clip_low, &clip_high,
+            &mask_max);
+        break;
+      default:  // >=16
+        highbd_blend_a64_d16_mask_subw1_subh1_w16_sse4_1(
+            dst, dst_stride, src0, src0_stride, src1, src1_stride, mask,
+            mask_stride, h, w, &v_round_offset, shift, &clip_low, &clip_high,
+            &mask_max);
+        break;
+    }
+  } else {
+    // Sub-sampling in only one axis doesn't seem to happen very much, so fall
+    // back to the vanilla C implementation instead of having all the optimised
+    // code for these.
+    aom_highbd_blend_a64_d16_mask_c(dst8, dst_stride, src0, src0_stride, src1,
+                                    src1_stride, mask, mask_stride, w, h, subw,
+                                    subh, conv_params, bd);
+  }
+}
diff --git a/aom_dsp/x86/synonyms_avx2.h b/aom_dsp/x86/synonyms_avx2.h
index 3f69b12..4d6ee6a 100644
--- a/aom_dsp/x86/synonyms_avx2.h
+++ b/aom_dsp/x86/synonyms_avx2.h
@@ -67,6 +67,11 @@
   return yy_set_m128i(mhi, mlo);
 }
 
+static INLINE void yy_storeu2_128(void *hi, void *lo, const __m256i a) {
+  _mm_storeu_si128((__m128i *)hi, _mm256_extracti128_si256(a, 1));
+  _mm_storeu_si128((__m128i *)lo, _mm256_castsi256_si128(a));
+}
+
 static INLINE __m256i yy_roundn_epu16(__m256i v_val_w, int bits) {
   const __m256i v_s_w = _mm256_srli_epi16(v_val_w, bits - 1);
   return _mm256_avg_epu16(v_s_w, _mm256_setzero_si256());
diff --git a/test/blend_a64_mask_test.cc b/test/blend_a64_mask_test.cc
index 66ca6fc..7592533 100644
--- a/test/blend_a64_mask_test.cc
+++ b/test/blend_a64_mask_test.cc
@@ -86,6 +86,7 @@
     w_ = block_size_wide[block_size];
     h_ = block_size_high[block_size];
     run_times = run_times > 1 ? run_times / w_ : 1;
+    ASSERT_GT(run_times, 0);
     subx_ = subx;
     suby_ = suby;
 
@@ -248,13 +249,13 @@
 INSTANTIATE_TEST_CASE_P(SSE4_1, BlendA64MaskTest8B,
                         ::testing::Values(TestFuncs(
                             aom_blend_a64_mask_c, aom_blend_a64_mask_sse4_1)));
-#endif  // HAVE_AVX2
+#endif  // HAVE_SSE4_1
 
 #if HAVE_AVX2
 INSTANTIATE_TEST_CASE_P(AVX2, BlendA64MaskTest8B,
                         ::testing::Values(TestFuncs(aom_blend_a64_mask_sse4_1,
                                                     aom_blend_a64_mask_avx2)));
-#endif  // HAVE_SSE4_1
+#endif  // HAVE_AVX2
 
 //////////////////////////////////////////////////////////////////////////////
 // 8 bit _d16 version
@@ -482,6 +483,7 @@
   static const int kSrcMaxBitsMaskHBD = (1 << 16) - 1;
 
   void Execute(const uint16_t *p_src0, const uint16_t *p_src1, int run_times) {
+    ASSERT_GT(run_times, 0) << "Cannot run 0 iterations of the test.";
     ConvolveParams conv_params;
     conv_params.round_0 = (bit_depth_ == 12) ? ROUND0_BITS + 2 : ROUND0_BITS;
     conv_params.round_1 = COMPOUND_ROUND1_BITS;
@@ -566,11 +568,45 @@
     }
   }
 }
+TEST_P(BlendA64MaskTestHBD_d16, DISABLED_Speed) {
+  const int kRunTimes = 10000000;
+  for (int bsize = 0; bsize < BLOCK_SIZES_ALL; ++bsize) {
+    for (bit_depth_ = 8; bit_depth_ <= 12; bit_depth_ += 2) {
+      for (int i = 0; i < kBufSize; ++i) {
+        dst_ref_[i] = rng_.Rand12() % (1 << bit_depth_);
+        dst_tst_[i] = rng_.Rand12() % (1 << bit_depth_);
+
+        src0_[i] = rng_.Rand16();
+        src1_[i] = rng_.Rand16();
+      }
+
+      for (int i = 0; i < kMaxMaskSize; ++i)
+        mask_[i] = rng_(AOM_BLEND_A64_MAX_ALPHA + 1);
+
+      RunOneTest(bsize, 1, 1, kRunTimes);
+      RunOneTest(bsize, 0, 0, kRunTimes);
+    }
+  }
+}
 
 INSTANTIATE_TEST_CASE_P(
     C, BlendA64MaskTestHBD_d16,
     ::testing::Values(TestFuncsHBD_d16(aom_highbd_blend_a64_d16_mask_c, NULL)));
 
+#if HAVE_SSE4_1
+INSTANTIATE_TEST_CASE_P(
+    SSE4_1, BlendA64MaskTestHBD_d16,
+    ::testing::Values(TestFuncsHBD_d16(aom_highbd_blend_a64_d16_mask_c,
+                                       aom_highbd_blend_a64_d16_mask_sse4_1)));
+#endif  // HAVE_SSE4_1
+
+#if HAVE_AVX2
+INSTANTIATE_TEST_CASE_P(
+    AVX2, BlendA64MaskTestHBD_d16,
+    ::testing::Values(TestFuncsHBD_d16(aom_highbd_blend_a64_d16_mask_c,
+                                       aom_highbd_blend_a64_d16_mask_avx2)));
+#endif  // HAVE_AVX2
+
 // TODO(slavarnway): Enable the following in the avx2 commit. (56501)
 #if 0
 #if HAVE_AVX2