log1pとかexpm1がなければ自分で作れば良いじゃない

マリー=アントワネット様の素晴らしいお言葉に勇気付けられながら，IA-32 + MSVC++でのlog1pの1実装．
以下の実装はIA-32 + MSVC++な環境でのみ有効なコードです．

// log(1.0 + x)を計算する．特に|x| << 1.0で高精度であることが期待できる．
double log1p(double x)
{
  if(-1.0 + 1.4142135623730950488016887242097 / 2.0 < x
       && x < 1.0 - 1.4142135623730950488016887242097 / 2.0)
  {
    __asm{
      fldln2;
      fld qword ptr [x];
      fyl2xp1;
      fstp qword ptr [x];
    }
  }
  else{
    x = std::log(1.0 + x);
  }

  return x;
}

GCCではでlog1pが宣言されているっぽいので恐らく不要．実際に試してはいないけれど．もちろんにはlog1pがimportされていないのでC++で使う場合にはC(C99)のlog1pを直接引っ張ってくる必要がある．
expm1もほぼ同様．とゆ〜か，昨日貼ったexp(x)のアセンブリで最後に1足している部分を削るだけ．また，|x| < 1.0が対象なのでfrndintで(-1.0, 1.0)に縮小して，後でfscaleで元の大きさに拡大する変形（ $e^x = e^{x - [ x ]}\cdot e^{[ x ]}$ な変形）も不要になる．
以下の実装はIA-32 + MSVC++な環境でのみ有効なコードです．

// exp(x) - 1.0 を計算する．特に|x| << 1.0において高精度であることが期待できる．
double expm1(double x)
{
  if(-1.0 < x && x < 1.0){
    __asm{
      fld qword ptr [x];
      fldl2e;
      fmulp st(1), st;
      f2xm1;
      fstp qword ptr [x];
    }
  }
  else{
    x = std::exp(x) - 1.0;
  }

  return x;
}

log1pと同様，これもGCCではにexpm1があるようなので多分いらない．