図1: 浮動小数点数の表現形式と x86 の丸め精度
Java言語は strictfp というメソッド、クラスおよ びインタフェース修飾子を用意している。プログラマはこれを 用いて、そのメソッドやクラス中の浮動小数点演算が、どんな Java仮想マシンでもまったく同じ結果を導くこと、つまり bitwise exact reproducibility を要求できる。x86 プロセッ サでは、strictfp には Java仮想マシンのサポート が必要である。しかし、Sun Microsystems社 (Sun) の最新の Java仮想マシン (JDK 1.3beta for Win32/x86) ですら、 strictfp を実装していない。
筆者は (おそらく) 世界で初めて Javaの実行系に strictfp を実装した。Java Grande Forum の Numerics Working Group が推奨するオーバヘッドの小さい手 法を、JIT コンパイラ shuJIT に実装し、手法が有効であるこ とを実証し、実行効率を評価した。
本稿ではまず、Java に strictfp が導入された背景、 strictfp の意味を述べる。続いて、手法、実装、そ の評価について報告する。
Java 言語、仮想マシン (以下まとめて単に Java) に strictfp が導入された背景、また、筆者が実装した 背景を述べる。
浮動小数点演算について、(おそらく) あらゆる x86 用 Java 仮想マシンは言語仕様を満たしていない。
Java 言語 [1] および Java仮想マシン [2] は浮動小数点演算規格 IEEE 754 [3]を採用している。これは、単精度、 倍精度数の表現形式、加減乗除、剰余、平方根、比較、整数と の変換などを規定している。例えば Sun Microsystems 社 (Sun) の SPARC と Intel の x86 は、どちらもこの IEEE 規格に 準拠したプロセッサアーキテクチャである。しかし、SPARC と x86 は、ある同じ演算に対して異なる結果を返す。これは exact reproducibility を標榜していた Java にとって問題で あった。
SPARC と x86 のこの差異を、strictfp 導入前の Java言語仕様および仮想マシン仕様 (以下、旧仕様) はまった く許していなかった。つまり、どんな環境でも、同じ演算に対 してはまったく同じ結果が、換言すれば exact reproducibility が、求められていた。Java の仕様と合致し ていたのは SPARC の挙動の方である。x86 用仮想マシンはこの点において Javaの仕様に従っていなかった ことになる。
x86 上で、Java の旧仕様に沿った、SPARC と同じ挙動を実現 しようとすると、どうしてもプロセッサの演算性能を最大限引 き出すことはかなわない。x86 では Java プログラムの性能が 悪いということになると、Intel などの x86 メーカや、HPC (High Performance Computing: 高性能計算) に Java を利用 したい人々は困る。Intel や IBM、SuperCede 社は x86 や PowerPC の性能を Java で引き出せるように Sun に要求し、 それに対して Sun は Java の拡張案 [4] を提示し、 public review を求めた (1998年 5月)。
Java Grande Forum (http://www.javagrande.org/) の Numerics Working Group (http://math.nist.gov/javanumerics/) (以下 JGNWG) は、Sun の案は Java の predictability を損なう、として、対案 [5] を提出、発表した (1998年 9, 10月)。 この JGNWG 案は、Java を数値計算に応用する際の諸問題 [6] の一部ではあるが、複素数クラス、 効率的にアクセスできる行列クラスの新設、light weight ク ラスなど、特に要求の多い問題をカバーするものであった。
Sun は JGNWG 案のうち strictfp に関する部分を、 若干修整したのち、Java 言語および仮想マシンの仕様に採り 入れた[7] [8]。JGNWG 案の strictfp に 関する部分は、Java仮想マシンに新しい命令を導入しなければ 実現できない内容を含んでいて、確かに Sun としてはそのま までは受け入れ難かったのだろう。Sun が採用した Java への 修整は、x86 用 Java仮想マシンの演算性能を妨げないことの みを目的とした、最低限のものであった。PowerPC の fused mac (multiply-accumulate) 命令の許容や、複素数クラス、 light weight クラスなどの導入は見送られた。
strictfp はメソッ ド、インタフェースまたはクラスの修 飾子であり、その文脈では、どの JVM も同じ浮動小数点演算 に対してまったく同じ結果を返すべきである。x86 プ ロセッサの挙動は strictfp が要求するものと異な るので、x86 用 JVM、JIT コンパイラ は strictfp のセマンティクスを実装する必要があ る。
x86 の挙動が、SPARC や strictfp の要求する挙動 とどのように異なるのかを説明する。
IEEE 754 の正規化数 (normalized number) は次の形式で表される。
(-1)s 2E (1.b1b2…bp-1)ここで各記号の意味は次の通りである。
また、単精度、倍精度など各形式のパラメタは次の通りである。
s : 符号ビット (0 か 1) bi : 仮数部の 1ビット (0 か 1) p : 仮数部の精度 E : 指数部
x86 プロセッサ内部では、浮動小数点 数は常に、単精度であろうと倍精度であろうと、拡張精度と同じ 80 bit で保持される 。ただし、x86 に対して丸め精度 (rounding precision) を単精度、倍精度、または拡張精度に設定するこ とが可能である。浮動小数点数は、レジスタへ格納される際に、 仮数部がその設定に応じて丸められる。ここで注意すべきは、 x86 では丸め精度の設定が仮数部にし か影響を与えないということである。指数部は丸め精 度の影響を受けず、常に 15 bit の精度を持つ(図1)。これが x86 の仕様であり、 x86 用の最初の浮動小数点演算ユニット (FPU) から現在最新 の x86、また、Intel 社以外が設計した x86 互換プロセッサ に至るまで、この仕様を受け継いでいる。
形式長(ビット) p : 仮数部の精度 Eのビット数 Eの最大値 Eの最小値 単精度 32 24 8 +127 -126 倍精度 64 53 11 +1023 -1022 x86の拡張精度 80 64 15 +16383 -16382
strictfp のセマンティクスは、指数部も仮数部と同 様に、各表現形式に沿った精度、つまり単精度なら 8ビット、 倍精度なら 11ビットであることを前提としている。SPARC の 仕様はそのようになっている。x86 では、常に指数部の精度が 15ビットであるため、ある計算の結果が strictfp の要求や SPARC とは異なる。
x86 で strictfp を実現しようとした場合、指数部 が overflow、underflow を起こす境界が問題となる。単精度 なら 8ビット(-126〜+127)、倍精度なら 11ビット(-1022〜 +1023)で overflow、underflow が起きなければならないとこ ろ、15ビットで表現できる範囲(-16382〜+16383) まで起きない。
例えば次の計算では、途中で overflow が起きて +∞ が得ら れるべきであるが、x86 では overflow が起きずに 2+1023 が得られる。 この問題は、加減乗除で発生し得る。
例 1)
register = 2+1023
register = register + 2+1023
register = register - 2+1023
例 2)
register = 2+1023
register = register × 2.0
register = register × 0.5
現在、strictfp のセマンティ クスを実装した x86 用 Java仮想マシンは (おそらく) ひとつもない。Sun の最新の Java Development Kit (JDK) (JDK 1.3beta) も然りである。
Java の仕様に strictfp が採り入れられたことに応 じて、Java Development Kit (JDK) 1.2 (Java 2 SDK 1.2) か らは、Javaコンパイラ javac が修飾子 strictfp を扱えるようになった。しかしこれは単に、 javac がソースコード中の strictfp 修飾 子を認識して、クラスファイル中のアクセスフラグに ACC_STRICT として反映させることができるようになっ たに過ぎない。サポートされたのは strictfp のシ ンタクスのみであり、現行 (JDK 1.3beta 以前) の Java仮想 マシンも strictfp のセマンティクスをサポートし ていない。
さらに言えば、javac も strictfp のセマ ンティクスをサポートする必要があるが、現行の javac は対応していない (JDK 1.2 pre-release 2 for Linux で JDK 1.3beta の javac を用いて確認 した)。javac はコンパイル時に簡単な定数式なら値 を計算してしまうので、その計算の際に、定数式の文脈 --- strictfp かそうでないか --- を考慮する必要があ る。幸い、strictfp のスコープは静的に決まるので、 コンパイル時に定数式を計算してしまうことは原理的に可能で ある。いずれにせよ、javac は Java仮想マシン上で 動作するので、まず Java仮想マシンが strictfp を サポートする必要がある。
1999年 2月、IEEE 754 が浮動小数点数についてどこまでを規 定しているかを知らなかった筆者は、はたして IEEE 規格が演 算結果や丸めの方法まで規定しているものかどうかを調べる目 的で、SPARC と x86 で挙動が異なる計算を探した。もし挙動 が異なるとしたら、丸めについてだろう、と見当をつけて様々 な値について乗算、除算を試し、挙動が異なる計算を見つけた。 これは Java で問題になるはず、と考えつつ、スクリプト言語 Ruby のメイリングリストで報告したきり、放置していた。
7月、ふと SPARC と x86 の挙動の違いを思いだし、NetNews の fj.comp.arch および fj.comp.lang.java に「SPARC と x86 は両方とも IEEE規格に準拠しているのだろうか？それとも…」と投稿した。 この投稿に対して、x86 のアーキテクチャと Java についてさ まざまな議論、助言、調査があり、それらを元に、x86 の挙動 (仮数部が常に 15ビット) を理解するに至った。
Java に関する議論のなかで strictfp が話題となり、 8月、高木さんより「shuJITは世界初完全strictfpに なるか? :-)」という suggest :-) を頂き、shuJIT に strictfp のセマンティクスを実装することを決めた。 strictfp を実装するということは Java仮想マシン の挙動を変えることになる。そのため、Java インタプリタ、 JIT コンパイラ、Ahead-of-Time コンパイラなどを実装、改造 する必要がある。筆者が、自らの研究、実験のために JIT コンパイラ shuJIT [11] [10] を開発してきたことは幸いであった。
本応募は次の点について、新しいかもしくは有用である。
strictfp のセマンティクスを実現する手法として、 現在知られているおそらく最良のものは Golliver が提案した手法 [12] であり、JGNWG も、その手法の実装を x86 用 Java実行系の開発者に対して推奨 している [5]。Golliver の手法は文献 [5] 中でも述べられ ている。
x86 の問題は、FPU レジスタ上では常に指数部の精度が 15ビッ ト(-16382〜+16383)あることであり、加減乗除で発現する。 strictfp を実現するためには、単精度なら 8ビット (-126〜+127)、倍精度 なら 11ビット(-1022〜+1023)で表現できる範囲で overflow、underflow を起こさせる 必要がある。そのために次の手法を用いる。
<- 53(or 24)bit ->
<- 53(or 24)bit ->
0.000...0001xxx...xxx 1000...000 1
または 0xxx(1つ以上の 1)xxx...
× 2**(-1022 (or -126))
(`or' は単精度の場合)
そこで、レジスタ上で適切に underflow を起 こして非正規化数にしてしまうことで2度丸めを防ぐために、次の手法を用 いる。
まとめると、加減乗除算で store-reload を行い、乗算、除算で scale down and up を行うことで、 strictfp を実現できる。
筆者は JIT コンパイラ shuJIT [11] [10] に Golliverの手法 [5] [12] を実装 した。この節で、どのような実装を行ったのかを述べ る。
JIT コンパイラに strictfp を実装するということ は、strictfp の文脈 (メソッド、クラス) における 浮動小数点演算に対して、strictfp の要求を満たす ネイティブコードを生成する、ということである。JIT コンパ イラが生成するコードは次の処理を行う必要がある。
また、次の付加機能も実装した。
丸め精度は、単純に考えれば、倍精度数の演算の前には倍精度 に、単精度数の演算の前には単精度に設定すればよい。しかし、 これをそのまま実装すると、倍(単)精度数の演算のあと単(倍) 精度数の演算があった場合、丸め精度を単(倍)精度に設定し直 すことになる。最悪の場合、浮動小数点演算のたびに、丸め精 度を設定することになってしまう。丸め精度の設定はメモリ (キャッシュ)アクセスを要し、決して軽い処理ではない。
しかし実は、単精度数の演算の際も、 丸め精度は倍精度に設定しておけば充分である。つま り strictfp の文脈では丸め精度は倍精度に設定し ておけばよく、演算ごとに単精度、倍精度と設定しなおす必要 はない。
丸め精度を倍精度にしたまま単精度数の (strictfp) 演算を行うと、仮数部は次のタイミングで 2度丸められる。
乗算の場合、単精度数は 2進 24桁 (24ビット精度) なので、 単精度数どうしの積はたかだか 2進 48桁である (*1)。ゆえに、上記 2段階丸めの 1度目の (53ビットへの) 丸めの影響を受けない。
(*1) 2進 n桁どうしの積はたかだか 2n桁で ある。[証明] 2進 n桁の最大の整数である 2n - 1 (10進) の 2乗を考える。 (2n - 1)2 = 22n + 1 - 2n+1、これは 2進で 2n桁である。
除算の場合を考える。2度丸めの結果が 1回で丸められた場合 と異なってしまうのは、真の(丸められる前の)商の仮数部が次 の値(2進)である場合に限られる。
正規化数の場合:この仮数部は、1回で単精度(24ビット)に丸められ場合は切り上げ となるが、倍精度(53ビット)に丸められた後で単精度に丸めら れると切り捨てとなってしまう。<- 53bit -> <- 24bit -> 1.xxx.....xxx 1000...000 1 または 0xxx(1つ以上の 1)xxx...非正規化数の場合 (scale down and up によって 1回目倍精度に丸められた時点で非正規化数となる):<- 53bit -> <- 24bit -> 0.0...01x...x 1000...000 1 または 0xxx(1つ以上の 1)xxx... × 2**(-126)
ところが単精度数の除算でこの形式の商は得られない。もし得 られるとしたら次の関係が成り立つはずである。
除数(24ビット) × 商 = 被除数(24ビット)被除数は 24ビット精度なので、次の形式で表現できなければならない。 24ビット精度の除数に上記の商を乗じて、下記の形式の値を得ることはできない。
<- 24bit ->
1.xxx.....xxx 000...
または 111... …(1)
この記述では厳密な証明にはなっていないが、24ビットの除数 に上記の商を筆算で乗じて (1) の被除数を得られるか考える と…どうにも得られそうにない。上記の商の 25ビット目の `1' とその次の 54ビット目の `1' が 29ビット離れているためである。
丸め精度は倍精度に設定しさえすれば充分そうである。Java仮 想マシンが元から倍精度に設定していれば、あらためて設定し 直す必要もない。しかしすべての Java仮想マシンがはじめか ら丸め精度を倍精度に設定しているとは限らない。例えば、 (Blockdown による) Linux用 JDK 1.1.7 は Linux の初期設定 である拡張精度のままにしている。
このような場合、JITコンパイラの初期化時に倍精度に設定し てしまう方法も考えられるが、JITコンパイラの有無によって プログラムの実行結果に差が出かねない。それを避けるため今 回は、strictfp メソッドの 先頭に、倍精度に設定するコードを、末尾に元の丸め精度に戻 すコードを生成するようにした。
ただ、Java仮想マシンによる設定が元から倍精度であった場合、 再設定するのも無駄である。JIT初期化時に Java仮想マシンの 設定を調べて、それが倍精度でない場合に限って丸め精度を設 定するコードを生成するようにした。
前で述べたように、store-reload で正しく overflow, underflow を起こさせ、scale down and up によってレジスタ上で正しく underflow を起こして 2度丸 めを防ぐ。store-reload は加減乗除算で必要であり、scale down and up は乗算、除算で必要である。
今回 strictfp を実装した JIT コンパイラ shuJIT では、store-reload を (残念ながら) 元から常に行っている。 すなわち、浮動少数点数の加減乗除算では、結果は必ずメモリ へストアされる。今回 store-reload のために手を加える必要 はなかった。
scale down and up は新たに実装した。scale down and up は 2の巾数によるスケーリングなので、x86 では次の 2通りの実 装が考えられる。双方を実装して性能を比較した。
さらに、レジスタに置いておく場合でも、いくつかの選択肢がある。
今回、上記 3 は実装せず、レジスタに置いておく scale は、 スケーリングを行う方法によって次のようにした。
スケーリングを乗算で行う場合は 2n という形式 の scale が必要なために逆数 2-n の計算コスト が大きいので、4つの scale すべてをレジスタに置いておく。 それに対し、スケーリングに fscale命令を用いる場 合、scale は巾 n で表現できるので、逆数 -n の計算コスト (fchs命令) は小さい。必要に応じて -n を算出して もレジスタからロードしても計算時間にほとんど差がなかった ので、-n を算出する方法を選んだ。
乗算で行う場合: 4つの値すべて fscale命令で行う場合: 2つの値 (-(16383 - 1023) と -(16383 - 127))
strictfp の実装によって、どの程度のオーバヘッド が導入されたのかを調べた。実験環境は、Pentium with MMX technology / 233MHz, Linux 2.2.13pre7, JDK 1.1.7 version 3 (green thread) である。コンパイルは Linux 用 JDK 1.2 pre-relase 2 付属の javac コマンドで行った。
10回の乗算 (除算) を含んだ 式を 106 回計算 するのにかかった時間を計った (StrictfpBenchmark.java)。 ある程度の最適化を行う C コンパイラであれば乗算 (除算) をループの外へ出してしまうが、javac はそれを行 わないので、このコードでも乗算 (除算) のベンチマークにな る。 scale down and up のための scale の FPU レジスタへの事前 ロードの効果を知るため、事前ロードを行う版と、scale をレ ジスタに置かずに計算のたびにメモリ(キャッシュ)から読む版 とを比較した。単位はミリ秒である。
strictfp のオーバヘッドと事前ロードの効果が見て とれる。ただしここで strictfp 有無の差として表 れているオーバヘッドは scale down and up のものだけであ る。前で述べたように、shuJIT ではstrictfp の有無に関わらず、常に store-reload を行ってしまっているからである。 除算より乗算に時間がかかっている理由は不明である。
(ミリ秒) 乗算 除算 strictfp なし 4004 903 fscale でスケーリング, 事前ロードあり 5648 2749 乗算でスケーリング, 事前ロードあり 10434 1117 乗算でスケーリング, 事前ロードなし 11712 1976
JGNWG は Golliverの手法による性能低下は 2〜4倍程度と予想 している [5]。今回の実験では、 strictfp ありの場合、実行時間は 1.24倍〜 2.93倍 となっている。ただし今回計測されたペナルティは store-reload のものを含んでいない。store-reload も含める とどの程度のペナルティとなるのかはまだわかっていない。賢 い JIT コンパイラが生成するであろうコードを手で記述して ペナルティを見積もることは有意義であろう。
Java言語仕様に strictfp が導入され [7]、x86 における浮動少数点演算 のセマンティクス問題は解決したかというと…実装が追い付い ていないことはもとより、仕様の妥当性の問題も残されている。
Java言語仕様は、x86 用 Java仮想マシンの演算性能を妨げな い目的で改定された [7]。とこ ろが、この目的が充分に達成されたとはとても言えない。 SPARC や MIPS は、単精度数、倍精度数それぞれに対応した演 算命令 (例: fmuls と fmuld、 mul.s と mul.d) を持っているが、x86 は 持っていない。浮動小数点演算は常に 80ビットの内部実数型 どうしで行われ、丸め精度の設定に基づいて仮数部が丸められ る。
更新された言語仕様 [7] は、 strictfp の文脈でない限りは、ある範囲で単精度数、 倍精度数が (8, 11 ではなく) 15ビットの指数部を持つことを 許している。しかし飽くまで、許されたのは指数部だけであり、 仮数部は単精度なら 24ビット、倍精度なら 53ビットの精度で ある必要がある。単精度、倍精度用それぞれ用の演算命令を持 たない x86 でこの仕様に従うためには、次のどちらかを行う 必要がある。例えば単に、単精度数を倍精度数として扱ってし まったのでは、言語仕様に違反してしまう。
strictfp と、15ビットの指数を持つ単精度数 (float-extended-exponent)、倍精度数 (double-extended-exponent) は、Java 2 から導入された。こ のことは、言語仕様の改定版 [7] の題名「Updates to the Java Language Specification for JDK Release 1.2 Floating Point」からも判る。
では、JDK 1.1 や 1.0 はどうあるべ きだったのか。これが明確に されていない。strictfp 導入以前、Sun が Java の `exact reproducibility' をうたっていたということ は、常に strictfp の文脈だと考えるのが自然では ある。しかし、以前の言語仕様は「IEEE 754 に準拠」の一点 張りである。SPARC とは異なる x86 の挙動も「IEEE 754 に準 拠」である。x86 も言語仕様に沿っていたという解釈ができる。
JDK 1.1 の Java 実行系の開発者は困る。言語仕様の曖昧だっ た点が改定によって明確になった、と考えるなら、JDK 1.1 に も 15ビット指数部が許されていたことになる。しかし、 strictfp 導入の改定前後で、異なる 2つの言語仕様 があり、15ビット指数部は改定によってはじめて許されたのか もしれない。
前で述べたように、Java仮想マ シンが strictfp に対応した上で、Javaコンパイラ も strictfp に対応する必要がある。Javaコンパイ ラは、定数式を計算する際に、その式の文脈が strictfp かそうでないかを考慮しなければならない。
javac は JDK 1.2 から strictfp のシン タクスには対応した。しかし、上で述べたようなセマンティクスへの対応はしていない 。strictfp に対応させた Java仮想マシン (Linux 用 JDK 1.2 pre-relase 2) 上で JDK 1.3beta の javac (tools.jar 中のもの) を用いて確 認した。
strictfp の文脈かそうでないかで結果が異なる計算 の例を示す (StrictfpTest.java)。こ こで挙げる例は特に、2度丸め (この 節を参照) の影響を受けてしまう計算なので、2度丸めに 対応していないと strictfp の要求する結果になら ない。
`strictfp: ...' の行が strictfp の文脈 での計算結果、`default: ...' の行がデフォルトの 文脈での計算結果である。このテストは、他の x86 用 Java仮 想マシンでは `strictfp:' の結果が `default:' と同じ値になってしまう。
% java StrictfpTest shuJIT for Sun JVM/IA-32 Copyright 1998,1999 by SHUDO Kazuyuki 1.112808544714844E-308 (0x0008008000000000) * 1.0000000000000002 (0x3ff0000000000001) default : 1.112808544714844E-308 (0x8008000000000) strictfp: 1.112808544714845E-308 (0x8008000000001) 2.225073858507201E-308 (0x000fffffffffffff) / 0.9999999999999999 (0x3fefffffffffffff) default : 2.2250738585072014E-308 (0x10000000000000) strictfp: 2.225073858507201E-308 (0xfffffffffffff)
本稿で述べた内容は次の通りである。
問題、言語仕様改定の背景、Golliver の手法の解釈、そして、 実装した際の工夫を述べ、実装を評価した。また、残されてい る問題について言及した内容は、Java Grande Forum のメンバ や、Sun で言語仕様にたずさわっている者にとっても有効であ ろう。
本稿が、問題の理解、x86 用 Java実行系の実装、よりよい言 語仕様のために役立てば幸いである。
NetNews での議論、文献の紹介などでお世話になった伊藤さん <eiki ＠ prd.cs.fujitsu.co.jp>、片山さん <kate ＠ pfu.co.jp>、高木さん <takagi ＠ etl.go.jp>、 前田さん <maeda ＠ is.uec.ac.jp> (50音順) に感謝し ます。また、なかなか論文を書かない筆者を暖かく指導して下 さる村岡先生に感謝します。 あらゆる面から私を支え、今日も午前様となるあいかわらず遅 い帰宅にもじっと耐えてくれた妻の麻里に深く感謝します。