algorithm - 短长度循环卷积的优化

Question

我有两个 8 个无符号字节的序列，我需要计算它们的循环卷积，这会产生 8 个无符号的 19 位整数。当我重复这一百万次时，我想进行优化。

直接的方法需要 64 个 MAC 操作。我已经知道如何通过 SSE/AVX 指令来加速它，但这不是我所追求的。

是否有另一种方法，可能基于 FFT 或数论变换来减少运算次数或其他技术来获得一些加速？

其实我不需要这8个值:最大的一个和相应的移位就够了。

Answer 1

我在研究一个相关的问题，发现了这个问题，并对它产生了好奇。

乘法可以看作是具有额外步骤的卷积，如果考虑到适当的零填充以允许进位等，则不难将一个嵌入另一个。

在这个具体问题中，可以把8个字节的循环卷积看成乘法模，即C = A*B mod M，其中:

A = a0 + a1 * K + a1 * K^2 + ...
B = b0 + b7 * K + b6 * K^2 + ...
C = c0 + c1 * K + c2 * K^2 + ...
M = K^8 - 1

选择的 K > 8 * 255 * 255 足够大以表示卷积每个点的最大可能值，而进位不会溢出下一个值。

在实践中，我们将始终选择 K=2^k，以便嵌入和提取单个值对于使用移位来说是微不足道的，因此 k>=19。

有多种方法可以在现代硬件中高效地实现这种任意精度的产品。在现代 CPU 上实现 64 位整数乘法在延迟和吞吐量方面与较短长度的乘法一样快，因此我们可以尝试将我们的产品塞入尽可能少的 64 位字中，在本例中为 N =3，从 64N >= 8k 到 19<=k<=24。

执行此操作的简单方法是使用 9 MULQ(每个 MULQ 仅比现代硬件上的常规 64 位乘法慢一点)，而 Toom-Cook 乘法仅使用 5 MULQ。理论上有一种方法可以使用 Winograd 循环卷积对此进行改进，方法是将 M 拆分为适当的因子，计算对这些因子取模的乘积并重建解决方案。原则上，这可以将操作数进一步减少到仅 4 或 3 个 MULQ，但代价是显着增加了位操作技巧以避免代价高昂的除法。

下面的代码给出了预期的结果:

#include <iostream>
#include <cinttypes>
#include <cassert>
using namespace std;

static const int NN     = 8; // length of the convolution
static const int LGNN   = 3; // ceil(log2(NN))
static const int PAD    = 2*(8*sizeof(uint8_t)) + LGNN; // == 19
static const int BPW    = (8*sizeof(uint64_t)) / PAD; // == 3 // bytes per word
static const int NW     = (NN+BPW-1)/BPW ; // number of 64-bit words to store the sequence


// 5-point toom-cook linear convolution
void cyclic5 (uint32_t c[], const uint8_t a[], const uint8_t b[]) {

    static const uint64_t WMASK = (1ull<<BPW*PAD) -1; //0x01FFFFFFFFFFFFFF;
    static const uint32_t CMASK = (1ul <<    PAD) -1; //0x0007FFFF;

    uint64_t A[NW]={0,0,0}, B[NW]={0,0,0};
    __uint128_t Q[2*NW-1];


    for (int k=0; k<NN; ++k) A[k/BPW] |= uint64_t(a[k])<<((k%BPW)*PAD);
    for (int k=0; k<NN; ++k) B[k/BPW] |= uint64_t(b[(NN-k)%NN])<<((k%BPW)*PAD);

    // 5-point toom-cook
    // 5 mul, 12 add, 6 shift
    {
        // x==0
        Q[0] = (__uint128_t)(A[0])               * (__uint128_t)(B[0]);
        // x==2
        Q[1] = (__uint128_t)(A[0]+2*A[1]+4*A[2]) * (__uint128_t)(B[0]+2*B[1]+4*B[2]);
        // x==1
        Q[2] = (__uint128_t)(A[0]+  A[1]+  A[2]) * (__uint128_t)(B[0]+  B[1]+  B[2]);
        // x==1/2
        Q[3] = (__uint128_t)(A[2]+2*A[1]+4*A[0]) * (__uint128_t)(B[2]+2*B[1]+4*B[0]);
        // x==infinity
        Q[4] = (__uint128_t)(A[2])               * (__uint128_t)(B[2]);
    }

    // recover the partial products by polynomial interpolation
    {
        // 6 add, 4 shift
        Q[1] -= Q[0]; Q[1] -= 16*Q[4]; Q[1]>>=1;
        Q[3] -= Q[4]; Q[3] -= 16*Q[0]; Q[3]>>=1;
        Q[2] -= Q[0]; Q[2] -=    Q[4];

        // 5 add, 2 shift
        // 2 add, 2 shift, 1 di!
        {
            Q[2]  = 5*Q[2]  - Q[1] - Q[3];
            Q[1] -= 2*Q[2];
            Q[3] -= 2*Q[2];

            __uint128_t Q1 = (4*Q[3]-Q[1])/15; //  division by 15 can be replaced by a few additions and shifts
            __uint128_t Q3 = Q[3] - 4*Q1;

            Q[1] = Q1;
            Q[3] = Q3;
        }

    }

    uint64_t P[2*NW-1 +1];

    for (int w=0; w<2*NW-1; ++w) P[w  ]  = (uint64_t)(Q[w]) & WMASK;
    for (int w=0; w<2*NW-2; ++w) P[w+1] += (uint64_t)(Q[w] >> (BPW*PAD));


    // wrap around and extract
    for (int k=0; k<2*NN-1; ++k) {
        uint32_t v = (P[k/BPW] >> (k%BPW)*PAD) & CMASK;
        if (k<NN) c[k]     = v;
        else      c[k%NN] += v;
    }
}


