int n = (int)floor(x) - 1;
 float d = x - n;

 e = exp(n - 1.0) * exp(d0 + 1) * exp(d1 - 1);

 e = exp((n - 1.0) / 8) * exp((d0 + 1) / 8) * exp((d1 - 1) / 8);

  fi.i =0x3f800000 | α // 0x3f800000は1.0fのバイナリ表現
  d1 = fi.f;

  fi.i =0x3f800000 | (α >> 3);
  d1 = fi.f;

    for (float x = 0; x <= 20; x+= 1e-5f) {
        sum += exp(x);
    }

// __m128 d;
float e;
*(int*)&e = _mm_cvtsi128_si32(_mm_castps_si128(d));

>最後の2行はx86ではfloatはFPUレジスタに入れて返すのでこうするしかない.

float e = _mm_cvtss_f32(d);

float Exp2(float x)
{
static const float one = 1.0f;
static const unsigned int mask2 = (mask(8) << 23) | mask(11);

// float xの値をx1に代入
__m128 x1 = _mm_load_ss(&x);
// floor(x1)を求めて1を引く
int n = _mm_cvttss_si32(x1) - 1;

__m128 d = _mm_setzero_ps();
d = _mm_sub_ss(x1, _mm_cvtsi32_ss(d, n));

__m128i di = _mm_castps_si128(d);
unsigned int d0i = (unsigned int)_mm_cvtsi128_si32(di) >> 11 & mask(12);

d = _mm_and_ps(d, _mm_load_ss((float*)&mask2));
d = _mm_mul_ss(d, _mm_load_ss(&expTbl0[n + 1]));
d = _mm_mul_ss(d, _mm_load_ss(&expTbl1[d0i]));
float e;
*(int*)&e = _mm_cvtsi128_si32(_mm_castps_si128(d));

• 下記の方法で分割したd3は仮数の上位16bitは0(下位7bitしか残ってない)なので非常に1に近い. d3 - 1 < 2^-16. これならx << 1のときのe^x = 1 + x + x²が使える. 今の場合x²は無視できるのでexp(d3 - 1) ～ 1 + (d3 - 1) = d3.
• 仮数部は23bitしかないのだから2^-12より小さいと2乗の項を無視できる. 逆に言うとfは上位12bitしかいらない. そのぐらいなら1回でテーブル引きしてしまってもよいかもしれない.
• ここまでの考察を実装して試してみる.

float Exp(float x)
{
    bool isPositive;
    if (x >= 0) {
        isPositive = true;
    } else {
        x = -x;
        isPositive = true;
    }

    int n = (int)floor(x) - 1;
    float d = x - n;

    union fi {
        float f;
        unsigned int i;
    };
    fi fi;
    fi.f = d;
    unsigned int di = fi.i; // これで di は d のビットパターンになる

    fi.i = di & ((mask(8) << 23) | (mask(12) << 11));
    float d0 = fi.f;

    fi.i = di & ((mask(8) << 23) | mask(11));
    float d1 = fi.f;

    float e = exp(n) * exp(d0) * exp(d1) / exp(1.0);

    if (!isPositive) {
        e = 1 / e;
    }
    return e;
}

static float expTbl0[128 + 1];
const size_t d0MaskLen = 12;
static float expTbl1[1 << d0MaskLen];
#define NUM_OF_ARRAY(x) (sizeof(x) / sizeof(*x))
void InitExp()
{
    for (int i = 0; i < NUM_OF_ARRAY(expTbl0); i++) {
        expTbl0[i] = exp((float)i - 1);
    }
    for (int i = 0; i < (1U << d0MaskLen); i++) {
        fi fi;
        fi.i = (127 << 23) | (i << 11);
        expTbl1[i] = exp(fi.f);
    }
}

    int d0i = (di >> 11) & mask(12);

    float e = expTbl0[n + 1] * expTbl1[d0i] * d1;

// VC9(32bit)
max=4.768372e-007, ave=8.847215e-008
sum=19910788.588494, org exp time=2.000000
sum=19910788.038545, my  exp time=1.453000

// gcc 4.3.4(32bit)
max=4.768372e-07, ave=9.258090e-08
sum=19910788.588641, org exp time=3.150000
sum=19910788.034149, my  exp time=1.490000

• xをどううまく分割するかという話になった. 整数部を除くと精度は23bitなので8bit x 3のテーブルで処理できそうと考える. x = n(整数) + (2進数8桁) + (2進数8桁) + (2進数7桁).
• ただ小数部を0～1の範囲にすると指数部も変わるためテーブル引きできない.
• 指数部を固定し, 仮数部のみの処理にしたいので小数部の正規化を考える. そのために0～1の範囲ではなく1～2の範囲にしよう. つまり x = n + d(n∈Z, d∈[1, 2)) とする.
• そうすると, dを複数の小数の和に分解することが, 仮数部を適当な固定幅で切った論理和に置き換えられる.

f = f1 | f2 | f3
f1 = f & 0b11111111000000000000000
f2 = f & 0b00000000111111110000000
f3 = f & 0b00000000000000011111111

• e^-x = 1 / e^xだからx >= 0とする.
• e^x = Σxⁱ / i!を使うのだけど収束に時間がかかるのでx = n + d(n∈Z, d∈[0, 1))と分解しよう. e^x = eⁿ e^d. n～128程度でfloatの最大値になるからそちらはテーブル引きすればいい.
• その場合11!=39916800～4.0e7だから1e-7の精度を出すにはΣはi = 11ぐらいまでとる必要がある. これは大きすぎる. d >= 0.5の場合はd = 0.5 + d1, d1∈[0, 0.5)と更に分解してみよう. e^x = eⁿ e^0.5 e^d1
• この場合0.5⁸/8!～1e-6ぐらい. うーん, まだ大きい.d1 = 0.25 + d2, d2∈[0, 0.25)ともう一つ分けるか. それでも0.25⁶/6!～3e-7. いまいち. 別の方法を考える.
• d1, d2は要は2進分数なんだからそちらもテーブル引きすればよさそう.
• x = n + (2進数8桁のd1) + (次の2進数8桁のd2) + ...として全部テーブル引きでいいか. なら級数展開は不要だなあ.

  ptime from_time_t(std::time_t t)
  {
    ptime start(gregorian::date(1970,1,1));
    return start + seconds(static_cast<long>(t));
  }