x86 上两个 128 位整数的高效乘法/除法（无 64 位）

编译器：明威/海湾合作委员会
Issues:不允许使用 GPL/LGPL 代码（GMP 或任何 bignum 库对于这个问题来说都太过分了，因为我已经实现了该类）。

我已经构建了自己的128-bit固定大小的大整数类（旨在用于游戏引擎，但可以推广到任何使用情况），我发现当前乘法和除法运算的性能非常糟糕（是的，我已经对它们进行了计时，见下文） ， 和我想改进（或更改）进行低级数字运算的算法。

当涉及乘法和除法运算符时，与类中的其他所有运算符相比，它们慢得难以忍受。

这些是相对于我自己的计算机的近似测量值：

Raw times as defined by QueryPerformanceFrequency:
1/60sec          31080833u
Addition:              ~8u
Subtraction:           ~8u
Multiplication:      ~546u
Division:           ~4760u (with maximum bit count)

正如您所看到的，仅进行乘法比加法或减法慢很多很多倍。除法比乘法慢大约 10 倍。

我想提高这两个运算符的速度，因为每帧可能需要进行大量计算（点积、各种碰撞检测方法等）。

结构（省略方法）看起来有点像：

class uint128_t
{
    public:
        unsigned long int dw3, dw2, dw1, dw0;
  //...
}

乘法目前使用典型的长乘法方法（在汇编中，以便我可以捕获EDX输出），同时忽略超出范围的单词（也就是说，我只做了 10mull与 16 相比）。

Division使用移位减法算法（速度取决于操作数的位数）。然而，它不是在装配中完成的。我发现这有点太难了，决定让编译器优化它。

我在谷歌上搜索了几天，查看描述算法的页面，例如唐叶乘法 http://en.wikipedia.org/wiki/Karatsuba_algorithm，高基数除法，以及牛顿-拉夫逊分部 http://en.wikipedia.org/wiki/Division_%28digital%29#Newton.E2.80.93Raphson_division但数学符号有点超出我的理解范围。我想使用其中一些高级方法来加速我的代码，但我必须首先将“希腊语”翻译成可以理解的内容。

对于那些可能认为我的努力“过早优化”的人；我认为这段代码是一个瓶颈，因为非常基本的数学运算本身变得很慢。我可以忽略对更高级别代码的此类优化，但该代码将被足够多的调用/使用以使其发挥作用。

我想要关于应该使用哪种算法来改进乘法和除法（如果可能）的建议，以及关于建议算法如何工作的基本（希望易于理解）解释highly赞赏。

编辑：乘以改进

我能够通过将代码内联到运算符 *= 中来改进乘法运算，并且它看起来尽可能快。

Updated raw times:
1/60sec          31080833u
Addition:              ~8u
Subtraction:           ~8u
Multiplication:      ~100u (lowest ~86u, highest around ~256u)
Division:           ~4760u (with maximum bit count)

这里有一些简单的代码供您检查（请注意，我的类型名称实际上不同，为了简单起见，对其进行了编辑）：

//File: "int128_t.h"
class int128_t
{
    uint32_t dw3, dw2, dw1, dw0;

    // Various constrctors, operators, etc...

