So-net無料ブログ作成
検索選択

SML/NJの浮動小数点処理のバグ(だと思う) [SML]

円周率は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だけなのだ。
nice!(0)  コメント(2)  トラックバック(0) 

nice! 0

コメント 2

mmsaito

どこで最後の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) 

ether

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

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

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

トラックバック 0

この記事のトラックバックURL:
※言及リンクのないトラックバックは受信されません。

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。