Computing alto 64 bits de un producto int 64x64 en C

Question

Me gustaría que mi función C calcule eficientemente los 64 bits altos del producto de dos entradas firmadas de 64 bits. Sé cómo hacer esto en el ensamblaje x86-64, con imulq y sacando el resultado de% rdx. Pero no sé cómo escribir esto en C, y mucho menos convencer al compilador para que lo haga de manera eficiente.

¿Alguien tiene alguna sugerencia para escribir esto en C? Esto es sensible al rendimiento, por lo que los "métodos manuales" (como Russian Peasant, o bignum libraries) están fuera.

Esta función ensamblado en línea dorky escribí funciona y es más o menos la codegen que busco:

static long mull_hi(long inp1, long inp2) { 
    long output = -1; 
    __asm__("movq %[inp1], %%rax;" 
      "imulq %[inp2];" 
      "movq %%rdx, %[output];" 
      : [output] "=r" (output) 
      : [inp1] "r" (inp1), [inp2] "r" (inp2) 
      :"%rax", "%rdx"); 
    return output; 
} 

Answer 1

La respuesta general es que x * y se puede dividir en (a + b) * (c + d), donde a y c son el alto orden partes.

En primer lugar, se amplía a ac + ad + bc + bd

Ahora, se multiplican los términos como números de 32 bits almacenados como long long (o mejor aún, uint64_t), y que acaba de recordar que cuando se multiplica por el número de orden superior, es necesario escala por 32 bits. Luego haces las adiciones, recordando detectar carry. Mantenga un registro del signo. Naturalmente, necesitas hacer las adiciones en pedazos.

Answer 2

Espere, usted tiene una solución de montaje en perfecto estado, ya optimizado trabajando para esto, y desea realizar una copia hacia fuera y tratar de escribirlo en un entorno que no admite 128 bits de matemáticas? No estoy siguiendo.

Como es obvio, esta operación es una instrucción simple en x86-64. Obviamente, nada de lo que haga va a hacer que funcione mejor. Si realmente quieres una C portátil, deberás hacer algo como el código de DigitalRoss arriba y esperar que tu optimizador descubra qué estás haciendo .

Si necesita portabilidad arquitectura pero está dispuesto a limitarse a las plataformas gcc, hay __int128_t (y __uint128_t) en los tipos intrínsecos del compilador que va a hacer lo que quiere.

Answer 3

Si está utilizando una relativamente reciente GCC en x86_64:

int64_t mulHi(int64_t x, int64_t y) { 
    return (int64_t)((__int128_t)x*y >> 64); 
} 

En -O1 y superior, esto se compila a lo que quiere:

_mulHi: 
0000000000000000 movq %rsi,%rax 
0000000000000003 imulq %rdi 
0000000000000006 movq %rdx,%rax 
0000000000000009 ret 

creo que sonido metálico y VC++ también tiene soporte para el tipo __int128_t, por lo que esto también debería funcionar en esas plataformas, con las advertencias habituales sobre intentarlo usted mismo.

Answer 4

Con respecto a su solución de ensamblaje, ¡no codifique las instrucciones mov! Deje que el compilador lo haga por usted. Aquí está una versión modificada de su código:

static long mull_hi(long inp1, long inp2) { 
    long output; 
    __asm__("imulq %2" 
      : "=d" (output) 
      : "a" (inp1), "r" (inp2)); 
    return output; 
} 

referencia útil: Machine Constraints

Answer 5

Desde que hizo un muy buen trabajo para resolver su propio problema con el código máquina, que pensé que merecía un poco de ayuda con la versión portátil.Dejaría un ifdef en donde solo use el ensamblado si está en gnu en x86.

De todos modos, aquí hay una implementación ... Estoy bastante seguro de que esto es correcto, pero no hay garantías, acabo de decir esto anoche ... probablemente deberías deshacerte de las estáticas positive_result [] y result_negative, esos son solo artefactos de mi unidad de prueba ...

#include <stdlib.h> 
#include <stdio.h> 

// stdarg.h doesn't help much here because we need to call llabs() 

typedef unsigned long long uint64_t; 
typedef signed long long int64_t; 

#define B32 0xffffffffUL 

static uint64_t positive_result[2]; // used for testing 
static int result_negative;   // used for testing 

static void mixed(uint64_t *result, uint64_t innerTerm) 
{ 
    // the high part of innerTerm is actually the easy part 

    result[1] += innerTerm >> 32; 

    // the low order a*d might carry out of the low order result 

    uint64_t was = result[0]; 

    result[0] += (innerTerm & B32) << 32; 

    if (result[0] < was) // carry! 
     ++result[1]; 
} 


static uint64_t negate(uint64_t *result) 
{ 
    uint64_t t = result[0] = ~result[0]; 
    result[1] = ~result[1]; 
    if (++result[0] < t) 
    ++result[1]; 
    return result[1]; 
} 

uint64_t higherMul(int64_t sx, int64_t sy) 
{ 
    uint64_t x, y, result[2] = { 0 }, a, b, c, d; 

    x = (uint64_t)llabs(sx); 
    y = (uint64_t)llabs(sy); 

    a = x >> 32; 
    b = x & B32; 
    c = y >> 32; 
    d = y & B32; 

    // the highest and lowest order terms are easy 

    result[1] = a * c; 
    result[0] = b * d; 

    // now have the mixed terms ad + bc to worry about 

    mixed(result, a * d); 
    mixed(result, b * c); 

    // now deal with the sign 

    positive_result[0] = result[0]; 
    positive_result[1] = result[1]; 
    result_negative = sx < 0^sy < 0; 
    return result_negative ? negate(result) : result[1]; 
}