// "naive" 3x3-term multiplication
void cyclic9 (uint32_t c[], const uint8_t a[], const uint8_t b[]) {

    static const uint64_t WMASK = (1ull<<BPW*PAD) -1; //0x01FFFFFFFFFFFFFF;
    static const uint32_t CMASK = (1ul <<    PAD) -1; //0x0007FFFF;

    uint64_t A[NW]={0,0,0}, B[NW]={0,0,0};
    uint64_t P[2*NW-1 +1] = {0,0,0,0,0,0};

    for (int k=0; k<NN; ++k) A[k/BPW] |= uint64_t(a[k])<<((k%BPW)*PAD);
    for (int k=0; k<NN; ++k) B[k/BPW] |= uint64_t(b[(NN-k)%NN])<<((k%BPW)*PAD);

    // linear convolution aka product, with the standard multiplication algorithm
    for (int i=0; i<NW; ++i)
        for (int j=0; j<NW; ++j) {
            // will compile to a single MULQ in x86-64 assembly
            __uint128_t Q = (__uint128_t)A[i] * (__uint128_t)B[j];

            P[i+j  ] += (uint64_t)(Q) & WMASK;
            P[i+j+1] += (uint64_t)(Q >> (BPW*PAD));
        }

    // wrap around and extract
    for (int k=0; k<2*NN-1; ++k) {
        uint32_t v = (P[k/BPW] >> (k%BPW)*PAD) & CMASK;
        if (k<NN) c[k]     = v;
        else      c[k%NN] += v;
    }
}

// basic 64x 16-bit FMA
void cyclic64 (uint32_t c[], const uint8_t a[], const uint8_t b[]) {

    for (int k=0; k<NN; ++k) c[k] = 0;

    // linear convolution, aka product
    for (int i=0; i<NN; ++i)
        for (int j=0; j<NN; ++j)
            c[(NN+i-j)%NN] += uint16_t(a[i])*uint16_t(b[j]);
}

int main() {

    assert (PAD==19);
    assert (BPW==3);
    assert (NW==3);       

    uint8_t B[] = {0x92, 0x16, 0x5e, 0xf1, 0x29, 0xe9, 0xbb, 0xcd, 0x8e, 0xd2, 0xa8, 0xb8, 0xae, 0xc1, 0x4a, 0x60};

    uint32_t R[NN];

    cyclic64(R, B, B+NN);
    for (int i=0; i<NN; ++i) cout << R[i] << "  "; cout << endl;

    cyclic9(R, B, B+NN);
    for (int i=0; i<NN; ++i) cout << R[i] << "  "; cout << endl;

    cyclic5(R, B, B+NN);
    for (int i=0; i<NN; ++i) cout << R[i] << "  "; cout << endl;

    return 0;
}