Use AVX512F and AVX512IFMA to speed up integer multiplication

[user202729]

There was a paper "Accelerating Big Integer Arithmetic Using Intel IFMA Extensions" (see below) that uses AVX512F and AVX512IFMA to multiply small integers. I (asked AI to) benchmark it against flint and get:

  flint_mpn_mul_n_1024      95.94 ns  +/-   0.44 ns
  mul1024_avx512           100.64 ns  +/-   0.82 ns
  mul1024_vpmadd            45.60 ns  +/-   0.43 ns

  flint_mpn_mul_n_2048     294.76 ns  +/-   4.59 ns
  mul2048_avx512           275.32 ns  +/-   3.53 ns
  mul2048_vpmadd           116.30 ns  +/-   1.96 ns

  flint_mpn_mul_n_3072     649.84 ns  +/-  11.00 ns
  mul3072_avx512           595.43 ns  +/-   7.82 ns
  mul3072_vpmadd           192.20 ns  +/-   2.27 ns

  flint_mpn_mul_n_4096     962.09 ns  +/-  14.97 ns
  mul4096_avx512           899.81 ns  +/-  11.22 ns
  mul4096_vpmadd           607.16 ns  +/-  10.44 ns

Looks quite competitive. Note that the suffix is the number of bits, so for example 1024 means 16 limbs input, 32 limbs output.

The paper only provide some fixed sizes, however. More work is needed to adapt it to all sizes.

Note that the repo use GPL v3, but flint is LGPL v3. Might have issue here.

There's a sharp spike at mul4096_vpmadd. No doubt that can be improved.

References:

• https://github.com/vkrasnov/vpmadd
• https://www.semanticscholar.org/paper/Accelerating-Big-Integer-Arithmetic-Using-Intel-Gueron-Krasnov/fd2f45cb7e2631d34874dde307f317a859d59a68
• http://acsel-lab.com/arithmetic/arith23/data/1616a032.pdf

[fredrik-johansson]

• Looks like the AVX512F versions aren't really worth the trouble.

• I think we could figure out how to generate these ourselves.

• We should also think about how to do Karatsuba. It is probably best to stay in a redundant representation so that the additions and subtractions can be done on SIMD vectors rather than with mpn_add_n / mpn_sub_n. Ideally we should only need schoolbook basecases up 8-16 limbs or so; something is off if Karatsuba doesn't win for larger sizes.

• Note that this is highly related to SIMD+Karatsuba polynomial multiplication #2650, and it is possible that mpn and nmod_poly could share most of the multiplication code (differing only in the input splitting and output recombination).

[user202729]

I spent several hours trying to convert the assembly code to C in such a way that both GCC and clang generate efficient assembly code. Not sure how useful this is, since [1] nothing prevents one from writing assembly directly, and [2] it is highly dependent on the compiler. I guess at least I understand the code better this way. (Might also allow inlining and easier code reuse.) The function to propagate carry is quite neat though, I wonder if it's well-known.

The code

#include <immintrin.h>
#include <stdint.h>
#include <string.h>

static const uint16_t permMask3072[32] __attribute__((aligned(64))) = {
    0,1,2,3, 3,4,5,6, 6,7,8,9, 9,10,11,12,
    13,14,15,16, 16,17,18,19, 19,20,21,22, 22,23,24,25
};

/* The radix-conversion shuffle extracts eight consecutive base-2^52 limbs from
   a 52-byte window.  For window w and lane l, the result is
       floor(X / 2^(52*(8*w + l))) mod 2^52.
   The shifts remove the low 0, 4, 8, or 12 bits introduced by the overlapping
   byte-window loads. */
static const uint64_t shiftMask3072[8] __attribute__((aligned(64))) =
    {0,4,8,12,0,4,8,12};