    int128_t& operator*=(const int128_t&  rhs) __attribute__((always_inline))
    {
        int128_t Urhs(rhs);
        uint32_t lhs_xor_mask = (int32_t(dw3) >> 31);
        uint32_t rhs_xor_mask = (int32_t(Urhs.dw3) >> 31);
        uint32_t result_xor_mask = (lhs_xor_mask ^ rhs_xor_mask);
        dw0 ^= lhs_xor_mask;
        dw1 ^= lhs_xor_mask;
        dw2 ^= lhs_xor_mask;
        dw3 ^= lhs_xor_mask;
        Urhs.dw0 ^= rhs_xor_mask;
        Urhs.dw1 ^= rhs_xor_mask;
        Urhs.dw2 ^= rhs_xor_mask;
        Urhs.dw3 ^= rhs_xor_mask;
        *this += (lhs_xor_mask & 1);
        Urhs += (rhs_xor_mask & 1);

        struct mul128_t
        {
            int128_t dqw1, dqw0;
            mul128_t(const int128_t& dqw1, const int128_t& dqw0): dqw1(dqw1), dqw0(dqw0){}
        };

        mul128_t data(Urhs,*this);
        asm volatile(
        "push      %%ebp                            \n\
        movl       %%eax,   %%ebp                   \n\
        movl       $0x00,   %%ebx                   \n\
        movl       $0x00,   %%ecx                   \n\
        movl       $0x00,   %%esi                   \n\
        movl       $0x00,   %%edi                   \n\
        movl   28(%%ebp),   %%eax #Calc: (dw0*dw0)  \n\
        mull             12(%%ebp)                  \n\
        addl       %%eax,   %%ebx                   \n\
        adcl       %%edx,   %%ecx                   \n\
        adcl       $0x00,   %%esi                   \n\
        adcl       $0x00,   %%edi                   \n\
        movl   24(%%ebp),   %%eax #Calc: (dw1*dw0)  \n\
        mull             12(%%ebp)                  \n\
        addl       %%eax,   %%ecx                   \n\
        adcl       %%edx,   %%esi                   \n\
        adcl       $0x00,   %%edi                   \n\
        movl   20(%%ebp),   %%eax #Calc: (dw2*dw0)  \n\
        mull             12(%%ebp)                  \n\
        addl       %%eax,   %%esi                   \n\
        adcl       %%edx,   %%edi                   \n\
        movl   16(%%ebp),   %%eax #Calc: (dw3*dw0)  \n\
        mull             12(%%ebp)                  \n\
        addl       %%eax,   %%edi                   \n\
        movl   28(%%ebp),   %%eax #Calc: (dw0*dw1)  \n\
        mull              8(%%ebp)                  \n\
        addl       %%eax,   %%ecx                   \n\
        adcl       %%edx,   %%esi                   \n\
        adcl       $0x00,   %%edi                   \n\
        movl   24(%%ebp),   %%eax #Calc: (dw1*dw1)  \n\
        mull              8(%%ebp)                  \n\
        addl       %%eax,   %%esi                   \n\
        adcl       %%edx,   %%edi                   \n\
        movl   20(%%ebp),   %%eax #Calc: (dw2*dw1)  \n\
        mull              8(%%ebp)                  \n\
        addl       %%eax,   %%edi                   \n\
        movl   28(%%ebp),   %%eax #Calc: (dw0*dw2)  \n\
        mull              4(%%ebp)                  \n\
        addl       %%eax,   %%esi                   \n\
        adcl       %%edx,   %%edi                   \n\
        movl   24(%%ebp),  %%eax #Calc: (dw1*dw2)   \n\
        mull              4(%%ebp)                  \n\
        addl       %%eax,   %%edi                   \n\
        movl   28(%%ebp),   %%eax #Calc: (dw0*dw3)  \n\
        mull               (%%ebp)                  \n\
        addl       %%eax,   %%edi                   \n\
        pop        %%ebp                            \n"
        :"=b"(this->dw0),"=c"(this->dw1),"=S"(this->dw2),"=D"(this->dw3)
        :"a"(&data):"%ebp");

        dw0 ^= result_xor_mask;
        dw1 ^= result_xor_mask;
        dw2 ^= result_xor_mask;
        dw3 ^= result_xor_mask;
        return (*this += (result_xor_mask & 1));
    }
};

至于除法，检查代码是毫无意义的，因为我需要更改数学算法才能看到任何实质性的好处。唯一可行的选择似乎是高基数除法，但我还没有解决（在我看来）how它会起作用的。

我不会太担心乘法。你正在做的事情看起来非常有效。我并没有真正遵循希腊语的唐叶乘法，但我的感觉是，只有比你正在处理的数字大得多时，它才会更有效。

我的一个建议是尝试使用最小的内联汇编块，而不是在汇编中编码逻辑。你可以写一个函数：

struct div_result { u_int x[2]; };
static inline void mul_add(int a, int b, struct div_result *res);

该函数将在内联汇编中实现，您将从 C++ 代码中调用它。它应该与纯汇编一样高效，并且更容易编码。

至于分工，我不知道。我看到的大多数算法都谈论渐近效率，这可能意味着它们仅对非常高的位数有效。