ビット演算でダブル型数値を変換―浮動小数点乗算と加算のみで実現（2020年）

トム・レイザーの言葉でいうと、
「これは完全に無意味ですが、いつかあなた方のうち誰かが少し奇妙な状況で役立つかもしれません。」

---

問題設定

以下のような環境を想定します。

• 使える型は IEEE‑754 の倍精度浮動小数点（「double」）だけです。
• そのビット表現に直接アクセスできる操作（C++ のビットキャスト、JavaScript の

  DataView

  等）はありません。

このとき、ある double を 2 個の符号なし 32‑bit 整数（いずれも double として保持）へ変換したり、逆にそれらを組み合わせて元の double に戻したりすることが必要になります。

こうした問題は時折出てきます。そこで興味深かったのは、「浮動小数点演算だけで（加減・乗除・ビット操作・分岐・指数計算・除算・比較を使わずに）」実装できるかということです。

実際、倍精度から整数への変換は「二分探索」的手法で行えます。
まず指数部をビット単位で確定し、その後小数部（仮数）も同様に取得します。ただし、乗算と加算だけで実現できるか？ という疑問は最初は非現実的に思えていましたが、結局 はい。ただし、いくつかの制約があります。

主な制約

• NaN, Infinity, -Infinity は変換不可（ビット操作なしでは区別できないため）。
• 負ゼロと正ゼロは区別不可。加減・乗算だけでは符号情報を保持できません。
• コードは IEEE‑754 の「四捨五入、偶数への丸め」 を前提にしています。
  シェーダーコンパイラなどで

  /fp:fast

  といった高速化フラグが有効になっている環境では動作しません。

以下は実際の実装例です（冗談半分）。
ライブラリに埋め込まれたら面白くなるかもしれません。どんな言語でもポート可能です。

function double_as_uint32s(double) {
  // NaN, Infinity, -Infinity は扱えない。-0 は 0 と同じ扱い。
  var a = double, b, c, d, e, f, g, h, i, j, k, l, m,
      low, high;

  f=2.2250738585072014e-308+a; j=5e-324; b=j+f; b-=f; m=-5e-324; d=m+b;
  b=4.4989137945431964e+161; d=b*d; d=b*d; g=d*d;
  d=1.0; g=d-g; h=m+f; f=h-f; f=j+f; f=b*f; f=b*f; f*=f; f=d-f; f*=g;
  g=-2.2250738585072014e-308+a; h=j+g; h-=g; h=m+h;
  h=b*h; h=b*h; h*=h; h=d-h; c=m+g; c-=g; c=j+c; c=b*c; c=b*c; c*=c; c=d-c; c*=h;
  k=c*f; c=5.562684646268003e-309*a;
  g=j+c; g-=c; g=m+g; g=b*g; g=b*g; g*=g; g=d-g; h=m+c; h-=c; h=j+h; h=b*h; h=b*h; h*=h; h=d-h;
  g=h*g; h=a*g; g=d-g;
  c=g*c; g=1024.0*g; f=2.0+g; c+=h; h=7.458340731200207e-155*c; l=1.0000000000000002;
  g=l*h; g=m+g; e=j+g; e-=g; e=b*e; e=b*e; c=e*c; e=d-e; g=e*h; c=g+c;
  e=512.0*e; g=8.636168555094445e-78*c; e+=f; f=l*g; f=m+f; h=j+f; f=h-f; f=b*f;
  f=b*f; c=f*c; f=d-f; g=f*g; f=256.0*f; c=g+c; e=f+e; f=2.938735877055719e-39*c;
  g=l*f; g=m+g; h=j+g; g=h-g; g=b*g; g=b*... [snipped for brevity] ...

  return [low, high];
}

function uint32s_as_double(low, high) {
  var a = low, b = high, c, d, e, f, g, h, i, j, k, l, m;
  b=9.5367431640625e-07*b; f=-0.4999999998835847; c=f+b; g=4503599627370497.0;
  ... [snipped for brevity] ...
  return a;
}