/* Masks and shifts used to convert normalized base-2^52 accumulator limbs back
   to little-endian base-2^64 limbs.  Each pack step gathers bytes spanning one
   or two neighboring 52-bit vectors, shifts them to their final bit positions
   inside 64-bit lanes, and then adds the partial vectors together.  The 255
   entries zero bytes selected from the inactive source in the corresponding
   vpermb/permutex operation. */
static const uint8_t fixMask0[64] __attribute__((aligned(64))) = {
     0, 1, 2, 3, 4, 5, 6, 7,
     8, 9,10,11,12,13,14,15,
    16,17,17,19,20,21,22,23,
    24,25,26,27,28,29,30,31,
    38,39,255,255,255,255,255,48,
    48,49,50,51,52,53,54,55,
    56,57,58,59,60,61,62,63,
    64,65,66,67,68,69,70,71
};

static const uint8_t fixMask1[64] __attribute__((aligned(64))) = {
     8, 9,10,11,12,13,14,15,
    16,17,18,19,20,21,22,23,
    24,25,26,27,28,29,30,31,
    32,33,34,35,36,37,38,39,
    40,41,42,43,44,45,46,47,
    47,255,255,255,255,56,57,58,
    64,65,66,67,68,69,70,71,
    72,73,74,75,76,77,78,79
};

static const uint64_t fixShift0[8] __attribute__((aligned(64))) =
    {0,12,24,36,0,8,20,32};

static const uint64_t fixShift1[8] __attribute__((aligned(64))) =
    {52,40,28,16,4,4,32,20};

static const uint8_t fixMask2[64] __attribute__((aligned(64))) = {
    13,14,15,255,255,255,255,24,
    24,25,26,27,28,29,30,31,
    32,33,34,35,36,37,38,39,
    40,41,42,43,44,45,46,47,
    48,49,50,51,52,53,54,55,
    255,255,255,255,255,255,255,255,
    255,255,255,255,255,255,255,255,
    255,255,255,255,255,255,255,255
};

static const uint8_t fixMask3[64] __attribute__((aligned(64))) = {
    16,17,18,19,20,21,22,23,
    23, 7, 7, 7, 7, 7,32,33,
    40,41,42,43,44,45,46,47,
    48,49,50,51,52,53,54,55,
    56,57,58,59,60,61,62,64,
    255,255,255,255,255,255,255,255,
    255,255,255,255,255,255,255,255,
    255,255,255,255,255,255,255,255
};

static const uint64_t fixShift2[8] __attribute__((aligned(64))) =
    {4,4,16,28,40,64,64,64};

static const uint64_t fixShift3[8] __attribute__((aligned(64))) =
    {8,0,36,24,12,64,64,64};

static inline __attribute__((always_inline)) __m512i
alignr_epi64(__m512i a, __m512i b, int shift)
{
    switch (shift)
    {
    case 0: return b; /* GCC does not fold _mm512_alignr_epi64(a, b, 0) to b; clang does. */
    case 1: return _mm512_alignr_epi64(a, b, 1);
    case 2: return _mm512_alignr_epi64(a, b, 2);
    case 3: return _mm512_alignr_epi64(a, b, 3);
    case 4: return _mm512_alignr_epi64(a, b, 4);
    case 5: return _mm512_alignr_epi64(a, b, 5);
    case 6: return _mm512_alignr_epi64(a, b, 6);
    case 7: return _mm512_alignr_epi64(a, b, 7);
    default: __builtin_unreachable();
    }
}

/* Input: four vectors LO, LOB, HI, HIB, and R10B as the incoming carry bit.
   HI[5..7] and HIB[5..7] must be zero.
   Output: LO, HI, and R10B, where R10B is 0 or 1, such that

       (R10B, HI[0..4], LO) =
           (input HI[0..4], input LO) + (HIB[0..4], LOB) + input R10B

   as a little-endian 13-limb base-2^64 integer. */
