4

2 つの符号付き 64 ビット整数を乗算し、(128 ビット) 結果を符号付き 64 ビット整数にシフトする必要がありaますb。それを行う最も速い方法は何ですか？

私の 64 ビット整数は、実際には固定小数点数をfmt小数ビットで表しています。fmtオーバーフローしないように が選択さa * b >> fmtれます。たとえばabs(a) < 64<<fmt、abs(b) < 2<<fmtwithfmt==56は最終結果が 64 ビットでオーバーフローしない< 128<<fmtため、int64 に収まります。

私がそれをしたい理由は、固定小数点形式の形式の 5 次多項式を迅速かつ正確に評価するためです。すべての数値は、小数ビット((((c5*x + c4)*x + c3)*x + c2)*x + c1)*x + c0を持つ符号付き 64 ビット固定小数点数です。fmtそれを達成するための最も効率的な方法を探しています。

4

1 に答える 1

8

指摘された質問に対するコメンターとして、これは、移植可能なコードではなく、マシン依存のコードによって最も簡単に効率的に達成されます。質問者は、メイン プラットフォームが x86_64 であり、64 × 64 → 128 ビットの乗算を実行するための組み込み命令があると述べています。これは、インライン アセンブリの小さな部分を使用して簡単にアクセスできます。インライン アセンブリの詳細はコンパイラによって多少異なる場合があることに注意してください。以下のコードはインテル 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;
}

上記のコードに相当する移植可能な 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;
}
于 2015-07-27T20:42:27.473 に答える