SML/NJの浮動小数点処理のバグ（だと思う）

円周率はIEEE754の64ビット浮動小数点数で表すと「40 09 21 fb 54 44 2d 18」というバイト列に相当する。
SMLで、PackRealストラクチャを使わずに、この最後の2バイト「2d 18」を取り出したい。そこで次のようなコードを書く。

fun main () =
  let
    fun println s = print (s ^ "\n")
    val pi' = #man (Real.toManExp Math.pi) * Math.pow (2.0, 53.0)
    val a = Math.pow (2.0, 16.0)
    val b = Real.rem (pi', a)
    val c = Real.toInt IEEEReal.TO_NEAREST b
  in
    println ("Real.precision=" ^ Int.toString Real.precision);
    println ("c=" ^ Int.toString c)
  end

val _ = main ()

cは11544=0x2d18が入ることが期待されている。さまざまなSML実装で試した結果は次の通り。環境はすべてx86 Linuxである。

SML/NJ

- use "divreal.sml";
[opening divreal.sml]
[autoloading]
[library $SMLNJ-BASIS/basis.cm is stable]
[autoloading done]
Real.precision=53
c=11541
val main = fn : unit -> unit
val it = () : unit

Poly/ML

> use "divreal.sml";
Real.precision=53
c=11544
val main = fn : unit -> unit
val it = () : unit

SML#

# use "divreal.sml";
Real.precision=53
c=11544
val main = fn : unit -> unit

MLton

Real.precision=53
c=11544

SML/NJだけ仲間はずれであり、今回の目的にもそぐわない。

バグか？というとバグだと断言できるほど浮動小数点計算の仕様に通じていないのだが、素人考えとしては、ここで行っている計算は2の階乗の掛け算と2の階乗による余りの算出だけなので、それぞれ指数部に対する足し算と、仮数部に対するビット演算だけでできそうだ。そういうのは誤差が出ないでくれればいいのではないかと思う。

こんなことやる場合あるのか？というと、SMLの仕様の範囲でrealのビット表現を直接得る方法はPackRealストラクチャしかない。しかしPackRealは必須ではなく、提供しているのは上記のSML実装のうちMLtonだけなのだ。

コメント 2

どこで最後の2bitが欠けてしまうのですね。ちょっと怖いですね。
PackRealは知らなかったので、バイナリライトしたいときはC言語インターフェースをつかって取り出しています。
それによると、さすがにMath.piが狂っていることはなさそうですね。

CM.make "$c/c.cm";
val a = C.new C.T.double;
val _ = C.Set.double (C.rw a, Math.pi);
val pa = C.Ptr.cast (C.T.pointer C.T.uchar) (C.Ptr.inject (C.Ptr.|&| a));
val img = Vector.tabulate(8, fn i => C.Get.uchar (C.Ptr.sub (pa,i)));
val img = #[0wx18,0wx2D,0wx44,0wx54,0wxFB,0wx21,0wx9,0wx40]
by mmsaito (2014-01-18 04:45) 

コメントありがとうございます。

C言語レベルでプローブしないと調べようがないと思ってたんですがSML/NJのCインターフェイスはいまいち資料が乏しくてわからないなあと思っていたところでした。
by ether (2014-01-19 00:23)