static inline __attribute__((always_inline)) void
carry_next_pair(__m512i *lo, __m512i lob, __m512i *hi, __m512i hib,
        uint8_t *r10b)
{
    __m512i minus_one = _mm512_set1_epi64(-1);

    *lo = _mm512_add_epi64(*lo, lob);
    *hi = _mm512_add_epi64(*hi, hib);
    uint8_t klo_c = (uint8_t) _mm512_cmplt_epu64_mask(*lo, lob);
    uint8_t klo_e = (uint8_t) _mm512_cmpeq_epu64_mask(*lo, minus_one);
    uint8_t khi_c = (uint8_t) _mm512_cmplt_epu64_mask(*hi, hib);
    uint8_t khi_e = (uint8_t) _mm512_cmpeq_epu64_mask(*hi, minus_one);
    uint32_t g = (uint32_t) klo_c | ((uint32_t) khi_c << 8);
    uint32_t e = (uint32_t) klo_e | ((uint32_t) khi_e << 8);
    uint32_t s = (g << 1) + (uint8_t) *r10b + e;
    *r10b = (uint8_t) (s >> 13);
    klo_c = (uint8_t) s ^ klo_e;
    khi_c = (uint8_t) (s >> 8) ^ khi_e;
    *lo = _mm512_mask_sub_epi64(*lo, (__mmask8) klo_c, *lo, minus_one);
    *hi = _mm512_mask_sub_epi64(*hi, (__mmask8) khi_c, *hi, minus_one);
}

/* Input: two vectors LO, LOB, and R10B as the incoming carry bit.
   Output: LO, such that

       LO = input LO + LOB + input R10B

   as a little-endian 8-limb base-2^64 integer modulo 2^512.  There is no
   output carry because the 3072-by-3072 product fits exactly in 96 limbs. */
static inline __attribute__((always_inline)) void
carry_last(__m512i *lo, __m512i lob, uint8_t r10b)
{
    __m512i minus_one = _mm512_set1_epi64(-1);

    *lo = _mm512_add_epi64(*lo, lob);
    uint8_t klo_c = (uint8_t) _mm512_cmplt_epu64_mask(*lo, lob);
    uint8_t klo_e = (uint8_t) _mm512_cmpeq_epu64_mask(*lo, minus_one);
    uint32_t s = ((uint32_t) klo_c << 1) + (uint8_t) r10b + klo_e;
    klo_c = (uint8_t) s ^ klo_e;
    *lo = _mm512_mask_sub_epi64(*lo, (__mmask8) klo_c, *lo, minus_one);
}

/* Repack two vectors of consecutive 64-bit accumulator lanes at radix-2^52
   offsets.

   Input: LO and HI hold x[0]..x[15], with x[15] < 2^52, representing
       S = sum x[i] * 2^(52*i).

   Output: LO, HI, LOB, and HIB are overwritten so that
       S = limbs(LO, HI[0..4]) + limbs(LOB, HIB[0..4])  (mod 2^(13*64)),
   where limbs(u, v) is the little-endian integer formed from the eight lanes of
   u followed by the five lanes of v.  HI[5..7] and HIB[5..7] are zero. */
