Самый быстрый способ умножить два 64-битных int на 128-битные, а затем >> на 64-битные?

Мне нужно умножить два подписанных 64-битных целых числа a и b вместе, а затем сдвинуть (128-разрядный) результат на подписанное 64-битное целое число. Какой самый быстрый способ сделать это?

Мои 64-битные целые числа фактически представляют числа с фиксированной точкой с дробными битами fmt . fmt выбирается так, чтобы a * b >> fmt не переполнялся, например abs(a) < 64<<fmt и abs(b) < 2<<fmt с fmt==56 никогда не переполнялся бы в 64-битах, поскольку конечный результат будет < 128<<fmt и, следовательно, поместится в int64.

Причина, по которой я хочу это сделать, – быстро и точно оценить квинтичные полиномы вида ((((c5*x + c4)*x + c3)*x + c2)*x + c1)*x + c0 в неподвижной точке формат, при каждом номере – подписанное 64-битное число с фиксированной точкой с дробными битами fmt . Я ищу наиболее эффективный способ добиться этого.

    Как отметил комментатор по указанному вопросу, это наиболее легко выполняется с помощью машинного кода, а не портативного кода. Заявитель утверждает, что основная платформа – x86_64 и имеет встроенную инструкцию для выполнения 64 ✕ 64 → 128 бит умножения. Это легко получить, используя небольшой кусок встроенной сборки. Обратите внимание, что детали встроенной сборки могут несколько отличаться от компилятора, код ниже был построен с помощью компилятора Intel C / C ++.

     #include  /* 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; } 

    Ниже приведена переносная C99, эквивалентная приведенному выше коду. Я тестировал это по сравнению с встроенной версией сборки, и никаких несоответствий не было найдено.

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