Multiplicación / división eficiente de dos enteros de 128 bits en x86 (no 64 bits)

Compilador MinGW / GCC
Cuestiones No se permite ningún código GPL / LGPL (GMP o cualquier biblioteca bignum para el caso, es excesivo para este problema, ya que ya tengo implementada la clase).

He construido mi propia 128 bits clase entera grande de tamaño fijo (destinada para su uso en un motor de juego, pero puede generalizarse a cualquier caso de uso) y considero que el rendimiento de las operaciones actuales de multiplicar y dividir es bastante abismal (sí, las he cronometrado, ver a continuación ) y Me gustaría mejorar (o cambiar) los algoritmos que hacen el cálculo de números de bajo nivel.

Cuando se trata de los operadores de multiplicar y dividir, son insoportablemente lentos en comparación con casi todo lo demás en la clase.

Estas son las medidas aproximadas relativas a mi propia computadora:

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

Como puede ver, hacer la multiplicación es mucho, mucho más lento que sumar o restar. La división es aproximadamente 10 veces más lenta que la multiplicación.

Me gustaría mejorar la velocidad de estos dos operadores porque puede haber una gran cantidad de cálculos por cuadro (productos de puntos, varios métodos de detección de colisiones, etc.).

La estructura (métodos omitidos) se parece a:

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

Multiplicació se realiza actualmente utilizando el típico multiplicación larga método (en ensamblado para que pueda atrapar elEDX output) mientras se ignoran las palabras que están fuera de rango (es decir, solo estoy haciendo 10mull 's en comparación con 16).

Divisió usa el shift-restar algoritmo (la velocidad depende del recuento de bits de los operandos). Sin embargo, no se realiza en conjunto. Me pareció un poco demasiado difícil de reunir y decidí dejar que el compilador lo optimice.

He estado en Google por varios días mirando páginas que describen algoritmos comoKaratsuba Multiplication, división de radix alta y División Newton-Rapson pero los símbolos matemáticos están demasiado lejos de mi cabeza. Me gustaría usar algunos de estos métodos avanzados para acelerar mi código, pero primero tendría que traducir el "griego" en algo comprensible.

Para aquellos que puedan considerar mis esfuerzos "optimización prematura"; Considero que este código es un cuello de botella porque las operaciones matemáticas básicas se vuelven lentas. Puedo ignorar estos tipos de optimización en el código de nivel superior, pero este código se llamará / usará lo suficiente como para que importe.

Quisiera sugerencias sobre qué algoritmo debería usar para mejorar la multiplicación y división (si es posible), y una explicación básica (con suerte fácil de entender) sobre cómo funciona el algoritmo sugerido seríaaltament apreciado.

EDIT: Mejoras de multiplicación

Pude mejorar la operación de multiplicación incorporando código en el operador * = y parece lo más rápido posible.

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

Aquí hay un código básico para que lo examine (tenga en cuenta que mis nombres de tipo son realmente diferentes, esto fue editado por simplicidad):

//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));
    }
};

En cuanto a la división, examinar el código no tiene sentido, ya que tendré que cambiar el algoritmo matemático para ver algún beneficio sustancial. La única opción factible parece ser la división de radix alta, pero todavía tengo que resolver (en mi opinión) solocóm funcionará

Respuestas a la pregunta(4)

Su respuesta a la pregunta