static inline __attribute__((always_inline)) void
pack_pair(__m512i *lo, __m512i *hi, __m512i *lob, __m512i *hib)
{
    __m512i fix1 = _mm512_load_si512((const void *) fixMask1);
    __m512i fix2 = _mm512_load_si512((const void *) fixMask2);
    __m512i fix3 = _mm512_load_si512((const void *) fixMask3);
    __m512i fs0 = _mm512_load_si512((const void *) fixShift0);
    __m512i fs1 = _mm512_load_si512((const void *) fixShift1);
    __m512i fs2 = _mm512_load_si512((const void *) fixShift2);
    __m512i fs3 = _mm512_load_si512((const void *) fixShift3);
    __m512i t0 = _mm512_load_si512((const void *) fixMask0);

    *lob = _mm512_maskz_permutex2var_epi8((__mmask64) 0xFFFFFF83FFFFFFFFull, *lo, t0, *hi);
    *lo = _mm512_maskz_permutex2var_epi8((__mmask64) 0xFFFFE1FFFFFFFFFFull, *lo, fix1, *hi);
    *lob = _mm512_srlv_epi64(*lob, fs0);
    *lo = _mm512_sllv_epi64(*lo, fs1);
    *lo = _mm512_mask_srli_epi32(*lo, (__mmask16) 0x400, *lo, 8);
    t0 = _mm512_maskz_permutexvar_epi8((__mmask64) 0xFFFFFFFFFFFFFF87ull, fix2, *hi);
    __m512i t1 = _mm512_maskz_permutexvar_epi8((__mmask64) 0xFFFFFFFFFFFFC1FFull, fix3, *hi);
    *hi = _mm512_srlv_epi64(t0, fs2);
    *hib = _mm512_sllv_epi64(t1, fs3);
    *hi = _mm512_mask_slli_epi32(*hi, (__mmask16) 0x2, *hi, 8);
}

/* Repack the final single accumulator vector.  This is the live low half of
   pack_pair(LO, zero, LOB, scratch), with the unused high half omitted. */
static inline __attribute__((always_inline)) void
pack_last(__m512i *lo, __m512i *lob)
{
    __m512i fix1 = _mm512_load_si512((const void *) fixMask1);
    __m512i fs0 = _mm512_load_si512((const void *) fixShift0);
    __m512i fs1 = _mm512_load_si512((const void *) fixShift1);
    __m512i t0 = _mm512_load_si512((const void *) fixMask0);

    *lob = _mm512_maskz_permutexvar_epi8((__mmask64) 0xFFFFFF83FFFFFFFFull, t0, *lo);
    *lo = _mm512_maskz_permutexvar_epi8((__mmask64) 0xFFFFE1FFFFFFFFFFull, fix1, *lo);
    *lob = _mm512_srlv_epi64(*lob, fs0);
    *lo = _mm512_sllv_epi64(*lo, fs1);
    *lo = _mm512_mask_srli_epi32(*lo, (__mmask16) 0x400, *lo, 8);
}

/* Shift the current A-window by one 52-bit limb while keeping eight live
   vectors.  Lane 7 of A[i - 1] becomes lane 0 of A[i], lanes 0..6 of A[i]
   move up by one lane, and zero supplies the vacated low lane of A[0]. */
static inline __attribute__((always_inline)) void
align_a8(__m512i A[9])
{
    __m512i zero = _mm512_setzero_si512();

    for (int i = 7; i > 0; i--)
        A[i] = _mm512_alignr_epi64(A[i], A[i - 1], 7);
    A[0] = _mm512_alignr_epi64(A[0], zero, 7);
}

/* Shift the current A-window by one 52-bit limb when the ninth vector is live;
   this is the same operation as align_a8, extended through A[8]. */
static inline __attribute__((always_inline)) void
align_a9(__m512i A[9])
{
    __m512i zero = _mm512_setzero_si512();

    A[8] = _mm512_alignr_epi64(A[8], A[7], 7);
    align_a8(A);
}

/* Expand a 3072-bit input from 48 little-endian 64-bit limbs into 64 base-2^52
   limbs.  If
       X = sum_{i=0}^{47} in[i] * 2^(64*i),
   then for 0 <= i < 64 this stores
       out[i] = floor(X / 2^(52*i)) mod 2^52.
   Limbs beyond bit 3071 are therefore zero except for the high partial limb
   that contains the top bits of X. */