---

背景：IEEE‑754 の倍精度

ここでは、後述のアルゴリズムを理解するために必要な double の構造と挙動を簡潔に説明します。
既にご存じであればスキップしてください。

double

は 64 ビットです。最上位ビットから順に

1. 符号（1 bit）
2. 指数部（11 bit）
3. 仮数部（52 bit）

を持ちます。このビット列は「特別値」か「数値」に解釈されます。以下の擬似コード（無限精度、

**

はべき乗）で示します。

if (sign == 1) s = -1; else s = 1;

if (exponent == 0x7FF) {
    if (fraction == 0)   return NaN;
    else if (sign == 1)  return -Infinity;
    else                 return Infinity;
} else if (exponent == 0) {        // subnormal
    return s * (0.0 + fraction * (2 ** -52)) * (2 ** -1022);
} else {                           // normal
    return s * (1.0 + fraction * (2 ** -52)) * (2 ** (exponent-1023));
}

重要な性質

項目	内容
符号付き最上位ビット	0/–0 と NaN を除き、すべての数値に対して唯一の表現を保証します。
サブノーマル（非正規化）	零に近づくにつれて表現可能な数の間隔が減少し、ULP が常に可算です。
整数ビット表現との順序性	正の数では、64‑bit 整数表現を大きいほど浮動小数点値も大きくなる（単調増加）。

演算は「正確な数学的操作」として定義され、結果が表現できない場合に丸められます。
IEEE‑754 では複数の丸めモードがありますが、本稿で扱うのは 「最も近い値へ、偶数への丸め」（round to nearest, ties to even）です。

---

2.1 論理演算

加減乗算だけで次の論理プリミティブを構築できます。

constexpr double double_and(double a, double b) { return a * b; }
constexpr double double_not(double a)           { return 1 - a; }
constexpr double double_or (double a, double b){ return a + b - a*b; }

constexpr double select(double cond, double t, double f)
{ return cond*t + double_not(cond)*f; }

これらはすべて有限数に対して機能します。
Infinity や NaN を生成しないよう注意が必要です。

---

2.2 比較の構築

adding_smallest_is_precise(x)

最小正サブノーマル（

p2(-1074)

）を加えると変化するかどうかを判定します。
結果は 0.0 または 1.0 の真偽値として返ります。

is_exp_0_or_1(x)

前述のテストを

x

と

-x

に対して行い、指数部が 0（サブノーマル）か 1（正規化で最小指数）かを判定します。

is_zero(x)

非常に小さな魔法数を足し引きすることで、ゼロ を検出します。

cmp_sub(x, y)

x - y

の符号を保持したまま、Infinity になることなくスケール調整された値を返します。

is_equal(x, y)

cmp_sub

と

is_zero

を組み合わせて等価判定を行います。

---

2.3 指数部の抽出

constexpr double get_encoded_exponent(double v)
{
    double tmp = v;
    double e   = 0;

    for (int test = 1024; test >= 1; test /= 2) {
        double trial      = tmp * p2(-test);
        double too_small  = is_exp_0_or_1(trial);

        tmp = select(too_small, tmp, trial);
        e   += select(too_small, 0.0, (double)test);
    }

    return select(is_exp_0_or_1(v), double_not(is_exp_0(v)), e + 2);
}

このループは指数部をビット単位で減算し、最終的に 0, 1, 2 のいずれかになるまで続きます。

e

に減算量を蓄積しておくことで、実際の指数値が得られます。

---

2.4 符号と仮数部

is_less_than_zero(x)

正規化後に大きな定数を足すことで、負か非負かを判別します。

make_exp_0_or_1(v)

v

を指数が 0 または 1（あるいは 2）になるまでスケールダウンします。
指数が 2 の場合は半分にします。

get_fraction(v)

