Necesito multiplicar dos enteros a y b de 64 bits con signo, luego cambiar el resultado (128 bits) a un entero de 64 bits con signo. ¿Cuál es la forma más rápida de hacer eso?
Mis números enteros de 64 bits en realidad representan números de punto fijo con bits fraccionarios fmt . fmt se elige para que a * b >> fmt no se desborde, por ejemplo abs(a) < 64<<fmt y abs(b) < 2<<fmt con fmt==56 nunca se desbordaría en 64 bits como el el resultado final sería < 128<<fmt y por lo tanto cabría en un int64.
La razón por la que quiero hacer eso es para evaluar de forma rápida y precisa los polinomios quínticos de la forma ((((c5*x + c4)*x + c3)*x + c2)*x + c1)*x + c0 en punto fijo formato, con cada número un número de punto fijo de 64 bits con signo con bits fraccionarios fmt . Estoy buscando la forma más eficiente de lograrlo.
Como señaló un comentarista de la pregunta, esto se logra de manera más fácil y eficiente mediante un código dependiente de la máquina, en lugar de un código portátil. El autor de la pregunta afirma que la plataforma principal es x86_64 y que tiene una instrucción integrada para realizar una multiplicación de 64 × 64 → 128 bits. Se accede fácilmente a esto usando una pequeña pieza de ensamblaje en línea. Tenga en cuenta que los detalles del ensamblaje en línea pueden diferir un poco con el compilador, el siguiente código se creó con el compilador Intel C/C++.
#include <stdint.h> /* compute mul_wide (a, b) >> s, for s in [0,63] */ int64_t mulshift (int64_t a, int64_t b, int s) { int64_t res; __asm__ ( "movq %1, %%rax;\n\t" // rax = a "movl %3, %%ecx;\n\t" // ecx = s "imulq %2;\n\t" // rdx:rax = a * b "shrdq %%cl, %%rdx, %%rax;\n\t" // rax = int64_t (rdx:rax >> s) "movq %%rax, %0;\n\t" // res = rax : "=rm" (res) : "rm"(a), "rm"(b), "rm"(s) : "%rax", "%rdx", "%ecx"); return res; }A continuación se muestra un C99 portátil equivalente al código anterior. He probado esto extensamente contra la versión de ensamblaje en línea y no se encontraron discrepancias.
void umul64wide (uint64_t a, uint64_t b, uint64_t *hi, uint64_t *lo) { uint64_t a_lo = (uint64_t)(uint32_t)a; uint64_t a_hi = a >> 32; uint64_t b_lo = (uint64_t)(uint32_t)b; uint64_t b_hi = b >> 32; uint64_t p0 = a_lo * b_lo; uint64_t p1 = a_lo * b_hi; uint64_t p2 = a_hi * b_lo; uint64_t p3 = a_hi * b_hi; uint32_t cy = (uint32_t)(((p0 >> 32) + (uint32_t)p1 + (uint32_t)p2) >> 32); *lo = p0 + (p1 << 32) + (p2 << 32); *hi = p3 + (p1 >> 32) + (p2 >> 32) + cy; } void mul64wide (int64_t a, int64_t b, int64_t *hi, int64_t *lo) { umul64wide ((uint64_t)a, (uint64_t)b, (uint64_t *)hi, (uint64_t *)lo); if (a < 0LL) *hi -= b; if (b < 0LL) *hi -= a; } /* compute mul_wide (a, b) >> s, for s in [0,63] */ int64_t mulshift (int64_t a, int64_t b, int s) { int64_t res; int64_t hi, lo; mul64wide (a, b, &hi, &lo); if (s) { res = ((uint64_t)hi << (64 - s)) | ((uint64_t)lo >> s); } else { res = lo; } return res; }