static __attribute__((noinline)) void
convert_64_to_52(uint64_t out[64], const uint64_t in[48])
{
    __m512i p = _mm512_load_si512((const void *) permMask3072);
    __m512i s = _mm512_load_si512((const void *) shiftMask3072);
    __m512i x[8];

    /* Load overlapping 52-byte windows.  The last window is partial because a
       3072-bit value contains exactly 384 bytes. */
    for (int i = 0; i < 7; i++)
        x[i] = _mm512_loadu_si512((const void *) ((const char *) in + 52*i));
    x[7] = _mm512_maskz_loadu_epi16((__mmask32) 0x3ff, (const void *) ((const char *) in + 52*7));

    /* Shuffle and shift each byte window so each 64-bit lane contains one
       52-bit radix limb. */
    for (int i = 0; i < 8; i++)
    {
        x[i] = _mm512_srlv_epi64(_mm512_permutexvar_epi16(p, x[i]), s);
        _mm512_storeu_si512((void *) (out + 8*i), x[i]);
    }
}

/* Multiply two 3072-bit integers in 48 little-endian 64-bit limbs.

   Let
       A = sum_{i=0}^{47} a[i] * 2^(64*i)
       B = sum_{i=0}^{47} b[i] * 2^(64*i).
   On return,
       res[i] = floor(A * B / 2^(64*i)) mod 2^64, 0 <= i < 96.
   All input limbs are converted to local vector/radix-52 storage before the
   first store to res, so res may overlap a and/or b without changing the
   result.

   The vector code is equivalent to this scalar reference:

       ar[i] = floor(A / 2^(52*i)) mod 2^52, 0 <= i < 64
       br[i] = floor(B / 2^(52*i)) mod 2^52, 0 <= i < 64
       c[0..126] = convolution ar * br in radix 2^52
       carry-propagate c in radix 2^52
       convert sum c[i] * 2^(52*i) to 96 radix-2^64 output limbs.

   The implementation keeps the convolution in 15 AVX-512 accumulators and uses
   madd52lo/madd52hi pairs for the low and high 52-bit halves of each product.
   The final section repacks the 52-bit-spaced accumulator stream into 64-bit
   limbs and performs the carries introduced by that repacking. */