上記で得た正規化された値から暗黙の 1 を除き、52 ビット分だけ左シフトして仮数部を取得します。

---

2.5 Double → Unsigned 32‑bit 整数 2 個

struct low_high_doubles { double low; double high; };

constexpr low_high_doubles constexpr_double_as_ints(double v)
{
    double sign     = is_less_than_zero(v);
    double exponent = get_encoded_exponent(v);
    double fraction = get_fraction(v);

    // 上位 32 ビット: 符号 (bit 31)、指数 (bits 20‑30),
    // および仮数部の上位 11 ビット
    double high_fraction = small_positive_floor(fraction * p2(-32));
    double high = sign*p2(31) + exponent*p2(20) + high_fraction;

    // 下位 32 ビット: 残りの仮数ビット
    double low = fraction - high_fraction*p2(32);
    return {low, high};
}

---

2.6 Unsigned 32‑bit 整数 2 個 → Double

constexpr double double_from_sign_exp_fraction(
        double sign, double exponent, double fraction)
{
    // 指数がゼロでないか？
    double exp_is_non_zero = double_not(is_zero(exponent));

    // 仮数部を指数 0 に合わせる
    double v = fraction * p2(-1074);

    // 必要なら暗黙の 1 を足す
    v += select(exp_is_non_zero, p2(-1022), 0.0);

    // シフト量計算
    double e = select(exp_is_non_zero, exponent - 1, 0.0);
    e *= p2(-10);                // ビット抽出用に調整

    // 二分探索的指数復元
    for (int test = 1024; test >= 1; test >>= 1) {
        double cond = small_positive_floor(e);
        e = (e - cond)*2;
        if (test == 1024)
            v *= select(cond, p2(512), 1.0); // 未表現の p2(1024) を扱う
        else
            v *= select(cond, p2(test), 1.0);
    }

    // NaN の生成: 指数全ビットが 1 で仮数部がゼロでない場合
    double is_nan = double_and(is_equal(exponent, 2047),
                               double_not(is_zero(fraction)));
    v *= double_not(is_nan);

    // 符号適用
    v *= select(sign, -1.0, 1.0);
    return v;
}

最後に、2 個の整数から符号・指数・仮数部を抽出して上記関数へ渡します。

constexpr double constexpr_ints_as_double(double l, double h)
{
    double exp_and_sign = small_positive_floor(h * p2(-20));
    double sign         = small_positive_floor(h * p2(-31));
    double exponent     = exp_and_sign - sign*p2(11);
    double fraction     = (h - exp_and_sign*p2(20)) * p2(32) + l;

    return double_from_sign_exp_fraction(sign, exponent, fraction);
}

---

3. 最適化と「汚い」Floor

floor

関数は多くの命令を必要としますが、以下のように丸め誤差を利用して下げることが可能です。
0 〜 2³² の範囲で小数点以下 32 ビットまでの精度がある値に対して有効です。

constexpr double dirty_floor(double v)
{
    return v - (0.5 - p2(-33)) + (p2(52)+1) - (p2(52)+1);
}

その他、

is_exp_around_0_or_1

や

unsafe_is_exp_0_or_1

などのヘルパーも簡略化できます。

---

4. まとめ

上記実装は 「加算・乗算だけで」 任意の有限 double 同士の関数を表現できることを示しています。
ただし、以下の点に留意してください。

1. NaN, Infinity, -Infinity は扱えない（ビット操作なしでは区別不可）。
2. 負ゼロと正ゼロは区別できない。
3. IEEE‑754 の「四捨五入、偶数への丸め」 が前提です。
4. シェーダーコンパイラ等で高速化フラグが有効な環境では動作しません。

それでも、多くの制限付き言語や constexpr 文脈などで「ビット表現を取得したい」という場面に対して、加減乗算のみで実装できるという点は興味深い発見です。
ぜひ、ご自身のプロジェクトに応用してみてください。