void mul3072_vpmadd_c(uint64_t res[96], const uint64_t a[48], const uint64_t b[48])
{
    __m512i zero = _mm512_setzero_si512();
    __m512i T[1];
    __m512i A[9];
    __m512i acc[15];
    uint64_t B[64] __attribute__((aligned(64)));

    /* Convert b once to radix 2^52.  a is converted below into the sliding
       vector window used by the multiplication schedule. */
    convert_64_to_52(B, b);

    __m512i p = _mm512_load_si512((const void *) permMask3072);
    __m512i s = _mm512_load_si512((const void *) shiftMask3072);

    /* Convert all 64 radix-2^52 limbs of a into eight vectors.  A[i] contains
       ar[8*i] through ar[8*i + 7] in its lanes.  A[8] is a zero extension used
       after the sliding window grows to nine vectors. */
    for (int i = 0; i < 7; i++)
        A[i] = _mm512_srlv_epi64(_mm512_permutexvar_epi16(p, _mm512_loadu_si512((const void *) ((const char *) a + 52*i))), s);
    A[7] = _mm512_maskz_loadu_epi16((__mmask32) 0x3ff, (const void *) ((const char *) a + 52*7));
    A[7] = _mm512_srlv_epi64(_mm512_permutexvar_epi16(p, A[7]), s);
    A[8] = zero;

    for (int i = 0; i < 15; i++)
        acc[i] = zero;

    int off = 0;

    /* First phase of the convolution schedule.

       For the current offset off, the j loop broadcasts br[off + 8*j] for
       j = 0..7 and multiplies it by the eight-vector A window.  The low IFMA
       half is accumulated first.  Except for the final pass of this phase, the
       A window is then shifted by one 52-bit limb and the high IFMA half for
       the same br[off + 8*j] values is accumulated.  This covers the columns
       whose active A window fits in eight vectors. */
# pragma GCC unroll 1
    for (int phase = 0; phase < 5; phase++)
    {
# pragma GCC unroll 8
        for (int j = 0; j < 8; j++)
        {
            T[0] = _mm512_set1_epi64((long long) B[off + j*8]);
# pragma GCC unroll 8
            for (int k = 0; k < 8; k++)
                acc[j + k] = _mm512_madd52lo_epu64(acc[j + k], T[0], A[k]);
        }
        if (phase == 4)
            break;
        align_a8(A);
# pragma GCC unroll 8
        for (int j = 0; j < 8; j++)
        {
            T[0] = _mm512_set1_epi64((long long) B[off + j*8]);
# pragma GCC unroll 8
            for (int k = 0; k < 8; k++)
                acc[j + k] = _mm512_madd52hi_epu64(acc[j + k], T[0], A[k]);
        }
        off++;
    }

    /* Second phase of the convolution schedule.

       The remaining columns need a nine-vector A window because the one-limb
       shifts have moved nonzero input through A[8].  Each low-half step uses
       br[off + 8*j] for j = 0..6 and accumulates products with A[0]..A[8]. */
# pragma GCC unroll 1
    for (int phase = 0; phase < 4; phase++)
    {
        if (phase != 0)
        {
# pragma GCC unroll 7
            for (int j = 0; j < 7; j++)
            {
                T[0] = _mm512_set1_epi64((long long) B[off + j*8]);
# pragma GCC unroll 9
                for (int k = 0; k < 9; k++)
                    acc[j + k] = _mm512_madd52lo_epu64(acc[j + k], T[0], A[k]);
            }
        }

        /* Complete the second-phase step by shifting the nine-vector A window by
           one radix limb and accumulating the corresponding high IFMA halves for
           the same br[off + 8*j] values. */
        align_a9(A);
# pragma GCC unroll 7
        for (int j = 0; j < 7; j++)
        {
            T[0] = _mm512_set1_epi64((long long) B[off + j*8]);
# pragma GCC unroll 9
            for (int k = 0; k < 9; k++)
                acc[j + k] = _mm512_madd52hi_epu64(acc[j + k], T[0], A[k]);
        }
        off++;
    }

    /* Each pack_pair call handles the low 832 bits of 16 accumulator lanes at
       radix-2^52 offsets.  The high bits of lanes 0..14 overlap another lane
       inside that same pair, but bits 52..63 of lane 15 are explicit bits in
       the next pair's first radix-52 position.  Move those bits before packing;
       any later carry out from the packed 64-bit additions is handled by
       carry_next_pair through r10b. */
    for (int i = 0; i < 7; i++)
    {
        A[i] = _mm512_srli_epi64(_mm512_alignr_epi64(zero, acc[2*i + 1], 7), 52);
        acc[2*i + 2] = _mm512_add_epi64(acc[2*i + 2], A[i]);
    }

    uint8_t r10b = 0;
    __m512i stage, out;

    /* Repack accumulator pairs 0/1 through 12/13.  Each pair contains 16
       radix-2^52 lanes, i.e. 13 output limbs.  pair_limb_start = 13*i gives
       the output limb offset for pair i; lane_off tells where that 13-limb
       run starts inside the current 8-limb output vector. */
# pragma GCC unroll 7
    for (int i = 0; i < 7; i++)
    {
        const int ai = 2*i;
        const int pair_limb_start = 13*i;
        const int out_vec = pair_limb_start/8;
        const int lane_off = pair_limb_start & 7;
        const int first_take = 8 - lane_off;
        const int remaining = 13 - first_take;

        pack_pair(&acc[ai], &acc[ai + 1], &A[0], &A[1]);
        carry_next_pair(&acc[ai], A[0], &acc[ai + 1], A[1],
                &r10b);

        if (lane_off == 0)
        {
            _mm512_storeu_si512((void *) (res + 8*out_vec), acc[ai]);
            stage = alignr_epi64(acc[ai + 1], acc[ai + 1], 5);
        }
        else
        {
            const int shift = 8 - lane_off;

            out = alignr_epi64(acc[ai], stage, shift);
            _mm512_storeu_si512((void *) (res + 8*out_vec), out);

            if (remaining >= 8)
            {
                out = alignr_epi64(acc[ai + 1], acc[ai], shift);
                _mm512_storeu_si512((void *) (res + 8*(out_vec + 1)), out);
                stage = alignr_epi64(acc[ai + 1], acc[ai + 1], 5);
            }
            else
            {
                stage = alignr_epi64(acc[ai + 1], acc[ai], 5);
            }
        }
    }

    /* Repack the single remaining accumulator acc[14], propagate the final
       carry mask, merge it with the staged low part in acc[11], and write the
       top product limbs res[88]..res[95]. */
    pack_last(&acc[14], &A[0]);
    carry_last(&acc[14], A[0], r10b);
    out = alignr_epi64(acc[14], stage, 5);
    _mm512_storeu_si512((void *) (res + 8*11), out);
}

With some AI assistance. Not entirely sure how much the AI helps versus messes things up though.

---

Anyway, using the same idea to implement addition gives a 2x speed up over the current mpz_add for large operand size (which is what fmpz_add uses under the hood).

The addition kernel

static inline mp_limb_t
add_n_avx512_16(mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
{
    const __m512i minus_one = _mm512_set1_epi64((long long) -1);
    mp_limb_t carry = 0;
    mp_size_t i = 0;

    for (; i + 16 <= n; i += 16)
    {
        __m512i alo = _mm512_loadu_si512((const void *) (ap + i));
        __m512i ahi = _mm512_loadu_si512((const void *) (ap + i + 8));
        __m512i blo = _mm512_loadu_si512((const void *) (bp + i));
        __m512i bhi = _mm512_loadu_si512((const void *) (bp + i + 8));
        __m512i rlo = _mm512_add_epi64(alo, blo);
        __m512i rhi = _mm512_add_epi64(ahi, bhi);

        uint32_t glo = (uint32_t) _mm512_cmplt_epu64_mask(rlo, blo);
        uint32_t ghi = (uint32_t) _mm512_cmplt_epu64_mask(rhi, bhi);
        uint32_t elo = (uint32_t) _mm512_cmpeq_epu64_mask(rlo, minus_one);
        uint32_t ehi = (uint32_t) _mm512_cmpeq_epu64_mask(rhi, minus_one);
        uint32_t g = glo | (ghi << 8);
        uint32_t e = elo | (ehi << 8);
        uint32_t s = (g << 1) + (uint32_t) carry + e;
        uint32_t carry_mask = s ^ e;

        carry = s >> 16;
        rlo = _mm512_mask_sub_epi64(rlo, (__mmask8) carry_mask, rlo, minus_one);
        rhi = _mm512_mask_sub_epi64(rhi, (__mmask8) (carry_mask >> 8), rhi, minus_one);

        _mm512_storeu_si512((void *) (rp + i), rlo);
        _mm512_storeu_si512((void *) (rp + i + 8), rhi);
    }

    for (; i < n; i++)
    {
        mp_limb_t a = ap[i];
        mp_limb_t b = bp[i];
        mp_limb_t s = a + b;
        mp_limb_t c1 = s < a;
        mp_limb_t r = s + carry;
        mp_limb_t c2 = r < s;
        rp[i] = r;
        carry = c1 | c2;
    }

    return carry;
}

I tried using 4 of __m512i as well as 2, doesn't make much difference. Benchmark also shows prefetching can improve the performance of the ADC chain beyond 2048 limbs (this one doesn't need new processor).

[fredrik-johansson]

Is that the same algorithm for addition as the one described here: https://www.numberworld.org/y-cruncher/internals/addition.html?

[user202729]

If you look closer at the iteration over the C and M bits, you'll notice that each step is equivalent to a full-adder. In other words, the entire iteration is just an integer addition of p-bit integers.

Yes, looks correct. Nice article.